From 10087691a7f123182fcd41535fe0f595ac1f4b07 Mon Sep 17 00:00:00 2001
From: Alex Austregesilo <aaustreg@jlab.org>
Date: Fri, 26 Feb 2021 13:21:45 -0500
Subject: [PATCH] a cut action to select PID with beta

---
 libraries/DSelector/DCutActions.cc | 78 ++++++++++++++++++++++++++++++
 libraries/DSelector/DCutActions.h  | 21 ++++++++
 2 files changed, 99 insertions(+)

diff --git a/libraries/DSelector/DCutActions.cc b/libraries/DSelector/DCutActions.cc
index 206c859..fb0220d 100644
--- a/libraries/DSelector/DCutActions.cc
+++ b/libraries/DSelector/DCutActions.cc
@@ -252,6 +252,84 @@ bool DCutAction_PIDDeltaT::Perform_Action(void)
 	return true;
 }
 
+string DCutAction_PIDBeta::Get_ActionName(void) const
+{
+	ostringstream locStream;
+	locStream << DAnalysisAction::Get_ActionName() << "_" << dPID << "_" << dSystem << "_" << dMinBetaCut  << "_" << dMaxBetaCut;
+	return locStream.str();
+}
+
+void DCutAction_PIDBeta::Initialize(void)
+{
+	if (dFunc_BetaCut_SelectPositive == nullptr){
+		dFunc_BetaCut_SelectPositive = new TF1("dFunc_BetaCut_SelectPositive","pol0",0.0, 12.0);
+		dFunc_BetaCut_SelectPositive->SetParameter(0,dMaxBetaCut);
+	}
+	else cout << "Using user function for positive PID beta cut " << endl;
+
+	if (dFunc_BetaCut_SelectNegative == nullptr){
+		dFunc_BetaCut_SelectNegative = new TF1("dFunc_BetaCut_SelectNegative","pol0",0.0, 12.0);
+		dFunc_BetaCut_SelectNegative->SetParameter(0,dMinBetaCut);
+	}
+	else cout << "Using user function for negative PID beta cut " << endl;
+}
+
+bool DCutAction_PIDBeta::Perform_Action(void)
+{
+	for(size_t loc_i = 0; loc_i < dParticleComboWrapper->Get_NumParticleComboSteps(); ++loc_i)
+	{
+		DParticleComboStep* locComboWrapperStep = dParticleComboWrapper->Get_ParticleComboStep(loc_i);
+
+		//final particles
+		for(size_t loc_j = 0; loc_j < locComboWrapperStep->Get_NumFinalParticles(); ++loc_j)
+		{
+			DKinematicData* locKinematicData = locComboWrapperStep->Get_FinalParticle(loc_j);
+			if(locKinematicData == NULL)
+				continue; //e.g. a decaying or missing particle whose params aren't set yet
+
+			//-2 if detected, -1 if missing, > 0 if decaying (step where it is the parent)
+			int locDecayStepIndex = locComboWrapperStep->Get_DecayStepIndex(loc_j);
+			if(locDecayStepIndex != -2)
+				continue; //not measured
+
+			if((dPID != Unknown) && (locKinematicData->Get_PID() != dPID))
+				continue;
+
+			// determine detector system
+			DetectorSystem_t locSystem = SYS_NULL;
+			double locBeta_Timing = 0.0;
+			if(ParticleCharge(locKinematicData->Get_PID()) != 0)
+			{
+				const DChargedTrackHypothesis* locChargedTrackHypothesis = dynamic_cast<const DChargedTrackHypothesis*>(locKinematicData);
+				if(locChargedTrackHypothesis != NULL)
+				  {
+					locSystem = locChargedTrackHypothesis->Get_Detector_System_Timing();
+					locBeta_Timing = dUseKinFitFlag ? locChargedTrackHypothesis->Get_Beta_Timing() : locChargedTrackHypothesis->Get_Beta_Timing_Measured();
+				  }
+			}
+			else
+			{
+				const DNeutralParticleHypothesis* locNeutralParticleHypothesis = dynamic_cast<const DNeutralParticleHypothesis*>(locKinematicData);
+				if(locNeutralParticleHypothesis != NULL)
+				  {
+					locSystem = locNeutralParticleHypothesis->Get_Detector_System_Timing();
+					locBeta_Timing = dUseKinFitFlag ? locNeutralParticleHypothesis->Get_Beta_Timing() : locNeutralParticleHypothesis->Get_Beta_Timing_Measured();
+				  }
+			}
+
+			if((dSystem != SYS_NULL) && (locSystem != dSystem))
+				continue;
+
+			TLorentzVector locP4 = dUseKinFitFlag ? locKinematicData->Get_P4() : locKinematicData->Get_P4_Measured();
+			if(  locBeta_Timing > dFunc_BetaCut_SelectPositive->Eval(locP4.P()) || locBeta_Timing < dFunc_BetaCut_SelectNegative->Eval(locP4.P()))
+				return false;
+
+		} //end of particle loop
+	} //end of step loop
+
+	return true;
+}
+
 string DCutAction_NoPIDHit::Get_ActionName(void) const
 {
 	ostringstream locStream;
diff --git a/libraries/DSelector/DCutActions.h b/libraries/DSelector/DCutActions.h
index 11604ed..70382e8 100644
--- a/libraries/DSelector/DCutActions.h
+++ b/libraries/DSelector/DCutActions.h
@@ -99,6 +99,27 @@ class DCutAction_PIDDeltaT : public DAnalysisAction
 		double dTargetCenterZ;
 };
 
+class DCutAction_PIDBeta : public DAnalysisAction
+{
+	public:
+		DCutAction_PIDBeta(const DParticleCombo* locParticleComboWrapper, bool locUseKinFitFlag, double locMinBetaCut, double locMaxBetaCut, Particle_t locPID = Unknown, DetectorSystem_t locSystem = SYS_NULL, string locActionUniqueString = "") :
+			DAnalysisAction(locParticleComboWrapper, "Cut_PIDBeta", locUseKinFitFlag, locActionUniqueString),
+			dFunc_BetaCut_SelectPositive(nullptr), dFunc_BetaCut_SelectNegative(nullptr), dMinBetaCut(locMinBetaCut), dMaxBetaCut(locMaxBetaCut), dPID(locPID), dSystem(locSystem) {}
+
+		string Get_ActionName(void) const;
+		void Initialize(void);
+		void Reset_NewEvent(void){}
+		bool Perform_Action(void);
+		TF1 *dFunc_BetaCut_SelectPositive;
+                TF1 *dFunc_BetaCut_SelectNegative;
+
+	private:
+		double dMinBetaCut;
+		double dMaxBetaCut;
+		Particle_t dPID;
+		DetectorSystem_t dSystem;
+};
+
 class DCutAction_NoPIDHit : public DAnalysisAction
 {
 	public: