Difference between revisions of "Analysis TTreeFormat"

#include "ANALYSIS/DEventWriterROOT.h"
//In plugin brun():
const DEventWriterROOT* locEventWriterROOT = NULL;
locEventLoop->GetSingle(locEventWriterROOT);
locEventWriterROOT->Create_DataTrees(locEventLoop); //creates TTrees for all output-enabled DReactions
locEventWriterROOT->Create_ThrownTree("tree_b1pi_thrownmc.root"); //optional: create a ttree containing only the thrown data //string is output file name

//In plugin evnt()
const DEventWriterROOT* locEventWriterROOT = NULL;
locEventLoop->GetSingle(locEventWriterROOT);
locEventWriterROOT->Fill_DataTrees(locEventLoop, "b1pi_hists"); //string is the DReaction factory tag that the DReactions were created in

//In plugin evnt()
#include "ANALYSIS/DEventWriterROOT.h"
vector<const DEventWriterROOT*> locEventWriterROOTVector;
locEventLoop->Get(locEventWriterROOTVector); //creates the TTrees for all DReactions upon first call
locEventWriterROOTVector[0]->Fill_Tree(locEventLoop, locReaction, locParticleCombos);

//In plugin evnt()
const DEventWriterROOT* locEventWriterROOT = NULL;
locEventLoop->GetSingle(locEventWriterROOT);
locEventWriterROOT->Fill_ThrownTree(locEventLoop);

MyTree->Draw("PiMinus1__P4_Measured->Theta()"); //draws all particle combinations

GetEntry(entry);
 
//this doesn't take into account double-counting!!
TLorentzVector locPiZeroP4 = *Gamma1__P4_KinFit + *Gamma2__P4_KinFit;
dPiZeroMassHist->Fill(locPiZeroP4.M());
 
TLorentzVector locOmegaP4 = locPiZeroP4 + *PiPlus2__P4_KinFit + *PiMinus2__P4_KinFit;
dOmegaMassHist->Fill(locOmegaP4.M());
 
TLorentzVector locB1PlusP4 = locOmegaP4 + *PiPlus1__P4_KinFit;
dB1PlusMassHist->Fill(locB1PlusP4.M());
 
TLorentzVector locX2000P4 = locB1PlusP4 + *PiMinus1__P4_KinFit;
dX2000MassHist->Fill(locX2000P4.M());

//member variables added to your TSelector
vector<vector<int> > sPreviousObjectIDs_Beam; //1st dimension is cut-index
vector<vector<int> > sPreviousObjectIDs_PiMinus1; //1st dimension is cut-index
vector<vector<int> > sPreviousObjectIDs_PiMinus2; //1st dimension is cut-index
vector<vector<int> > sPreviousObjectIDs_PiPlus1; //1st dimension is cut-index
vector<vector<int> > sPreviousObjectIDs_PiPlus2; //1st dimension is cut-index
vector<vector<int> > sPreviousObjectIDs_Gamma1; //1st dimension is cut-index
vector<vector<int> > sPreviousObjectIDs_Gamma2; //1st dimension is cut-index

//sNumCuts is the number of cuts you're placing in your TSelector
if((RunNumber != sPreviousRunNumber) || (EventNumber != sPreviousEventNumber))
{
  sPreviousObjectIDs_Beam.clear();
  sPreviousObjectIDs_PiMinus1.clear();
  sPreviousObjectIDs_PiMinus2.clear();
  sPreviousObjectIDs_PiPlus1.clear();
  sPreviousObjectIDs_PiPlus2.clear();
  sPreviousObjectIDs_Gamma1.clear();
  sPreviousObjectIDs_Gamma2.clear();
  sPreviousObjectIDs_Beam.resize(sNumCuts + 1);
  sPreviousObjectIDs_PiMinus1.resize(sNumCuts + 1);
  sPreviousObjectIDs_PiMinus2.resize(sNumCuts + 1);
  sPreviousObjectIDs_PiPlus1.resize(sNumCuts + 1);
  sPreviousObjectIDs_PiPlus2.resize(sNumCuts + 1);
  sPreviousObjectIDs_Gamma1.resize(sNumCuts + 1);
  sPreviousObjectIDs_Gamma2.resize(sNumCuts + 1);
}

size_t locCutIndex = 0;
if(Gamma1__Energy_BCAL < 0.5)
  return End_Entry(locCutIndex);
++locCutIndex;
 
