epic-analysis/src/AnalysisEpic.cxx

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Gregory Matousek, Christopher Dilks
0003
0004 #include "AnalysisEpic.h"
0005
0006 AnalysisEpic::AnalysisEpic(TString infileName_, TString outfilePrefix_)
0007   : Analysis(infileName_, outfilePrefix_)
0008 {};
0009
0010 AnalysisEpic::~AnalysisEpic() {};
0011
0012 void AnalysisEpic::Execute()
0013 {
0014   // setup
0015   Prepare();
0016
0017   // read EventEvaluator tree
0018   auto chain = std::make_unique<TChain>("events");
0019   for(Int_t idx=0; idx<infiles.size(); ++idx) {
0020     for(std::size_t idxF=0; idxF<infiles[idx].size(); ++idxF) {
0021       // std::cout << "Adding " << infiles[idx][idxF] << std::endl;
0022       chain->Add(infiles[idx][idxF].c_str());
0023     }
0024   }
0025   chain->CanDeleteRefs();
0026   auto listOfBranches = chain->GetListOfBranches();
0027
0028   TTreeReader tr(chain.get());
0029
0030   TTreeReaderArray<Int_t> hepmcp_status(tr, "GeneratedParticles.type");
0031   TTreeReaderArray<Int_t> hepmcp_PDG(tr,    "GeneratedParticles.PDG");
0032   TTreeReaderArray<Float_t> hepmcp_E(tr,      "GeneratedParticles.energy");
0033   TTreeReaderArray<Float_t> hepmcp_psx(tr,    "GeneratedParticles.momentum.x");
0034   TTreeReaderArray<Float_t> hepmcp_psy(tr,    "GeneratedParticles.momentum.y");
0035   TTreeReaderArray<Float_t> hepmcp_psz(tr,    "GeneratedParticles.momentum.z");
0036
0037
0038
0039   // All true particles (including secondaries, etc)
0040   TTreeReaderArray<Int_t> mcpart_PDG(tr,      "MCParticles.PDG");
0041   TTreeReaderArray<Int_t> mcpart_genStat(tr,         "MCParticles.generatorStatus");
0042   TTreeReaderArray<Int_t> mcpart_simStat(tr,         "MCParticles.simulatorStatus");
0043   TTreeReaderArray<Double_t> mcpart_m(tr,         "MCParticles.mass");
0044   TTreeReaderArray<Float_t> mcpart_psx(tr,       "MCParticles.momentum.x");
0045   TTreeReaderArray<Float_t> mcpart_psy(tr,       "MCParticles.momentum.y");
0046   TTreeReaderArray<Float_t> mcpart_psz(tr,       "MCParticles.momentum.z");
0047
0048
0049   // Reco tracks
0050   TTreeReaderArray<Int_t> recparts_type(tr,  "ReconstructedChargedParticles.type"); // needs to be made an int eventually in actual EE code
0051   TTreeReaderArray<Float_t> recparts_e(tr, "ReconstructedChargedParticles.energy");
0052   TTreeReaderArray<Float_t> recparts_p_x(tr, "ReconstructedChargedParticles.momentum.x");
0053   TTreeReaderArray<Float_t> recparts_p_y(tr, "ReconstructedChargedParticles.momentum.y");
0054   TTreeReaderArray<Float_t> recparts_p_z(tr, "ReconstructedChargedParticles.momentum.z");
0055   TTreeReaderArray<Int_t> recparts_PDG(tr,  "ReconstructedChargedParticles.PDG");
0056   TTreeReaderArray<Float_t> recparts_CHI2PID(tr,  "ReconstructedChargedParticles.goodnessOfPID");
0057
0058   // RecoAssociations
0059   std::string assoc_branch_name = "ReconstructedChargedParticleAssociations";
0060   if(listOfBranches->FindObject(assoc_branch_name.c_str()) == nullptr)
0061     assoc_branch_name = "ReconstructedChargedParticlesAssociations"; // productions before 23.5
0062   TTreeReaderArray<UInt_t> assoc_simID(tr, (assoc_branch_name+".simID").c_str());
0063   TTreeReaderArray<UInt_t> assoc_recID(tr, (assoc_branch_name+".recID").c_str());
0064   TTreeReaderArray<Float_t> assoc_weight(tr, (assoc_branch_name+".weight").c_str());
0065
0066   // calculate Q2 weights
0067   CalculateEventQ2Weights();
0068
0069   // counters
0070   Long64_t numNoBeam, numEle, numNoEle, numNoHadrons, numProxMatched;
0071   numNoBeam = numEle = numNoEle = numNoHadrons = numProxMatched = 0;
0072
0073
0074
0075   cout << "begin event loop..." << endl;
0076   tr.SetEntriesRange(1,maxEvents);
0077   do{
0078     if(tr.GetCurrentEntry()%10000==0) cout << tr.GetCurrentEntry() << " events..." << endl;
0079     total_events[0]++; // increment number of events
0080     // resets
0081     kin->ResetHFS();
0082     kinTrue->ResetHFS();
0083
0084
0085     double maxP = 0;
0086     int genEleID = -1;
0087     bool foundBeamElectron = false;
0088     bool foundBeamIon = false;
0089
0090     // Index maps for particle sets
0091     std::map <double,int> genidmap; // <pz, index_gen>
0092     std::map <int,int> mcidmap; // <index_mc, index_gen>
0093     std::map <int,int> recidmap; // <index, index_mc>
0094     std::map <int,int> trackstatmap; // <index, genstatus>
0095
0096     // Particles vectors
0097     // The index of the vectors correspond to their for loop idx
0098     std::vector<Particles> genpart;    // mcID --> igen
0099     std::vector<Particles> mcpart;     // mcID --> imc
0100     std::vector<Particles> recpart;  // mcID --> (imc of matching mcpart) or (-1 if no match is found)
0101
0102     /*
0103       GenParticles loop
0104     */
0105
0106     for(int igen=0; igen<hepmcp_PDG.GetSize(); igen++) {
0107
0108       int pid_ = hepmcp_PDG[igen];
0109
0110       double px_ = hepmcp_psx[igen];
0111       double py_ = hepmcp_psy[igen];
0112       double pz_ = hepmcp_psz[igen];
0113       double e_  = hepmcp_E[igen];
0114
0115       double p_ = sqrt(pow(hepmcp_psx[igen],2) + pow(hepmcp_psy[igen],2) + pow(hepmcp_psz[igen],2));
0116       double mass_ = (fabs(pid_)==211)?constants::pimass:(fabs(pid_)==321)?constants::kmass:(fabs(pid_)==11)?constants::emass:(fabs(pid_)==13)?constants::mumass:(fabs(pid_)==2212)?constants::pmass:0.;
0117
0118       // Add to genpart
0119       Particles part;
0120
0121       part.pid=pid_;
0122       // part.charge = // TODO; not used yet
0123       part.mcID=igen;
0124       part.vecPart.SetPxPyPzE(px_,py_,pz_,e_);
0125       genpart.push_back(part);
0126
0127       genidmap.insert({pz_,igen});
0128
0129     }
0130
0131     /*
0132       MCParticles loop
0133     */
0134
0135     for(int imc=0; imc < mcpart_PDG.GetSize(); imc++){
0136
0137       int pid_ = mcpart_PDG[imc];
0138       double px_ = mcpart_psx[imc];
0139       double py_ = mcpart_psy[imc];
0140       double pz_ = mcpart_psz[imc];
0141       double m_ = mcpart_m[imc];
0142       double e_ = sqrt(px_*px_+py_*py_+pz_*pz_+m_*m_);
0143
0144       // Add to mcpart
0145       Particles part;
0146
0147       part.pid=pid_;
0148       // part.charge = // TODO; not used yet
0149       part.mcID=imc;
0150       part.vecPart.SetPxPyPzE(px_,py_,pz_,e_);
0151       mcpart.push_back(part);
0152
0153       int igen=-1;
0154       if(auto search = genidmap.find(pz_); search != genidmap.end())
0155     igen=search->second; //index of the GeneratedParticle
0156
0157       mcidmap.insert({imc,igen});
0158
0159     }
0160
0161     /*
0162       ReconstructedParticles loop
0163       - Add all particles to the std::vector<> of particles
0164       - Look up associated MC particle
0165     */
0166
0167
0168
0169     for(int irec=0; irec < recparts_PDG.GetSize(); irec++){
0170
0171       int pid_ = recparts_PDG[irec];
0172       double px_ = recparts_p_x[irec];
0173       double py_ = recparts_p_y[irec];
0174       double pz_ = recparts_p_z[irec];
0175       double e_ = recparts_e[irec];
0176
0177       // Add to recpart
0178       Particles part;
0179
0180       part.pid=pid_;
0181       // part.charge = // TODO; not used yet
0182       part.vecPart.SetPxPyPzE(px_,py_,pz_,e_);
0183
0184       double m_ = part.vecPart.M();
0185
0186       /*
0187     Read through Associations to match particles
0188     By default, we assume no association, so mcID --> -1
0189     assoc_recID --> irec (index of the RecoParticle)
0190     assoc_simID --> imc (index of the MCParticle)
0191       */
0192
0193
0194       part.mcID=-1;
0195       for(int iassoc = 0 ; iassoc < assoc_simID.GetSize() ; iassoc++){
0196     int idx_recID = assoc_recID[iassoc];
0197     int idx_simID = assoc_simID[iassoc];
0198     if(irec==idx_recID){ // This track has an association
0199       part.mcID=idx_simID;
0200       break; // Only one association per particle
0201     }
0202       }
0203
0204       recpart.push_back(part);
0205       recidmap.insert({irec,part.mcID});
0206
0207     }
0208
0209
0210
0211     /*
0212       With the GeneratedParticles, MCParticles, and ReconstructedParticles filled,
0213       we can begin to search for the beam particles and hadronic final state (HFS)
0214       This is done for both the Truth and Reconstructed Particles
0215     */
0216
0217     /*
0218       Loop over MCParticles
0219     */
0220
0221     for(auto mcpart_: mcpart){
0222
0223       int imc = mcpart_.mcID;
0224       /* Beam particles have a MCParticles.generatorStatus of 4 */
0225       int genStat_ = mcpart_genStat[imc];
0226       if(mcpart_.pid==11 && genStat_ == 4){
0227     foundBeamElectron=true;
0228     kinTrue->vecEleBeam = mcpart_.vecPart;
0229       }
0230       else if(mcpart_.pid==2212 && genStat_ == 4){
0231     foundBeamIon=true;
0232     kinTrue->vecIonBeam = mcpart_.vecPart;
0233       }
0234       else if(genStat_==4){
0235     cout << "Warning...unknown beam particle with generatorStatus == 4 found...Continuing..." << endl;
0236       }
0237
0238       /* Assume the scattered electron is the pid==11 final state particle with the most energy */
0239       if(mcpart_.pid==11 && genStat_ == 1 && mcpart_.vecPart.P() > maxP)
0240     {
0241       maxP=mcpart_.vecPart.P();
0242       kinTrue->vecElectron = mcpart_.vecPart;
0243       genEleID = mcpart_.mcID;
0244     }
0245
0246       /*
0247     Only append MCParticles to the HFS if they are matched with a GeneratedParticle
0248        */
0249
0250       else if(genStat_ == 1 && mcidmap[mcpart_.mcID]>-1){
0251     kinTrue->AddToHFS(mcpart_.vecPart);
0252       }
0253
0254     }
0255
0256     //check beam finding
0257     if(!foundBeamElectron || !foundBeamIon) { numNoBeam++; continue;};
0258
0259     /*
0260       Loop over RecoParticles
0261     */
0262
0263     int irec = 0;
0264     bool recEleFound=false;
0265     for(auto recpart_ : recpart){
0266       // If there is a matching MCParticle
0267       if(recidmap[irec]!=-1){
0268     // If the recidmap is linked to the genEleID (generated scattered electron), identify this reco particle as the electron
0269     if(recidmap[irec]==genEleID){
0270       recEleFound=true;
0271       kin->vecElectron= recpart_.vecPart;
0272     }
0273     // Add the final state particle to the HFS
0274     kin->AddToHFS(recpart_.vecPart);
0275
0276     // Add reconstructed particle and true info to HFSTree
0277     if( writeHFSTree ){
0278       kin->AddToHFSTree(recpart_.vecPart, pid);
0279     }
0280       }
0281       irec++; // Increment to next particle
0282     }
0283
0284     // Skip event if the reco scattered electron was missing
0285     if(recEleFound==false){
0286       numNoEle++;
0287       continue;
0288     }
0289     else
0290       numEle++;
0291
0292     // subtrct electron from hadronic final state variables
0293     kin->SubtractElectronFromHFS();
0294     kinTrue->SubtractElectronFromHFS();
0295
0296     // skip the event if there are no reconstructed particles (other than the
0297     // electron), otherwise hadronic recon methods will fail
0298     if(kin->countHadrons == 0) {
0299       numNoHadrons++;
0300       continue;
0301     };
0302
0303     // calculate DIS kinematics
0304     if(!(kin->CalculateDIS(reconMethod))) continue; // reconstructed
0305     if(!(kinTrue->CalculateDIS(reconMethod))) continue; // generated (truth)
0306
0307     // Get the weight for this event's Q2
0308     auto Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]);
0309
0310     // fill inclusive histograms, if only `inclusive` is included in output
0311     // (otherwise they will be filled in track and jet loops)
0312     if(includeOutputSet["inclusive_only"]) {
0313       auto wInclusive = Q2weightFactor * weightInclusive->GetWeight(*kinTrue);
0314       wInclusiveTotal += wInclusive;
0315       FillHistosInclusive(wInclusive);
0316     }
0317
0318     /*
0319       Loop again over the reconstructed particles
0320       Calculate Hadron Kinematics
0321       Fill output data structures (Histos)
0322     */
0323
0324     int ipart = 0;
0325     for(auto part : recpart){
0326
0327       int pid_ = part.pid;
0328       int mcid_ = part.mcID;
0329
0330       // final state cut
0331       // - check PID, to see if it's a final state we're interested in for
0332       //   histograms; if not, proceed to next track
0333       auto kv = PIDtoFinalState.find(pid_);
0334       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
0335       if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
0336
0337       // calculate reconstructed hadron kinematics
0338       kin->vecHadron = part.vecPart;
0339       kin->CalculateHadronKinematics();
0340
0341       // add selected single hadron FS to HFS tree
0342       if( writeHFSTree ){
0343     kin->AddTrackToHFSTree(part.vecPart, part.pid);
0344       }
0345
0346
0347       // find the matching truth hadron using mcID, and calculate its kinematics
0348       if(mcid_ >= 0) {
0349     for(auto imc : mcpart) {
0350       if(mcid_ == imc.mcID) {
0351         kinTrue->vecHadron = imc.vecPart;
0352         // add tracks of interest for kinematic studies to HFSTree
0353         if( writeHFSTree ){
0354           kinTrue->AddTrackToHFSTree(imc.vecPart, imc.pid);
0355         }
0356         break;
0357       }
0358     }
0359       }
0360
0361       kinTrue->CalculateHadronKinematics();
0362
0363       if(includeOutputSet["1h"]) {
0364         // fill single-hadron histograms in activated bins
0365     auto wTrack = Q2weightFactor * weightTrack->GetWeight(*kinTrue);
0366         wTrackTotal += wTrack;
0367         FillHistos1h(wTrack);
0368         FillHistosInclusive(wTrack);
0369
0370     // fill simple tree
0371     // - not binned
0372     // - `IsActiveEvent()` is only true if at least one bin gets filled for this track
0373     if( writeSidisTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
0374       }
0375     } //hadron loop
0376
0377     if( writeHFSTree && kin->nHFS > 0) HFST->FillTree(Q2weightFactor);
0378
0379     /*
0380       Loop again over the reconstructed particles
0381       Fill output data structures (ParticleTree)
0382     */
0383
0384     // fill particle tree
0385     if ( writeParticleTree ) {
0386       int ipart = 0;
0387       for( auto part: recpart ){
0388       Particles mcpart_;
0389       int mcpart_idx=recidmap[ipart];     // Map idx to the matched MCParticle
0390       int genStat_ = -1;                    // Default Generator Status of MCParticle is -1 (no match)
0391       if(mcpart_idx>-1){                    // RecoParticle has an MCParticle match
0392         mcpart_ = mcpart.at(mcpart_idx);        // Get MCParticle
0393         int imc = mcpart_.mcID;
0394         genStat_ = mcpart_genStat[imc];         // Get Generator status of MCParticle
0395         if(imc==genEleID)                       // If MCParticle was scattered electron, set status to 2
0396           genStat_=2;
0397       }
0398       PT->FillTree(part.vecPart,      // Fill Tree
0399                mcpart_.vecPart,
0400                part.pid,
0401                genStat_);
0402       } // particle loop
0403       ipart++;
0404     }
0405
0406   } while(tr.Next());
0407
0408
0409   // finish execution
0410   Finish();
0411
0412   // final printout
0413   cout << "Total number of scattered electrons found: " << numEle << endl;
0414   if(numNoEle>0)
0415     cerr << "WARNING: skipped " << numNoEle << " events which had no reconstructed scattered electron" << endl;
0416   if(numNoHadrons>0)
0417     cerr << "WARNING: skipped " << numNoHadrons << " events which had no reconstructed hadrons" << endl;
0418   if(numNoBeam>0)
0419     cerr << "WARNING: skipped " << numNoBeam << " events which had no beam particles" << endl;
0420   if(numProxMatched>0)
0421     cerr << "WARNING: " << numProxMatched << " recon. particles were proximity matched to truth (when mcID match failed)" << endl;
0422 }