epic-analysis/src/AnalysisDelphes.cxx

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Christopher Dilks, Connor Pecar, Duane Byer, Matthew McEneaney, Brian Page
0003
0004 #include "AnalysisDelphes.h"
0005
0006 ClassImp(AnalysisDelphes)
0007
0008 using std::map;
0009 using std::vector;
0010 using std::cout;
0011 using std::cerr;
0012 using std::endl;
0013
0014 // constructor
0015 AnalysisDelphes::AnalysisDelphes(TString configFileName_, TString outfilePrefix_) :
0016   Analysis(configFileName_, outfilePrefix_)
0017 {
0018   // delphes-specific settings defaults
0019   /* ... none defined yet ... */
0020 };
0021
0022 // destructor
0023 AnalysisDelphes::~AnalysisDelphes() { };
0024
0025
0026 //=============================================
0027 // perform the analysis
0028 //=============================================
0029 void AnalysisDelphes::Execute() {
0030
0031   // setup
0032   Prepare();
0033
0034   // read delphes tree
0035   auto chain = std::make_unique<TChain>("Delphes");
0036   for(Int_t idx=0; idx<infiles.size(); ++idx) {
0037     for(std::size_t idxF=0; idxF<infiles[idx].size(); ++idxF) {
0038       // std::cout << "Adding " << infiles[idx][idxF] << std::endl;
0039       chain->Add(infiles[idx][idxF].c_str());
0040     }
0041   }
0042   chain->CanDeleteRefs();
0043   auto tr = std::make_unique<ExRootTreeReader>(chain.get());
0044   ENT = tr->GetEntries();
0045   if(maxEvents>0) ENT = std::min(maxEvents,ENT);
0046
0047   // branch iterators
0048   TObjArrayIter itTrack(tr->UseBranch("Track"));
0049   TObjArrayIter itElectron(tr->UseBranch("Electron"));
0050   TObjArrayIter itParticle(tr->UseBranch("Particle"));
0051   TObjArrayIter itEFlowTrack(tr->UseBranch("EFlowTrack"));
0052   TObjArrayIter itEFlowPhoton(tr->UseBranch("EFlowPhoton"));
0053   TObjArrayIter itEFlowNeutralHadron(tr->UseBranch("EFlowNeutralHadron"));
0054   TObjArrayIter itpfRICHTrack(tr->UseBranch("pfRICHTrack"));
0055   TObjArrayIter itDIRCepidTrack(tr->UseBranch("barrelDIRC_epidTrack"));
0056   TObjArrayIter itDIRChpidTrack(tr->UseBranch("barrelDIRC_hpidTrack"));
0057   TObjArrayIter itBTOFepidTrack(tr->UseBranch("BTOF_eTrack"));
0058   TObjArrayIter itBTOFhpidTrack(tr->UseBranch("BTOF_hTrack"));
0059   TObjArrayIter itdualRICHagTrack(tr->UseBranch("dualRICHagTrack"));
0060   TObjArrayIter itdualRICHcfTrack(tr->UseBranch("dualRICHcfTrack"));
0061
0062   // calculate Q2 weights
0063   CalculateEventQ2Weights();
0064
0065   // event loop =========================================================
0066   cout << "begin event loop..." << endl;
0067   for(Long64_t e=0; e<ENT; e++) {
0068     if(e>0&&e%10000==0) cout << (Double_t)e/ENT*100 << "%" << endl;
0069     tr->ReadEntry(e);
0070
0071     // electron loop
0072     // - finds max-momentum electron
0073     itElectron.Reset();
0074     maxEleP = 0;
0075     while(Electron *ele = (Electron*) itElectron()) {
0076       eleP = ele->PT * TMath::CosH(ele->Eta);
0077       if(eleP>maxEleP) {
0078         maxEleP = eleP;
0079         kin->vecElectron.SetPtEtaPhiM(
0080             ele->PT,
0081             ele->Eta,
0082             ele->Phi,
0083             Kinematics::ElectronMass()
0084             );
0085       };
0086     };
0087     if(maxEleP<0.001) continue; // no scattered electron found
0088
0089     // - repeat for truth electron
0090     itParticle.Reset();
0091     maxElePtrue = 0;
0092     bool found_elec = false;
0093     bool found_ion = false;
0094     while(GenParticle *part = (GenParticle*) itParticle()){
0095       if(part->PID == 11 && part->Status == 1){
0096         elePtrue = part->PT * TMath::CosH(part->Eta);
0097         if(elePtrue > maxElePtrue){
0098           maxElePtrue = elePtrue;
0099           kinTrue->vecElectron.SetPtEtaPhiM(
0100               part->PT,
0101               part->Eta,
0102               part->Phi,
0103               part->Mass
0104               );
0105         };
0106       };
0107       if(part->PID == 11 && part->Status == 4){
0108         if(!found_elec){
0109           found_elec = true;
0110           kinTrue->vecEleBeam.SetPtEtaPhiM(
0111               part->PT,
0112               part->Eta,
0113               part->Phi,
0114               part->Mass
0115               );
0116         }else{
0117           ErrorPrint("ERROR: Found two beam electrons in one event");
0118         };
0119       };
0120       if(part->PID != 11 && part->Status == 4){
0121         if(!found_ion){
0122           found_ion = true;
0123           kinTrue->vecIonBeam.SetPtEtaPhiM(
0124               part->PT,
0125               part->Eta,
0126               part->Phi,
0127               part->Mass
0128               );
0129         }else{
0130           ErrorPrint("ERROR: Found two beam ions in one event");
0131         };
0132       };
0133     };
0134     if(!found_elec){
0135       ErrorPrint("ERROR: Didn't find beam electron in event");
0136     };
0137     if(!found_ion){
0138       ErrorPrint("ERROR: Didn't find beam ion in event");
0139     }
0140
0141     // get hadronic final state variables
0142     kin->GetHFS(
0143         itTrack,
0144         itEFlowTrack,
0145         itEFlowPhoton,
0146         itEFlowNeutralHadron,
0147         itpfRICHTrack,
0148         itDIRCepidTrack, itDIRChpidTrack,
0149         itBTOFepidTrack, itBTOFhpidTrack,
0150         itdualRICHagTrack, itdualRICHcfTrack
0151         );
0152     kinTrue->GetTrueHFS(itParticle);
0153
0154     // calculate DIS kinematics
0155     if(!(kin->CalculateDIS(reconMethod))) continue; // reconstructed
0156     if(!(kinTrue->CalculateDIS(reconMethod))) continue; // generated (truth)
0157
0158     // Get the weight for this event's Q2
0159     auto Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]);
0160
0161     // fill inclusive histograms, if only `inclusive` is included in output
0162     // (otherwise they will be filled in track and jet loops)
0163     if(includeOutputSet["inclusive_only"]) {
0164       auto wInclusive = Q2weightFactor * weightInclusive->GetWeight(*kinTrue);
0165       wInclusiveTotal += wInclusive;
0166       FillHistosInclusive(wInclusive);
0167     }
0168
0169     // track loop - - - - - - - - - - - - - - - - - - - - - - - - - - - -
0170     itTrack.Reset();
0171     while(Track *trk = (Track*) itTrack()) {
0172       //cout << e << " " << trk->PID << endl;
0173
0174       // final state cut
0175       // - check PID, to see if it's a final state we're interested in for
0176       //   histograms; if not, proceed to next track
0177       // pid = trk->PID; //NOTE: trk->PID is currently not smeared so it just returns the truth-level PID
0178       pid = kin->GetTrackPID( // get smeared PID
0179           trk,
0180           itpfRICHTrack,
0181           itDIRCepidTrack, itDIRChpidTrack,
0182           itBTOFepidTrack, itBTOFhpidTrack,
0183           itdualRICHagTrack, itdualRICHcfTrack
0184           );
0185       auto kv = PIDtoFinalState.find(pid);
0186       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
0187       if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
0188
0189       // get parent particle, to check if pion is from vector meson
0190       GenParticle *trkParticle = (GenParticle*)trk->Particle.GetObject();
0191       TObjArray *brParticle = (TObjArray*)itParticle.GetCollection();
0192       GenParticle *parentParticle = (GenParticle*)brParticle->At(trkParticle->M1);
0193       int parentPID = (parentParticle->PID); // TODO: this is not used yet...
0194
0195       // calculate hadron kinematics
0196       kin->hadPID = pid;
0197       kin->vecHadron.SetPtEtaPhiM(
0198           trk->PT,
0199           trk->Eta,
0200           trk->Phi,
0201           trk->Mass /* TODO: do we use track mass here ?? */
0202           );
0203
0204       GenParticle* trkPart = (GenParticle*)trk->Particle.GetObject();
0205       kinTrue->hadPID = pid;
0206       kinTrue->vecHadron.SetPtEtaPhiM(
0207           trkPart->PT,
0208           trkPart->Eta,
0209           trkPart->Phi,
0210           trkPart->Mass /* TODO: do we use track mass here ?? */
0211           );
0212
0213       kin->CalculateHadronKinematics();
0214       kinTrue->CalculateHadronKinematics();
0215
0216       if( writeHFSTree ){
0217     kin->AddTrackToHFSTree(kin->vecHadron,kin->hadPID);
0218     kinTrue->AddTrackToHFSTree(kinTrue->vecHadron,kinTrue->hadPID);
0219       }
0220
0221       // asymmetry injection
0222       //kin->InjectFakeAsymmetry(); // sets tSpin, based on reconstructed kinematics
0223       //kinTrue->InjectFakeAsymmetry(); // sets tSpin, based on generated kinematics
0224       //kin->tSpin = kinTrue->tSpin; // copy to "reconstructed" tSpin
0225
0226       if(includeOutputSet["1h"]) {
0227         // fill single-hadron histograms in activated bins
0228         auto wTrack = Q2weightFactor * weightTrack->GetWeight(*kinTrue);
0229         wTrackTotal += wTrack;
0230         FillHistos1h(wTrack);
0231         FillHistosInclusive(wTrack);
0232
0233         // fill simple tree
0234         // - not binned
0235         // - `IsActiveEvent()` is only true if at least one bin gets filled for this track
0236         if( writeSidisTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
0237       }
0238
0239       // tests
0240       //kin->ValidateHeadOnFrame();
0241
0242     }; // end track loop
0243
0244     if( writeHFSTree ) HFST->FillTree(Q2weightFactor);
0245     // jet loop - - - - - - - - - - - - - - - - - - - - - - - - - - - -
0246     if(includeOutputSet["jets"]) {
0247
0248       // get vector of jets
0249       // TODO: should this have an option for clustering method?
0250       //kin->GetJets(itEFlowTrack, itEFlowPhoton, itEFlowNeutralHadron, itParticle);
0251       kinJet->GetJets(itEFlowTrack, itEFlowPhoton, itEFlowNeutralHadron, itParticle, jetAlg, jetRad, jetMin);
0252
0253       finalStateID = "jet";
0254
0255 #ifdef INCLUDE_CENTAURO
0256       if(useBreitJets) kinJet->GetBreitFrameJets(itEFlowTrack, itEFlowPhoton, itEFlowNeutralHadron, itParticle);
0257 #endif
0258
0259       auto wJet = Q2weightFactor * weightJet->GetWeight(*kinJetTrue); // TODO: should we separate weights for breit and non-breit jets?
0260       wJetTotal += wJet;
0261
0262       Int_t nJets;
0263       if(useBreitJets) nJets = kinJet->breitJetsRec.size();
0264       else      nJets = kinJet->jetsRec.size();
0265
0266       for(int i = 0; i < kinJet->jetsRec.size(); i++){
0267
0268         if(useBreitJets) {
0269 #ifdef INCLUDE_CENTAURO
0270           jet = kinJet->breitJetsRec[i];
0271           kinJet->CalculateBreitJetKinematics(jet);
0272 #endif
0273         } else {
0274           jet = kinJet->jetsRec[i];
0275           kinJet->CalculateJetKinematics(jet);
0276
0277       // Match Reco Jet to Nearest Truth Jet (Specify DR matching limit between jets)
0278       kinJet->CalculateJetResolution(jetMatchDR);
0279         };
0280
0281         // fill jet histograms in activated bins
0282         FillHistosJets(wJet);
0283         FillHistosInclusive(wJet);
0284
0285       };
0286     }; // end jet loop
0287
0288   };
0289   cout << "end event loop" << endl;
0290   // event loop end =========================================================
0291
0292
0293   // finish execution
0294   Finish();
0295   //cout << "DEBUG PID in HFS: nSmeared=" << kin->countPIDsmeared << "  nNotSmeared=" << kin->countPIDtrue << endl;
0296 };