Bool_t MySelector::End_Entry(size_t locCutIndex)
{
  sPreviousRunNumber = RunNumber;
  sPreviousEventNumber = EventNumber;
  for(size_t loc_i = 0; loc_i <= locCutIndex; ++loc_i)
  {
    sPreviousObjectIDs_Beam[loc_i].push_back(Beam__ObjectID);
    sPreviousObjectIDs_PiMinus1[loc_i].push_back(PiMinus1__ObjectID);
    sPreviousObjectIDs_PiMinus2[loc_i].push_back(PiMinus2__ObjectID);
    sPreviousObjectIDs_PiPlus1[loc_i].push_back(PiPlus1__ObjectID);
    sPreviousObjectIDs_PiPlus2[loc_i].push_back(PiPlus2__ObjectID);
    sPreviousObjectIDs_Gamma1[loc_i].push_back(Gamma1__ObjectID);
    sPreviousObjectIDs_Gamma2[loc_i].push_back(Gamma2__ObjectID);
  }
  return kTRUE;
}

bool locDuplicateFlag = false;
for(size_t loc_i = 0; loc_i < sPreviousObjectIDs_Gamma1[locCutIndex].size(); ++loc_i)
{
  if((sPreviousObjectIDs_Gamma1[locCutIndex][loc_i] == Gamma1__ObjectID) && (sPreviousObjectIDs_Gamma2[locCutIndex][loc_i] == Gamma2__ObjectID))
  {
    locDuplicateFlag = true;
    break;
  }
  if((sPreviousObjectIDs_Gamma1[locCutIndex][loc_i] == Gamma2__ObjectID) && (sPreviousObjectIDs_Gamma2[locCutIndex][loc_i] == Gamma1__ObjectID))
  {
    locDuplicateFlag = true; //g1 & g2 are switched but the invariant mass will be identical
    break;
  }
}
if(!locDuplicateFlag)
  MyPiZeroMassHist->Fill(locPiZeroMass);

// EVENT DATA
"RunNumber": UInt_t
"EventNumber": UInt_t
 
// # PARTICLES //these are the array sizes for many of the other branches
"NumThrown": UInt_t
"NumUnused": UInt_t
 
// RF
"RFTime_Thrown": Double_t
"RFTime_Measured": Double_t
"RFTime_KinFit": Double_t //only if kinematic fit performed
 
// KINEMATIC FIT
"ChiSq_KinFit": Double_t //only if kinematic fit performed
"NDF_KinFit": UInt_t //only if kinematic fit performed
 
//THROWN PARTICLES BY PID
"NumPIDThrown_FinalState": ULong64_t //the # of thrown final-state particles (+ pi0) of each type (multiplexed in base 10)
                                       //types (in order from 10^0 -> 10^15): g, e+, e-, nu, mu+, mu-, pi0, pi+, pi-, KLong, K+, K-, n, p, p-bar, n-bar
                                       //e.g. particles decaying from final-state particles are NOT included (e.g. photons from pi0, muons from pions, etc.)
                                     //is sum of #-of-PID * 10^ParticleMultiplexPower() (defined in libraries/include/particleType.h)
                                     //ParticleMultiplexPower() returns a different power of 10 for each final-state PID type. 
                                     //A value of 9 should be interpreted as >= 9.  
"PIDThrown_Decaying": ULong64_t //the types of the thrown decaying particles in the event (multiplexed in base 2)
                                //binary power of a PID is given by ParticleMultiplexPower() (defined in libraries/include/particleType.h)
                                //types: most Particle_t's that aren't final state (e.g. lambda, eta, phi, rho0, etc.) see ParticleMultiplexPower()

"X4_Production": TLorentzVector //only if beam particle is present in DReaction
"X4_LambdaDecay": TLorentzVector //only if decaying particles with detached vertices are present in DReaction //Lambda is just an example, could be KShort, etc.

//IDENTIFIERS
"ObjectID": UInt_t["NumThrown"] //each particle has its own #
"ParentID": UInt_t["NumThrown"] //the "ObjectID" of the particle this particle decayed from (0 if none (e.g. photoproduced))
"PID_PDG": Int_t["NumThrown"]
 
//KINEMATICS: THROWN  //At the production vertex 
"X4_Thrown": TClonesArray(TLorentzVector["NumThrown"])
"P4_Thrown": TClonesArray(TLorentzVector["NumThrown"])

//IDENTIFIER
"ObjectID": Int_t //each beam particle has its own #
 
