-
Notifications
You must be signed in to change notification settings - Fork 9
Class_structure
- The DSelector Class Structure
In the .h
header file you can find a class defined as
DSelector_program-name
that inherits from the DSelector
class. The
DSelector
class is defined in the GlueX ROOT analysis library and
provides all of the methods that are necessary to read the ROOT tree. It
inherits from the ROOT TSelector
class that provides the methods
Init(), Process()
and Finalize()
and the structure to loop over all
of the entries in the ROOT tree; these form the base of the .C
file
generated.
The DSelector_program-name
can be run in several ways. For testing
with just one input file, it can be run using ROOT interactively, eg
root mytreefile.root
.x $ROOT_ANALYSIS_HOME/scripts/Load_DSelector.C
name-of-the-root-tree->Process("DSelector_program-name.C+")
Note that you should supply either one or two plus signs at the end of
the DSelector
file name to increase the processing speed of your
analysis. If you specify +
, the DSelector will only be recompiled if
you modify either of the DSelector
files. The preferred option is
++
, which will force recompilation at every execution, which can be
help protect against errors that occur when switching between software
versions.
The DSelector_program-name
can be run over one or more files using a
root script that has the following form:
TChain chain("name-of-the-root-tree");
chain->Add("root-file-name");
gROOT->ProcessLine(".x $ROOT_ANALYSIS_HOME/scripts/Load_DSelector.C");
chain.Process("DSelector_program-name.C++");
Alternatively, the DSelector_program-name
can be run in parallel on
multiple cores with PROOF-Lite:
int NThreads = 16; // or as many as you want to use
TChain *chain = new TChain("name-of-the-root-tree");
chain.Add("root-file-name");
gROOT->ProcessLine(".x $ROOT_ANALYSIS_HOME/scripts/Load_DSelector.C");
DPROOFLiteManager *dproof = new DPROOFLiteManager();
dproof->Process_Chain(chain, "DSelector_program-name.C+", NThreads, "outputHistFileName", "outputTreeFileName");
Each piece of information is stored as a branch in the ROOT tree.
Several wrapper objects are implemented, to group the branches belonging
to the same reconstructed object (track/shower/combination of these)
together. These wrappers are provided by the DSelector and are defined
in ROOT_ANALYSIS_HOME/libraries/DSelector
. There are four basic types
of objects in the tree -- charged particle tracks, neutral particles,
combos and beam photons -- which leads to the definition of four
standard wrapper objects:
// OBJECT ARRAYS: RECONSTRUCTED
DChargedTrackHypothesis* dChargedHypoWrapper;
DNeutralParticleHypothesis* dNeutralHypoWrapper;
DBeamParticle* dBeamWrapper;
DParticleCombo* dComboWrapper;
The first three wrappers are used to access the reconstructed final
state particles, charged and neutral, and the reconstructed initial
state beam photons. The fourth wrapper represents the output of the
standard analysis library used by ReactionFilter, which combines the
initial and final state particles, with kinematic fitting by default, to
generate the full event. The resulting final state particles and
intermediate states (if defined) from the ReactionFilter are accessed
by additional wrappers specific to the reaction at hand. As an example,
an exclusive final state like
//Step 0
DParticleComboStep* dStep0Wrapper;
DBeamParticle* dComboBeamWrapper;
DChargedTrackHypothesis* dPiPlusWrapper;
DChargedTrackHypothesis* dPiMinusWrapper;
DChargedTrackHypothesis* dProtonWrapper;
//Step 1
DParticleComboStep* dStep1Wrapper;
DKinematicData* dDecayingPi0Wrapper;
DNeutralParticleHypothesis* dPhoton1Wrapper;
DNeutralParticleHypothesis* dPhoton2Wrapper;
For example, to loop over all charged particle tracks of an event that were used in forming all the available combos and perform a calculation using the value of the start counter energy-loss one might attempt something like the following:
// check if any charged track has a ST hit by looking at its dEdx value
int NChargedTracksHypos = (int)Get_NumChargedHypos();
int FoundStartCounterHits = 0;
for (int k=0; k<NChargedTracksHypos; k++){
dChargedHypoWrapper->Set_ArrayIndex(k);
float ST_dEdx = dChargedHypoWrapper->Get_dEdx_ST();
if (ST_dEdx>0.){
FoundStartCounterHits++; // NOTE THAT THIS CODE OVERCOUNTS THE NUMBER OF HITS!
}
}
The key to this loop is that one sets the pointer of the wrapper to the
right index by using the method Set_ArrayIndex(k)
. Now that the
wrapper points to the desired object one can access the desired value
from that object.
Note that this code will not do what the user desires for several
reasons. First, by default ReactionFilter will only write out the
initial and final state particles used in the combos saved to the ROOT
tree, unless the U1
option is specified. Second, each charged track is
reconstructed under several different mass hypotheses (DChargedTrackHypothesis::Get_TrackID()
function, so one
could in principle keep track of the number of unique track IDs, but in
general this is an issue to be very careful about.
This function is used to initialize histograms and any other class
variables which will be used by the analysis code for every event. This
method can and most likely will be called more than once, in particular
when using a TChain
with more than one root file or when using PROOF
in the context of multi-threading. This method is inherited from the
ROOT TSelector
class and overridden in the DSelector
class. The code
snippet in the Init()
method shown below will ensure that the code
that follows after will be executed only once:
bool locInitializedPriorFlag = dInitializedFlag; //save whether have been initialized previously
DSelector::Init(locTree); //This must be called to initialize wrappers for each new TTree
//gDirectory now points to the output file with name dOutputFileName (if any)
if(locInitializedPriorFlag)
return; //have already created histograms, etc. below: exit
Therefore any initialization of variables or histogram definitions intended for the whole analysis need to be done after this part of the code. This includes any definitions of Analysis Actions and Cut Actions as well as custom branches for an output ROOT tree.
The call to the method Get_ComboWrappers()
initializes all instances
of wrappers available for the analysis of the ROOT tree. This call will
configure a wrapper for each particle in the reaction. For example, if
the reaction is
//Step 0
dComboBeamWrapper = static_cast<DBeamParticle*>(dStep0Wrapper->Get_InitialParticle());
dPiPlusWrapper = static_cast<DChargedTrackHypothesis*>(dStep0Wrapper->Get_FinalParticle(1));
dPiMinusWrapper = static_cast<DChargedTrackHypothesis*>(dStep0Wrapper->Get_FinalParticle(2));
dProtonWrapper = static_cast<DChargedTrackHypothesis*>(dStep0Wrapper->Get_FinalParticle(3));
//Step 1
dStep1Wrapper = dComboWrapper->Get_ParticleComboStep(1);
dDecayingPi0Wrapper = dStep1Wrapper->Get_InitialParticle();
dPhoton1Wrapper = static_cast<DNeutralParticleHypothesis*>(dStep1Wrapper->Get_FinalParticle(0));
dPhoton2Wrapper = static_cast<DNeutralParticleHypothesis*>(dStep1Wrapper->Get_FinalParticle(1));
These wrapper objects, as shown in the example above, are instances of classes and are pointers to the corresponding instances of the particles for a given combo and information about the combo itself. All instances are defined in the DSelector library and provide methods to access any data within the tree for a given combo. So for example with the line
TLorentzVector xv4 = dPhoton1Wrapper->Get_X4_Shower();
you will access the Lorentz 4-vector for the calorimeter shower
associated with the first photon used in this example reaction for a
given combo within the combo loop. This will be discussed in more detail
below when describing the Process()
method.
The full particle combination corresponding to the reconstructed
reaction is described by a DParticleCombo
object, which contains one
or more DParticleComboStep
objects. Each DParticleComboStep
corresponds to an intermediate decay step in the overall reaction. For
example, the reaction Step0
), while
Step0
being the primary photoproduction reaction
and Step1
being the
Each final state particle is either charged or neutral and is accessed
through a wrapper of type
D
NeutralParticleHypothesis or D
ChargedTrackHypothesis. There are
many methods associated with these two classes that provide access to
all the relevant variables in the tree associated with these particles
within a combo. Intermediate particles that decay into some final state
particles will also be represented by a wrapper class of type
DKinematicData
but only if the mass of this intermediate state
particle is constrained in the kinematic fit when requested by the
reaction filter. If the intermediate particle mass is not constrained in
the kinematic fit, then its properties (e.g. 4-momentum) can be
constructed by the user based on the other reconstructed particles. All
beam photons are represented by the wrapper class DBeamParticle
. Note
that the beam photon energies are not altered by the kinematic fitter.
More details about these wrappers and their methods will be discussed
further below when discussing the loop over all combos of an event.
Histograms can be defined easily as in any root script. The user will
need to define histograms as needed for their analysis. Note that you
should define the histogram variables in the class definition in the
.h
file before creating them in the DSelector::Init()
function.
Analysis actions are a useful system to compartmentalize analysis code,
which is applied to every combination in the tree. This code can either
include histograms
(see Histogram Actions) or event selection cuts
(see Cut Actions). Every analysis action is registered
as an entry in the vector dAnalysisActions
, and they are executed in
the same order as specified in the Init()
method.
dAnalysisActions.push_back(new DHistogramAction_NAME(...));
//OR
dAnalysisActions.push_back(new DCutAction_NAME(...));
Three histograms are automatically generated and show the effect of each action on the data. Evidently, histogram actions will not alter the number of surviving events or combos.
-
NumEventsSurvivedAction
Shows the number of events that passed each analysis action. -
NumCombosSurvivedAction
Shows the number of combinations that passed each analysis action. This is a 2-dimensional histogram, as the number of initial combinations is different for each event. -
NumCombosSurvivedAction
Shows the total number of combinations that passed each analysis action, added for all events.
In order for the analysis actions to work as intended, the entire vector
of analysis actions has to be initialized in the Init()
method:
Initialize_Actions();
and reset for every event in the Process()
method:
Reset_Actions_NewEvent();
Within the loop over all combinations in the tree, the actions can be executed at any time with this command:
if(!Execute_Actions())
continue;
This will skip the rest of your DSelector in case the combination does
not satisfy all cuts, which may speed up the process. The IsComboCut
flag is automatically set for these combinations.
Several analysis actions are provided in the code and are summarized below. The reader should be encouraged to extend this list.
Explanation on how to write a completely new custom action, maybe a template generating script?
Currently available histogram actions listed in the following are defined in libraries/DSelector/DHistogramActions.cc.
Boolean variables always specify if the measured (false) or the fitted
four-vectors (true) are used. The optional string
is a unique
identifier, which is only needed if the same action is used multiple
times.
-
AnalyzeCutActions
see Analyze Cut Actions -
ParticleComboKinematics
Creates a set of histograms for the kinematic distributions of each particle in the reaction.DHistogramAction_ParticleComboKinematics(dComboWrapper, bool, "")
-
ParticleID
Creates a set of PID histograms for each particle in the reaction.DHistogramAction_ParticleID(dComboWrapper, bool, "")
-
InvariantMass
Plots the invariant mass of a particle or set of particles. 2D histograms as a function of the kinematic fit confidence level will be generated in linear and logarithmic scale. The step specifies the level in the decay chain, where 0 is the main reaction. In the example below, the mass is plotted in 100 bins between 1.0 and 2.0 GeV.DHistogramAction_InvariantMass(dComboWrapper, bool, Lambda, 100, 1.0, 2.0, "") //OR DHistogramAction_InvariantMass(dComboWrapper, bool, step, {Proton,Pi0}, 100, 1.0, 2.0, "")
-
MissingMass
This action can plot the total missing mass of the reaction or the recoil mass against a particle or set of particles. In the example below, the missing mass is plotted in 100 bins between -1.0 and 1.0 GeV.DHistogramAction_MissingMass(dComboWrapper, bool, 100, -1.0, 1.0, "") //OR DHistogramAction_MissingMass(dComboWrapper, bool, step, Proton, 100, -1.0, 1.0, "") //OR DHistogramAction_MissingMass(dComboWrapper, bool, step, {Proton,Pi0}, 100, -1.0, 1.0, "")
-
MissingMassSquared
This plots the square of the missing mass. The syntax is the same as above. -
MissingP4
This action plots missing energy, 3 missing momentum components$p_x, p_y, p_z$ and missing transverse momentum$p_T$ .DHistogramAction_MissingP4(dComboWrapper, bool, "")
-
2DInvariantMass
This action can be used to plot the invariant mass of one particle combination against another one. All particles have to exist in the same reaction step.DHistogramAction_2DInvariantMass(dComboWrapper, bool, step, {Pi0,Proton}, {Eta,Proton}, 100, 0, 2, 100, 0, 2, "")
-
Dalitz
This plots the square of two invariant mass combinations against each other. The syntax is exactly the same as above. -
vanHove
This action plots the vanHove plot for three particles or three sets of particles.DHistogramAction_vanHove(dComboWrapper, bool, Proton, Eta, Pi0, "") //OR DHistogramAction_vanHove(dComboWrapper, bool, {Proton,PiMinus}, Eta, Pi0, "")
-
vanHoveFour
Plots the great-circle plot for 4 particles or particle systems. -
KinFitResults
Plots the χ²/NDF and the confidence level of the kinematic fit, the latter also in a plot with a logarithmic horizontal axis.DHistogramAction_KinFitResults(dComboWrapper, "")
-
BeamEnergy
Plots the beam energy.DHistogramAction_BeamEnergy(dComboWrapper, bool, "")
-
Energy_UnusedShowers
Plot the total energy of the unused showers in the reaction.DHistogramAction_Energy_UnusedShowers(dComboWrapper, "")
Every Cut Action will skip the rest of the code if a combination does
not pass the selection and will set the flag IsComboCut
accordingly.
They are added to the dAnalysisActions
vector exactly like Histogram
Actions. Currently available cut actions are defined in gluex_root_analysis/libraries/DSelector/DCutActions.cc,
as listed here:
-
ChiSqOrCL
Special cut action to compare multiple reactions, see dedicated section KinFit Comparison between different reactions. -
PIDDeltaT
Cut on the timing of one type of particle in a certain detector.DCutAction_PIDDeltaT(dComboWrapper, bool, 0.5, PiPlus, SYS_TOF, "")
-
NoPIDHit
Cut specific particles that have no particle ID hit in a certain detector.DCutAction_NoPIDHit(dComboWrapper, PiPlus, SYS_BCAL, "")
-
dEdx
Predefined cut on the specific energy loss$dE/dx$ . Automatically distinguishes between heavy particles (mass >= proton mass) and light particles and uses predefined functions (might need tuning!). Only the cut on the CDC$dE/dx$ has been implemented.DCutAction_dEdx(dComboWrapper, bool, PiPlus, SYS_CDC, "")
-
KinFitFOM
Cut all combos with smaller kinematic fit confidence level.DCutAction_KinFitFOM(dComboWrapper, double locMinimumConfidenceLevel, "")
-
KinFitChiSq
Cut all combos with larger kinematic fit χ²/NDF.DCutAction_KinFitChiSq(dComboWrapper, double locMaximumChiSq, "")
-
MissingMass
This cuts the total missing mass of the reaction or the recoil mass against a particle or set of particles.DCutAction_MissingMass(dComboWrapper, bool, 0.1, 0.2, "") //OR DCutAction_MissingMass(dComboWrapper, bool, step, Pi0, 0.1, 0.2, "") //OR DCutAction_MissingMass(dComboWrapper, bool, step, {Proton, Pi0}, 0.5, 1.5, "")
-
MissingMassSquared
Cuts on the squared value of the missing mass. The syntax is the same as for the missing mass. -
MissingEnergy
Cuts on the missing energy of the reaction. Syntax is the same as for missing mass. -
InvariantMass
Cut on the invariant mass of a particle or a particle combination within the same reaction step.DCutAction_InvariantMass(dComboWrapper, bool, Pi0, 0.10, 0.15, "") //OR DCutAction_InvariantMass(dComboWrapper, bool, step, {Proton, PiMinus}, 1.0, 1.5, "")
-
InvariantMassVeto
Same as above. -
BeamEnergy
Cut on beam energy interval.DCutAction_BeamEnergy(dComboWrapper, bool, locMinBeamEnergy, locMaxBeamEnergy, "")
-
TrackShowerEOverP
Cut on$E/p$ above value for$e^+/e^-$ , below for all other particles.DCutAction_TrackShowerEOverP(dComboWrapper, bool, SYS_FCAL, Electron, 0.7, "")
-
ShowerQuality
Requires a minimal shower quality for all photons in the detector, here cut below 0.5 for FCal.DCutAction_ShowerQuality(dComboWrapper, SYS_FCAL, 0.5)
-
Kinematics
Cut on momentum, θ and φ of a particle. If the lower limit is larger than the higher limit for one variable, the cut on this variable is ignored.DCutAction_Kinematics(dComboWrapper, step, Proton, bool, MinP, MaxP, MinTheta = 1.0, MaxTheta = 0.0, MinPhi = 1.0, MaxPhi = 0.0)
-
TrackBCALPreshowerFraction
Cut on preshower fraction, weighted by$\sin\theta$ of the track. At the moment only implemented for$e^+/e^-$ .DCutAction_TrackBCALPreshowerFraction(dComboWrapper, bool, Electron, 0.1, "")
-
Energy_UnusedShowers
Cut on the summed energy of unused showers.DCutAction_Energy_UnusedShowers(dComboWrapper, locMaxEnergy_UnusedShowersCut, "")
-
NumUnusedTracks
Cut on the maximum number of unused tracks.DCutAction_NumUnusedTracks(dComboWrapper, locMaxUnusedTracks, "")
-
NumUnusedShowers
Cut on the number of unused showers.DCutAction_NumUnusedShowers(dComboWrapper, locMaxUnusedShowers, "")
-
VanHoveAngle
Select a sector in the vanHove plane.DCutAction_VanHoveAngle(dComboWrapper, bool, Proton, Eta, Pi0, MinAngle, MaxAngle, "") //OR DCutAction_VanHoveAngle(dComboWrapper, bool, {Proton, Pi0}, Eta, Pi0, MinAngle, MaxAngle)
-
VanHoveAngleFour
Select a part of the vanHove space for 4 particles, defined by 2 angles. Syntax similar as above.
The functions in this part is defined in the libraries/DSelector/DAnalysisUtilities.cc.
-
Get_IsPolarizedBeam
If the RCDB environment is set up, this function returnstrue
for polarized beam runs. It also sets a flag according toPARA/PERP
.Get_IsPolarizedBeam(locRunNumber, locIsPARAFlag)
-
Get_PolarizationAngle
If the RCDB environment is set up, this function returns the polarization angle.Get_PolarizationAngle(locRunNumber, locPolarizationAngle)
-
Get_CoherentPeak
If the RCDB environment is set up, this function returns the coherent peak energy for polarized runs, 0 for amorphous runs.Get_CoherentPeak(locRunNumber, locCoherentPeak, locIsPolarizedFlag)
-
Get_BeamBunchPeriod
If the RCDB environment is set up, this returns the beam bunch period.Get_BeamBunchPeriod(locRunNumber)
-
Get_AccidentalScalingFactor
If the RCDB environment is set up, this returns the accidental scaling factor.Get_AccidentalScalingFactor(locRunNumber, locBeamEnergy, locIsMC)
-
Get_AccidentalScalingFactorError
Caution: Uncertainty on the accidental scaling factor is not defined at this point.Get_AccidentalScalingFactorError(locRunNumber, locBeamEnergy, locIsMC)
-
Get_DeltaT_RF
Calculate the time difference with respect to the RF signal.Get_DeltaT_RF(locRunNumber, locBeamX4_Measured, dComboWrapper)
-
Get_RelativeBeamBucket
Calculate which relative beam bucket this combination is in.Get_RelativeBeamBucket(locRunNumber, locBeamX4_Measured, dComboWrapper)
-
Get_BeamEndpoint
Get the energy of the beam end point from CCDB.Get_BeamEndpoint(locRunNumber)
-
Get_ColumnTAGM
Get the tagger microscope column for this beam energy from CCDB.Get_ColumnTAGM(locRunNumber, locBeamEnergy)
-
Get_CounterTAGH
Get the tagger hodoscope counter for this beam energy from CCDB.Get_CounterTAGH(locRunNumber, locBeamEnergy)
-
Calc_ProdPlanePhi_Pseudoscalar
Calculates the φ angle for beam asymmetries of pseudoscalar meson photoproduction.Calc_ProdPlanePhi_Pseudoscalar(locBeamEnergy, Proton, locMesonP4)
-
Calc_DecayPlanePsi_Vector_2BodyDecay
Calculates decay asymmetry angle ψ for vector meson photoproduction, in the case where the vector meson decays into 2 pseudoscalar mesons (e.g. φ → 2K). Returns the scattering angle θ as well.Calc_DecayPlanePsi_Vector_2BodyDecay(locBeamEnergy, Proton, locBaryonP4, locMesonP4, locMesonProduct1P4, locDecayPlaneTheta)
-
Calc_DecayPlanePsi_Vector_3BodyDecay
Calculates decay asymmetry angle ψ for vector meson photoproduction, in the case where the vector meson decays into 3 pseudoscalar mesons (e.g. ω → 3π). Returns the scattering angle θ as well.Calc_DecayPlanePsi_Vector_3BodyDecay(locBeamEnergy, Proton, locBaryonP4, locMesonP4, locMesonProduct1P4, locMesonProduct2P4, locDecayPlaneTheta)
-
Calc_vanHoveCoord
Calculate the 2D vanHove coordinates from 3 four-vectors.Calc_vanHoveCoord(locXP4, locYP4, locZP4)
-
Calc_vanHoveCoordFour
Calculate the 3D vanHove coordinates from 4 four-vectors.Calc_vanHoveCoordFour(locVec1P4, locVec2P4, locVec3P4, locVec4P4)
This is the main event loop which is called for each event. There are three main sections in the default code which is provided. The first is for any event-level initialization or calculations. The second comprises the loop over all available combinations in the event. The third is for additional event-level calculations and output.
The data for the event is read from the tree with the following code section:
//CALL THIS FIRST
DSelector::Process(locEntry); //Gets the data from the tree for the entry
Note that this method Process()
is inherited from the ROOT class
TSelector
. This method is overridden by the DSelector
class and
calls the method Get_Entry()
, defined in the class DTreeInterface
,
which reads the event data from file into memory. Remember that if a
TChain
is used and a new tree file is opened, the method Init()
will
be executed again first.
If you are using analysis actions, then they need to be initialized for each event:
/*** SETUP UNIQUENESS TRACKING ****/
//ANALYSIS ACTIONS: Reset uniqueness tracking for each action
//For any actions that you are executing manually, be sure to call Reset_NewEvent() on them here
Reset_Actions_NewEvent();
dAnalyzeCutActions->Reset_NewEvent(); // manual action, must call Reset_NewEvent()
If you are filling histograms yourself, the combinatorics must be carefully handled to avoid double-counting, either by "uniqueness tracking" or some other method.
Specifically, the uniquess tracking is implemented with several types of containers to avoid double-counting with different criteria defined to distinguish one combo from the other. Before looping over different combos, the uniqueness tracking is set up in the following examples:
-
Best kinfit $\chi^2$ combo per event
For this criteria, the uniquess tracking container is simply aBool_t
flag. Once the best$\chi^2$ combo is found and filled in the histogram, the event is flagged and no more histogram filling will be made://EXAMPLE 0: Event-specific info: Bool_t locUsedSoFar_Event = false; // Flag used to mark if the best chi-squared combo is filled in the histogram
-
Best kinfit $\chi^2$ combo per beam-photon candidate
For this criteria, the uniquess tracking container is a set of beam photon candidate trackID. For every beam photon candidate, only the combo with the best$\chi^2$ combo will be filled in the histogram. All other combos with the same beam photon candidate and different final-state tracks/showers will be rejected by this container hence not make it into the histogram://EXAMPLE 1: Particle-specific info: set<Int_t> locUsedSoFar_BeamEnergy; //Int_t: Unique ID for beam particles. set: easy to use, fast to search. This container is used for the "hybrid" method dealing with combinatorics.
-
Best kinfit $\chi^2$ combo per combo with different beam photon and final-state tracks/showers
For this criteria, the uniquess tracking container is a set of maps including all <Particle Type, Set(trackID/showerID)> pairs asked by the criteria to distinguish one combo from the other:
//EXAMPLE 2: Combo-specific info:
//In general: Could have multiple particles with the same PID: Use a set of Int_t's
//In general: Multiple PIDs, so multiple sets: Contain within a map
//Multiple combos: Contain maps within a set (easier, faster to search)
set<map<Particle_t, set<Int_t> > > locUsedSoFar_MissingMass;
Later when looping over each combo, the histogram shall be filled inside of if
conditions using criteria defined with the uniqueness containers defined above. For specific code example please check later subsection Histogram Filling with Uniqueness Tracking.
A pre-loop is implemented to generate an sequence of <ComboIndex, chiSq>
pairs with ascending order of the Kinfit
// Vector to store combo information
std::vector<std::pair<UInt_t, Double_t>> loc_combos;
// Pre-loop to gather kinfit ComboIndex-chiSq pairing and sort by chiSq value ascendingly
for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i)
{
dComboWrapper->Set_ComboIndex(loc_i);
Double_t locChiSq = dComboWrapper->Get_ChiSq_KinFit("");
loc_combos.push_back(std::make_pair(loc_i, locChiSq));
Double_t logKinFitCL = TMath::Log10(dComboWrapper->Get_ConfidenceLevel_KinFit(""));
cout << logKinFitCL << endl;
}
// Sort the combos by ChiSq
std::sort(loc_combos.begin(), loc_combos.end(), [](const std::pair<UInt_t, Double_t>& a, const std::pair<UInt_t, Double_t>& b) {
return a.second < b.second;
});
At this point the loop over all combos for this event is started. The definition of a combo is combining a beam photon and the target particle with the reconstructed final state particles in the GlueX detector that can from the defined reaction. The four-momenta of all these particles are then used by the kinematic fitter in the DReaction plugin to evaluate the likely-hood of this combination of particles to be the reaction that triggered this event. Note that the total number of reconstructed final state particles in the GlueX detector, neutrals and charged tracks for a given event can be more than is required by the reaction. Therefore there exist methods to get the number of "Unused" neutral and charged particles for the combo (Get_NumUnusedShowers() and Get_NumUnusedTracks()) the reflect the number of particles that are not used by the kinematic fitter.
Note that even if there a unique combination of final state particles in
an event, this event still will result in multiple combos if there are
multiple beam photons available for this event, which is often the case
for the intense GlueX photon beam. Typically, combos are kept if the
reconstructed final state particles are in time with the tagged beam
photons within 3 or 4 of the beam bunches on either side of the prompt
beam bunch. For example, if you are looking at the reaction
dComboWrapper->Set_ComboIndex(UInt_t)
function, as shown below to initialize all variables within the loop
with the appropriate values for the given combo.
Note that the DSelector can be used not only to run over the output of
the ReactionFilter plugin, but also to write out ROOT trees which
contain a set of reduced events and/or branches containing additional
information. When these trees are created, if additional selections are
being applied on the combos being analyzed, an event will be written out
if least one combo in that event passes the selections. To keep track of
which combos have been rejected and which are kept for further analysis,
we use the "ComboCut" flag. This can be set using
dComboWrapper->Set_IsComboCut(bool)
and its value retrieved using
dComboWrapper->Get_IsComboCut()
.
An example beginning of the combo loop is shown below:
for(const auto& loc_combo : loc_combos)
{
UInt_t loc_i = loc_combo.first;
//Set branch array indices for combo and all combo particles
dComboWrapper->Set_ComboIndex(loc_i);
// Is used to indicate when combos have been cut
if(dComboWrapper->Get_IsComboCut()) // Is false when tree originally created
continue; // Combo has been cut previously
// *******************************************
// now from here on out you can do your stuff!
// *******************************************
A detailed discussion of the data available through the different particle wrappers is given in Section 3. By default, code is added to configure variables to store the unique ID numbers of each particle in the combination, along with their measured and kinematically-fitted momentum 4-vectors. Equipped with the data access methods described below and the wrappers that provide the access to these methods, we can get any data from within the tree. Continuing with the earlier example of three pions and a recoil proton in the final state, we can calculate the invariant mass of the 3-pion system using either the measured momentum 4-vectors or those determined by the kinematic fitter. In the example below, we also show how to access other information, in particular the shower quality factor for neutral particles and the z-position of the associated shower in the calorimeter:
// These quantities are from the kinematic fitter
TLorentzVector locBeamP4 = dComboBeamWrapper->Get_P4();
TLorentzVector locPiPlusP4 = dPiPlusWrapper->Get_P4();
TLorentzVector locPiMinusP4 = dPiMinusWrapper->Get_P4();
TLorentzVector locProtonP4 = dProtonWrapper->Get_P4();
TLorentzVector locPhoton1P4 = dPhoton1Wrapper->Get_P4();
TLorentzVector locPhoton2P4 = dPhoton2Wrapper->Get_P4();
// 3-Pion Final state based on kinematic fit quantities:
TLorentzVector ThreePiFinalState = locPiPlusP4 + locPiMinusP4 + locPhoton1P4 + locPhoton2P4;
// These quantities are the measured quantities
TLorentzVector locBeamP4_Measured = dComboBeamWrapper->Get_P4_Measured();
TLorentzVector locPiPlusP4_Measured = dPiPlusWrapper->Get_P4_Measured();
TLorentzVector locPiMinusP4_Measured = dPiMinusWrapper->Get_P4_Measured();
TLorentzVector locProtonP4_Measured = dProtonWrapper->Get_P4_Measured();
TLorentzVector locPhoton1P4_Measured = dPhoton1Wrapper->Get_P4_Measured();
TLorentzVector locPhoton2P4_Measured = dPhoton2Wrapper->Get_P4_Measured();
// 3-Pion Final state based on measured quantities:
TLorentzVector ThreePiFinalState_Measured = locPiPlusP4_M + locPiMinusP4_M + locPhoton1P4_M + locPhoton2P4_M;
// Get the shower quality factor of the first photon and find the z-position of the shower
float qf = dPhoton1Wrapper->Get_Shower_Quality();
TLorentzVector Gamma1_Shower = dPhoton1Wrapper->Get_X4_Shower();
float Shower1Zposition = Gamma1_Shower.Z();
The reconstructed final state particle combination is associated with the timing of a particular accelerator beam bunch (which we will call the "RF time"). The measured initial state beam photon time will either be consistent, or "in-time", with this selected beam bunch, or "out-of-time" with it. While signal events are generally always associated with in-time beam photons, there is also an in-time contribution due to combinations where the beam photon is not the real photon which created the reconstructed final state particles (for example, due to multiple in-time photons or due to the correct photon not being reconstructed). There are several methods which can be used to subtract this background contribution, the most common of which is to estimate it using out-of-time beam photon combos and subtract the background by filling histograms with appropriate weights, or propagating these weights to later analysis stages.
In this method, in-time beam photons are assigned a weight of 1, while
out-of-time beam photons are assigned a weight of dComboBeamWrapper()->Get_RFTime()
method. The following code shows an
example of how to determine the RF time of the beam photon used in a
given combo and how to define the weight for the combo based on the RF
time of the beam photon. In this particular example it is expected that
4 beam bunches on either side (excluding the two nearest out-of-time bunches) of the in-time beam bunch (at DAnalysisUtilities::Get_AccidentalScalingFactor()
.
TLorentzVector locVertex = dComboBeamWrapper->Get_X4_Measured();
float locRFTime = dComboWrapper->Get_RFTime();
TLorentzVector locBeamP4 = dComboBeamWrapper->Get_P4();
.....
float DT_RF = locVertex.T() - (locRFTime + (locVertex.Z() - dTargetCenter.Z())/29.9792458);
// Now an event weight can be assigned:
// if DT_RF = +/- 2ns within zero the beam photon is in time
// within +-4, +-8, +-12, ... the beam photon is out of time
double weight = 1.;
if (TMath::Abs(DT_RF)>2.){ // TOFIX: should actually pull the real beam bunch spacing for this
weight = -0.125; // -1/8
weight *= dAnalysisUtilities->Get_AccidentalScalingFactor(); // correct for correlation effects
}
// the weight then can be used to fill a histogram with the weight
// thereby automatically subtract accidental background
// for example, dEBeam->Fill(locBeamP4.E(), weight);
\end{lstlisting}
A more sophisticated block of code is now implemented in DSelector template made by the executable MakeDSelector
:
/********************* GET COMBO RF TIMING INFO *********************/
TLorentzVector locBeamX4_Measured = dComboBeamWrapper->Get_X4_Measured();
Double_t locBunchPeriod = dAnalysisUtilities.Get_BeamBunchPeriod(Get_RunNumber());
Double_t locDeltaT_RF = dAnalysisUtilities.Get_DeltaT_RF(Get_RunNumber(), locBeamX4_Measured, dComboWrapper);
Int_t locRelBeamBucket = dAnalysisUtilities.Get_RelativeBeamBucket(Get_RunNumber(), locBeamX4_Measured, dComboWrapper); // 0 for in-time events, non-zero integer for out-of-time photons
Int_t locNumOutOfTimeBunchesInTree = 4; //YOU need to specify this number
//Number of out-of-time beam bunches in tree (on a single side, so that total number out-of-time bunches accepted is 2 times this number for left + right bunches)
Bool_t locSkipNearestOutOfTimeBunch = true; // True: skip events from nearest out-of-time bunch on either side (recommended).
Int_t locNumOutOfTimeBunchesToUse = locSkipNearestOutOfTimeBunch ? locNumOutOfTimeBunchesInTree-1:locNumOutOfTimeBunchesInTree;
Double_t locAccidentalScalingFactor = dAnalysisUtilities.Get_AccidentalScalingFactor(Get_RunNumber(), locBeamP4.E(), dIsMC); // Ideal value would be 1, but deviations require added factor, which is different for data and MC.
Double_t locAccidentalScalingFactorError = dAnalysisUtilities.Get_AccidentalScalingFactorError(Get_RunNumber(), locBeamP4.E()); // Ideal value would be 1, but deviations observed, need added factor.
Double_t locHistAccidWeightFactor = locRelBeamBucket==0 ? 1 : -locAccidentalScalingFactor/(2*locNumOutOfTimeBunchesToUse) ; // Weight by 1 for in-time events, ScalingFactor*(1/NBunches) for out-of-time
if(locSkipNearestOutOfTimeBunch && abs(locRelBeamBucket)==1) { // Skip nearest out-of-time bunch: tails of in-time distribution also leak in
dComboWrapper->Set_IsComboCut(true);
continue;
}
In later stage of the analysis, the weight, locHistAccidWeightFactor
can be used for accidental subtracted yields when filling histograms.
There are several methods to output information in the DSelector
with
combinatoric issue dealt with by uniqueness tracking containers:
-
Hybrid Method: It uses the uniqueness tracking for The best kinfit
$\chi^2$ combo per beam photon candidate.
//Histogram beam energy (if haven't already)
if(locUsedSoFar_BeamEnergy.find(locBeamID) == locUsedSoFar_BeamEnergy.end())
{
dHist_BeamEnergy->Fill(locBeamP4.E(),locHistAccidWeightFactor); // Alternate version with accidental subtraction
locUsedSoFar_BeamEnergy.insert(locBeamID);
}
For cross-section measurement, this method should be implemented to produce the raw yield before doing flux normalization. The tagged flux used for the flux normalization can be produced using a tool available at JeffersonLab/hd_utilities/psflux.
-
Best
$\chi^2$ Method: For every event only keep the in-time combo with the best Kinfit$\chi^2$ . This is a fast straightforward method that allows the histogram/flat-tree entries numbers to be exactly the same as the number of events. Here is an example of filling a histogram with beam energy$E_{\gamma}$ using the best$\chi^2$ in-time combo for every event:
// Need to uncomment the section computing combo timing info before running this block of code
if(locUsedSoFar_Event == false)
{
//Fill the histogram only when the beam bunch is in-time.
if(!locRelBeamBucket)
{
dHist_BeamEnergy_BestChiSq->Fill(locBeamP4.E());
locUsedSoFar_Event = true;
}
}
Note that this method is not suitable for the cross section measurement, as in its yield the tagger accidental background is not properly subtracted for in-time photons. This will cause unstable result with large systematic uncertainty in flux-normalization where the yield is divided by the tagged flux, while the tagged flux has already being subtracted with the random accidental background.
The ReactionFilter ROOT trees are generated with no constraint on the
quality of the kinematic fit by default. This quantity can be very
powerful for separating signal from background. This information can be
accessed in one of two ways: by looking at the
// Kinematic fit confidence level and chi^2/d.o.f.
double chi2dof = dComboWrapper->Get_ConfidenceLevel_KinFit("") / dComboWrapper->Get_NDF_KinFit("");
double chi2dof = dComboWrapper->Get_ChiSq_KinFit("") / dComboWrapper->Get_NDF_KinFit("");
// reject events with a poor chi^2 for this reaction
if(chi2dof > 5.) {
dComboWrapper->Set_IsComboCut(true); // reject this combo for future consideration
continue;
}
As previously discussed, particle combinations can have multiple steps
in order to model the effect of additional decays of narrow resonances.
These decaying particles can be short-lived (
The cases in which the decaying particle is and is not mass constrained in the fit are handled differently. If no mass constraint is applied in the fit, first obtain the production X4 (vertex spatial position and time) and the X4 for the decaying particle. The pathlength of the decaying particle is the magnitude of the difference between the two vertices. The uncertainty of the pathlength is saved directly to the tree, but must be extracted by hand. The flight significance is then the pathlength divided by the uncertainty in the pathlength.
TLorentzVector locProdSpacetimeVertex = dComboBeamWrapper->Get_X4();
TLorentzVector locDecayingXiX4 = dTreeInterface->Get_TObject<TLorentzVector>("DecayingXiMinus__X4",loc_i); // change if the mass is constrained in the fit
TLorentzVector locDeltaSpacetimeXi = locProdSpacetimeVertex - locDecayingXiX4;
double locPathLengthXi = locDeltaSpacetimeXi.Vect().Mag();
float locPathLengthSigmaXi = Get_Fundamental<Float_t>("DecayingXiMinus__PathLengthSigma", loc_i);
double locPathLengthSignificanceXi = locPathLengthXi/locPathLengthSigmaXi;
If a mass constraint on the decaying particle was applied, then a wrapper for the decaying particle is configured and should be used to obtain the vertex position instead:
TLorentzVector locDecayingXiX4 = dDecayingXiMinusWrapper->Get_X4();
For a secondary detached vertex, for example the
TLorentzVector locDecayingLambX4 = dTreeInterface->Get_TObject<TLorentzVector>("DecayingLambda__X4",loc_i);
TLorentzVector locDeltaSpacetimeLamb = locDecayingXiX4 - locDecayingLambX4;
double locPathLengthLamb = locDeltaSpacetimeLamb.Vect().Mag();
float locPathLengthSigmaLamb = Get_Fundamental<Float_t>("DecayingLambda__PathLengthSigma", loc_i);
double locPathLengthSignificanceLamb = locPathLengthLamb/locPathLengthSigmaLamb;
When analyzing simulated data (MC), information about the particles
thrown by the physics event generator is also stored for each event. The
data is stored in objects of type DMCThrown
and is also accessible
through a wrapper interface. This information can be accessed in several
different ways. Each DChargedTrackHypothesis
and
DNeutralParticleHypothesis
has a member function named
Get_ThrownIndex()
, which can be passed as an argument to
dThrownWrapper->Set_ArrayIndex()
. The list of thrown particles can
also be studied directly, such as in the following example:
//Thrown beam: just use directly
if(dThrownBeam != NULL)
double locEnergy = dThrownBeam->Get_P4().E();
//Loop over throwns
for(UInt_t loc_i = 0; loc_i < Get_NumThrown(); ++loc_i)
{
//Set branch array indices corresponding to this particle
dThrownWrapper->Set_ArrayIndex(loc_i);
//Do stuff with the wrapper here ...
//Example - fill a histogram of the momentum distribution of protons
if(dThrownWrapper->Get_PID() == Proton)
dProtonPmag->Fill(dThrownWrapper->Get_P4().Vect().Mag());
}
The DSelector can be used to create TTrees for further analysis, either by another DSelector or by another program. Two types of TTree output are support: the standard ReactionFilter output in PART format, which stores data mostly in various arrays of objects, and a "flat" tree with a simpler structure.
The code to handle ROOT tree output is in several locations, and
generally works through the (thread-safe) DTreeInterface class. In the
Init()
function, besides defining the names of output files, you can
optionally add branches to store additional event/combo information. The
following shows an example of these features:
// In Init() function...
//USERS: SET OUTPUT FILE NAME //can be overriden by user in PROOF
dOutputFileName = "ggg.root"; //"" for none
dOutputTreeFileName = ""; //"" for none
dFlatTreeFileName = ""; //output flat tree (one combo per tree entry), "" for none
dFlatTreeName = ""; //if blank, default name will be chosen
/************************** EXAMPLE USER INITIALIZATION: CUSTOM OUTPUT BRANCHES - MAIN TREE *************************/
//EXAMPLE MAIN TREE CUSTOM BRANCHES (OUTPUT ROOT FILE NAME MUST FIRST BE GIVEN!!!! (ABOVE: TOP)):
//The type for the branch must be included in the brackets
//1st function argument is the name of the branch
//2nd function argument is the name of the branch that contains the size of the array (for fundamentals only)
dTreeInterface->Create_Branch_Fundamental<Int_t>("my_int"); //fundamental = char, int, float, double, etc.
dTreeInterface->Create_Branch_FundamentalArray<Int_t>("my_int_array", "my_int");
dTreeInterface->Create_Branch_FundamentalArray<Float_t>("my_combo_array", "NumCombos");
dTreeInterface->Create_Branch_NoSplitTObject<TLorentzVector>("my_p4");
dTreeInterface->Create_Branch_ClonesArray<TLorentzVector>("my_p4_array");
/************************** EXAMPLE USER INITIALIZATION: CUSTOM OUTPUT BRANCHES - FLAT TREE *************************/
//EXAMPLE FLAT TREE CUSTOM BRANCHES (OUTPUT ROOT FILE NAME MUST FIRST BE GIVEN!!!! (ABOVE: TOP)):
//The type for the branch must be included in the brackets
//1st function argument is the name of the branch
//2nd function argument is the name of the branch that contains the size of the array (for fundamentals only)
dFlatTreeInterface->Create_Branch_Fundamental<Int_t>("flat_my_int"); //fundamental = char, int, float, double, etc.
dFlatTreeInterface->Create_Branch_FundamentalArray<Int_t>("flat_my_int_array", "flat_my_int");
dFlatTreeInterface->Create_Branch_NoSplitTObject<TLorentzVector>("flat_my_p4");
dFlatTreeInterface->Create_Branch_ClonesArray<TLorentzVector>("flat_my_p4_array");
In the Process()
function, these trees can be filled either once per
event (i.e. before the combo loop)
/**************************************** EXAMPLE: FILL CUSTOM OUTPUT BRANCHES **************************************/
Int_t locMyInt = 7;
dTreeInterface->Fill_Fundamental<Int_t>("my_int", locMyInt);
TLorentzVector locMyP4(4.0, 3.0, 2.0, 1.0);
dTreeInterface->Fill_TObject<TLorentzVector>("my_p4", locMyP4);
for(int loc_i = 0; loc_i < locMyInt; ++loc_i)
dTreeInterface->Fill_Fundamental<Int_t>("my_int_array", 3*loc_i, loc_i); //2nd argument = value, 3rd = array index
or once per combo (i.e. inside the combo loop)
/**************************************** EXAMPLE: FILL CUSTOM OUTPUT BRANCHES **************************************/
TLorentzVector locMyComboP4(8.0, 7.0, 6.0, 5.0);
//for arrays below: 2nd argument is value, 3rd is array index
//NOTE: By filling here, AFTER the cuts above, some indices won't be updated (and will be whatever they were from the last event)
//So, when you draw the branch, be sure to cut on "IsComboCut" to avoid these.
dTreeInterface->Fill_Fundamental<Float_t>("my_combo_array", -2*loc_i, loc_i);
dTreeInterface->Fill_TObject<TLorentzVector>("my_p4_array", locMyComboP4, loc_i);
For the normal tree output, the tree should be written out at the end of
Process()
:
/************************************ EXAMPLE: FILL CLONE OF TTREE HERE WITH CUTS APPLIED ************************************/
Bool_t locIsEventCut = true;
for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) {
//Set branch array indices for combo and all combo particles
dComboWrapper->Set_ComboIndex(loc_i);
// Is used to indicate when combos have been cut
if(dComboWrapper->Get_IsComboCut())
continue;
locIsEventCut = false; // At least one combo succeeded
break;
}
if(!locIsEventCut && dOutputTreeFileName != "")
Fill_OutputTree();
For a flat tree, this format has one entry per combo, and so it should be filled and written out right at the end of the combo loop, for example:
/****************************************** FILL FLAT TREE (IF DESIRED) ******************************************/
//FILL ANY CUSTOM BRANCHES FIRST!!
Int_t locMyInt_Flat = 7;
dFlatTreeInterface->Fill_Fundamental<Int_t>("flat_my_int", locMyInt_Flat);
TLorentzVector locMyP4_Flat(4.0, 3.0, 2.0, 1.0);
dFlatTreeInterface->Fill_TObject<TLorentzVector>("flat_my_p4", locMyP4_Flat);
for(int loc_j = 0; loc_j < locMyInt_Flat; ++loc_j)
{
dFlatTreeInterface->Fill_Fundamental<Int_t>("flat_my_int_array", 3*loc_j, loc_j); //2nd argument = value, 3rd = array index
TLorentzVector locMyComboP4_Flat(8.0, 7.0, 6.0, 5.0);
dFlatTreeInterface->Fill_TObject<TLorentzVector>("flat_my_p4_array", locMyComboP4_Flat, loc_j);
}
//FILL FLAT TREE
Fill_FlatTree(); //for the active combo