//KINEMATICS: MEASURED //At the target center
"X4_Measured": TLorentzVector
"P4_Measured": TLorentzVector
 
//KINEMATICS: KINFIT //At the interaction vertex //only present if kinfit performed
"X4_KinFit": TLorentzVector
"P4_KinFit": TLorentzVector

//IDENTIFIER
"ObjectID": ULong64_t //each physical particle has its own # (to keep track of different pid hypotheses for the same particle)
 
//KINEMATICS: MEASURED //At the production vertex
"X4_Measured": TLorentzVector //t is the measured value in TOF/BCAL/FCAL projected back to Position_Measured
"P4_Measured": TLorentzVector
 
// KINEMATICS: END //at the reconstructed position of the BCAL/FCAL/TOF hit
"X4_End": TLorentzVector
"P4_End": TLorentzVector
 
//KINEMATICS: KINFIT //At the production vertex //only present if kinfit performed
"X4_KinFit": TLorentzVector
"P4_KinFit": TLorentzVector
 
// KINEMATICS: OTHER
"PathLength": Double_t //from dPosition_KinFit/Measured (KinFit if performed) to dPosition_End
 
// PID QUALITY:
"NDF_Tracking": UInt_t //for charged only
"ChiSq_Tracking": Double_t //for charged only
"NDF_Timing": UInt_t
"ChiSq_Timing_Measured": Double_t //using measured data
"ChiSq_Timing_KinFit": Double_t //using kinematic fit data //only present if kinfit performed
"NDF_DCdEdx": UInt_t //for charged only
"ChiSq_DCdEdx": Double_t //for charged only
 
// DEPOSITED ENERGY:
"dEdx_CDC": Double_t //for charged only
"dEdx_FDC": Double_t //for charged only
"dEdx_TOF": Double_t //for charged only
"dEdx_ST": Double_t //for charged only
"Energy_BCAL": Double_t
"Energy_FCAL": Double_t

//IDENTIFIERS
"ObjectID": ULong64_t["NumUnused"] //each physical particle has its own # (to keep track of different pid hypotheses for the same particle)
"PID": UInt_t["NumUnused"] //Particle_t value
 
//KINEMATICS: MEASURED  //At the production vertex 
"P4_Measured": TClonesArray(TLorentzVector["NumUnused"])
"X4_Measured": TClonesArray(TLorentzVector["NumUnused"]) //t is the measured value in TOF/BCAL/FCAL projected back to Position_Measured
 
// KINEMATICS: END //at the reconstructed position of the BCAL/FCAL/TOF hit
"P4_End": TClonesArray(TLorentzVector["NumUnused"])
"X4_End": TClonesArray(TLorentzVector["NumUnused"])
 
// KINEMATICS: OTHER
"PathLength": Double_t["NumUnused"] //from Position_Measured to Position_End
 
// PID QUALITY:
"NDF_Tracking": UInt_t["NumUnused"] //for charged only
"ChiSq_Tracking": Double_t["NumUnused"] //for charged only
"NDF_Timing": UInt_t["NumUnused"]
"ChiSq_Timing": Double_t["NumUnused"]
"NDF_DCdEdx": UInt_t["NumUnused"] //for charged only
"ChiSq_DCdEdx": Double_t["NumUnused"] //for charged only
 
// DEPOSITED ENERGY:
"dEdx_CDC": Double_t["NumUnused"] //for charged only
"dEdx_FDC": Double_t["NumUnused"] //for charged only
"dEdx_TOF": Double_t["NumUnused"] //for charged only
"dEdx_ST": Double_t["NumUnused"] //for charged only
"Energy_BCAL": Double_t["NumUnused"]
"Energy_FCAL": Double_t["NumUnused"]

class DTreeParticle : public TObject
{
  public:
    // PID:
    Particle_t dPID;
 
    // KINEMATICS: Measured
    TVector3 dPosition_Measured; //the position where the particle is produced
    Double_t dTime_Measured; //time of the track at dPosition_Measured, projected from the TOF/BCAL/FCAL hit
    TVector3 dMomentum_Measured; //momentum of the track at dPosition_Start
 
    // KINEMATICS: KinFit 
    TVector3 dPosition_KinFit; //the position where the particle is produced
    Double_t dTime_KinFit; //time of the track at dPosition_KinFit
    TVector3 dMomentum_KinFit; //momentum of the track at dPosition_KinFit
 
    // KINEMATICS: End
    TVector3 dPosition_End; //detected particles: the reconstructed position of the BCAL/FCAL/TOF hit; decaying particles: the point where it decays
    Double_t dTime_End; //time of the track at dPosition_End
    TVector3 dMomentum_End; //momentum of the track at dPosition_End
 
    // KINEMATICS: Other
    Double_t dPathLength; //from dPosition_KinFit/Measured (if kinfit performed/not-performed) to dPosition_End
 
    // PID QUALITY:
    UInt_t dNDF_Tracking;
    Double_t dChiSq_Tracking;
    UInt_t dNDF_Timing;
    Double_t dChiSq_Timing;
    UInt_t dNDF_DCdEdx;
    Double_t dChiSq_DCdEdx;
 
    // DEPOSITED ENERGY:
    map<DetectorSystem_t, Double_t> dDepositedEnergies; //Is dE/dx for all but BCAL/FCAL (shower energies)
 
    // DTREESTEP POINTERS:
    const DTreeStep* dProductionStep; //the step object in which this DTreeParticle is produced (is a final-state particle)
    const DTreeStep* dDecayStep; //the step object in which this DTreeParticle decays (is an initial-state particle) (will be null if not a decaying particle!)
 
    // CUSTOM VARIABLES:
    map<string, Double_t> dCustomVariables; //key is unique string, Double_t is value
    map<string, const TObject*> dCustomObjects; //key is unique string, TObject* is object
 
  ClassDef(DTreeParticle, 1)
};

class DTreeStep : public TObject
{
  public:
    // INITIAL PARTICLES:
    const DTreeParticle* dInitialParticle; //if is null: decaying or beam particle not yet set!
    const DTreeParticle* dTargetParticle; //NULL for no target
 
    // FINAL PARTICLES:
    deque<DParticle_t> dFinalParticleIDs; //separate in case particle is NULL (e.g. decaying resonance)
    deque<const DTreeParticle*> dFinalParticles; //particle may be NULL if it is decaying or missing (especially if no kinematic fit was performed!!)
 
    // CUSTOM VARIABLES:
    map<string, Double_t> dCustomVariables; //key is unique string, Double_t is value
    map<string, const TObject*> dCustomObjects; //key is unique string, TObject* is object
 
  ClassDef(DTreeStep, 1)
};

class DTreeCombo : public TObject
{
  public:
    // STEPS:
    deque<const DTreeStep*> dTreeSteps;
 
    // RUN, EVENT #'s:
    UInt_t dRunNumber;
    UInt_t dEventNumber;
 
    // RF:
    Double_t dRFTime;
    Double_t dRFTime_KinFit;
 
    // KINEMATIC FIT:
    DKinFitType dKinematicFitType; //Defined in DKinFitResults.h //d_NoFit if not performed
    Double_t dChiSq_KinematicFit; //NaN if not performed
    UInt_t dNDF_KinematicFit; //0 if not performed
 
    // UNUSED PARTICLES:
    vector<const DTreeParticle*> dUnusedDetectedParticles;
    vector<const DTreeShower*> dUnusedDetectedShowers;
 
    // THROWNS:
    Double_t dRFTime_Thrown;
    deque<const DTreeParticle*> dThrownParticles;
 
    // CUSTOM VARIABLES:
    map<string, Double_t> dCustomVariables; //key is unique string, Double_t is value
    map<string, const TObject*> dCustomObjects; //key is unique string, TObject* is object
 
  ClassDef(DTreeCombo, 1)
};

//Create the trees
vector<const DEventWriterROOT*> locEventWriterROOTVector;
locEventLoop->Get(locEventWriterROOTVector); //creates the TTrees for all DReactions upon first call
 
//Create the branch
Double_t* locMyVariable = new Double_t;
string locReactionName = locReaction->Get_ReactionName();
string locTreeName = locReactionName + string("_Tree");
japp->RootWriteLock(); //always acquire a lock before accessing the global ROOT file
{
  //get the tree
  gDirectory->cd("/");
  gDirectory->cd(locReactionName.c_str());
  TTree* locTree = (TTree*)gDirectory->Get(locTreeName.c_str());
  locTree->Branch("MyBranchName", locMyVariable, "D");
}
japp->RootUnLock();
 
//Set the data in the branch and fill
japp->RootWriteLock(); //lock: any thread has access to the tbranch
{
  *locMyVariable = 3.0; //or whatever you want to save
}
japp->RootUnLock();
locEventWriterROOTVector[0]->Fill_Trees(locEventLoop);