epic-analysis/src/AnalysisAthena.cxx

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Sanghwa Park, Christopher Dilks
0003
0004 #include "AnalysisAthena.h"
0005
0006 using std::map;
0007 using std::vector;
0008 using std::cout;
0009 using std::cerr;
0010 using std::endl;
0011
0012 // constructor
0013 AnalysisAthena::AnalysisAthena(TString configFileName_, TString outfilePrefix_) :
0014   Analysis(configFileName_, outfilePrefix_)
0015 { };
0016
0017 // destructor
0018 AnalysisAthena::~AnalysisAthena() { };
0019
0020
0021 //=============================================
0022 // perform the analysis
0023 //=============================================
0024 void AnalysisAthena::Execute()
0025 {
0026   // setup
0027   Prepare();
0028
0029   // read dd4hep tree
0030   auto chain = std::make_unique<TChain>("events");
0031   for(Int_t idx=0; idx<infiles.size(); ++idx) {
0032     for(std::size_t idxF=0; idxF<infiles[idx].size(); ++idxF) {
0033       // std::cout << "Adding " << infiles[idx][idxF] << std::endl;
0034       chain->Add(infiles[idx][idxF].c_str());
0035     }
0036   }
0037   chain->CanDeleteRefs();
0038
0039   TTreeReader tr(chain.get());
0040
0041   // Truth
0042   TTreeReaderArray<Int_t>    mcparticles_ID(tr,        "mcparticles.ID");
0043   TTreeReaderArray<Int_t>    mcparticles_pdgID(tr,     "mcparticles.pdgID");
0044   TTreeReaderArray<Double_t> mcparticles_psx(tr,       "mcparticles.ps.x");
0045   TTreeReaderArray<Double_t> mcparticles_psy(tr,       "mcparticles.ps.y");
0046   TTreeReaderArray<Double_t> mcparticles_psz(tr,       "mcparticles.ps.z");
0047   TTreeReaderArray<Int_t>    mcparticles_status(tr,    "mcparticles.status");
0048   TTreeReaderArray<Int_t>    mcparticles_genStatus(tr, "mcparticles.genStatus");
0049   TTreeReaderArray<Double_t> mcparticles_mass(tr,      "mcparticles.mass");
0050
0051   // Reco
0052   TTreeReaderArray<Int_t> ReconstructedParticles_pid(tr,   "ReconstructedParticles.pid");
0053   TTreeReaderArray<float> ReconstructedParticles_energy(tr,  "ReconstructedParticles.energy");
0054   TTreeReaderArray<float> ReconstructedParticles_p_x(tr,     "ReconstructedParticles.p.x");
0055   TTreeReaderArray<float> ReconstructedParticles_p_y(tr,     "ReconstructedParticles.p.y");
0056   TTreeReaderArray<float> ReconstructedParticles_p_z(tr,     "ReconstructedParticles.p.z");
0057   TTreeReaderArray<float> ReconstructedParticles_p(tr,       "ReconstructedParticles.momentum");
0058   TTreeReaderArray<float> ReconstructedParticles_th(tr,      "ReconstructedParticles.direction.theta");
0059   TTreeReaderArray<float> ReconstructedParticles_phi(tr,     "ReconstructedParticles.direction.phi");
0060   TTreeReaderArray<float> ReconstructedParticles_mass(tr,    "ReconstructedParticles.mass");
0061   TTreeReaderArray<short> ReconstructedParticles_charge(tr,  "ReconstructedParticles.charge");
0062   TTreeReaderArray<int>   ReconstructedParticles_mcID(tr,    "ReconstructedParticles.mcID.value");
0063
0064   // calculate Q2 weights
0065   CalculateEventQ2Weights();
0066
0067   // counters
0068   Long64_t numNoBeam, numEle, numNoEle, numNoHadrons, numProxMatched;
0069   numNoBeam = numEle = numNoEle = numNoHadrons = numProxMatched = 0;
0070
0071   // event loop =========================================================
0072   cout << "begin event loop..." << endl;
0073   tr.SetEntriesRange(1,maxEvents);
0074   do {
0075     if(tr.GetCurrentEntry()%10000==0) cout << tr.GetCurrentEntry() << " events..." << endl;
0076
0077     // resets
0078     kin->ResetHFS();
0079     kinTrue->ResetHFS();
0080
0081     // generated truth loop
0082     /* - add truth particle to `mcpart`
0083      * - add to hadronic final state sums (momentum, sigma, etc.)
0084      * - find scattered electron
0085      * - find beam particles
0086      */
0087     std::vector<Particles> mcpart;
0088     double maxP = 0;
0089     int genEleID = -1;
0090     bool foundBeamElectron = false;
0091     bool foundBeamIon = false;
0092     for(int imc=0; imc<mcparticles_pdgID.GetSize(); imc++) {
0093
0094       int pid_ = mcparticles_pdgID[imc];
0095       int genStatus_ = mcparticles_genStatus[imc]; // genStatus 4: beam particle,  1: final state
0096       double px_ = mcparticles_psx[imc];
0097       double py_ = mcparticles_psy[imc];
0098       double pz_ = mcparticles_psz[imc];
0099       double mass_ = mcparticles_mass[imc]; // in GeV
0100       double p_ = sqrt(pow(mcparticles_psx[imc],2) + pow(mcparticles_psy[imc],2) + pow(mcparticles_psz[imc],2));
0101
0102       if(genStatus_ == 1) { // final state
0103
0104         // add to `mcpart`
0105         Particles part;
0106         part.pid = pid_;
0107         part.vecPart.SetPxPyPzE(px_, py_, pz_, sqrt(p_*p_ + mass_*mass_));
0108         part.mcID = mcparticles_ID[imc];
0109         mcpart.push_back(part);
0110
0111         // add to hadronic final state sums
0112         kinTrue->AddToHFS(part.vecPart);
0113
0114         // identify scattered electron by max momentum
0115         if(pid_ == 11) {
0116           if(p_ > maxP) {
0117             maxP = p_;
0118             kinTrue->vecElectron.SetPxPyPzE(px_, py_, pz_, sqrt(p_*p_ + mass_*mass_));
0119             genEleID = mcparticles_ID[imc];
0120           }
0121         }
0122       }
0123
0124       else if(genStatus_ == 4) { // beam particles
0125         if(pid_ == 11) { // electron beam
0126           if(!foundBeamElectron) {
0127             foundBeamElectron = true;
0128             kinTrue->vecEleBeam.SetPxPyPzE(px_, py_, pz_, sqrt(p_*p_ + mass_*mass_));
0129           }
0130           else { ErrorPrint("ERROR: Found two beam electrons in one event"); }
0131         }
0132         else { // ion beam
0133           if(!foundBeamIon) {
0134             foundBeamIon = true;
0135             kinTrue->vecIonBeam.SetPxPyPzE(px_, py_, pz_, sqrt(p_*p_ + mass_*mass_));
0136           }
0137           else { ErrorPrint("ERROR: Found two beam ions in one event"); }
0138         }
0139       }
0140     } // end truth loop
0141
0142     // check beam finding
0143     if(!foundBeamElectron || !foundBeamIon) { numNoBeam++; continue; };
0144
0145
0146     // reconstructed particles loop
0147     /* - add reconstructed particle to `recopart`
0148      * - find the scattered electron
0149      *
0150      */
0151     std::vector<Particles> recopart;
0152     int recEleFound = 0;
0153     for(int ireco=0; ireco<ReconstructedParticles_pid.GetSize(); ireco++) {
0154
0155       int pid_ = ReconstructedParticles_pid[ireco];
0156       if(pid_ == 0) continue; // pid==0: reconstructed tracks with no matching truth pid
0157
0158       // add reconstructed particle `part` to `recopart`
0159       Particles part;
0160       part.pid = pid_;
0161       part.mcID = ReconstructedParticles_mcID[ireco];
0162       part.charge = ReconstructedParticles_charge[ireco];
0163       double reco_E = ReconstructedParticles_energy[ireco];
0164       double reco_px = ReconstructedParticles_p_x[ireco];
0165       double reco_py = ReconstructedParticles_p_y[ireco];
0166       double reco_pz = ReconstructedParticles_p_z[ireco];
0167       double reco_mass = ReconstructedParticles_mass[ireco];
0168       double reco_p = sqrt(reco_px*reco_px + reco_py*reco_py + reco_pz*reco_pz);
0169       part.vecPart.SetPxPyPzE(reco_px, reco_py, reco_pz, sqrt(reco_p*reco_p + reco_mass*reco_mass));
0170
0171       // add to `recopart` and hadronic final state sums only if there is a matching truth particle
0172       if(part.mcID > 0) {
0173         for(auto imc : mcpart) {
0174           if(part.mcID == imc.mcID) {
0175             recopart.push_back(part);
0176             kin->AddToHFS(part.vecPart);
0177         if( writeHFSTree ){
0178           kin->AddToHFSTree(part.vecPart, part.pid);
0179         }
0180             break;
0181           }
0182         }
0183       }
0184
0185       // find scattered electron, by matching to truth // TODO: not realistic... is there an upstream electron finder?
0186       if(pid_ == 11 && part.mcID == genEleID) {
0187         recEleFound++;
0188         kin->vecElectron.SetPxPyPzE(reco_px, reco_py, reco_pz, sqrt(reco_p*reco_p + reco_mass*reco_mass));
0189       }
0190
0191     } // end reco loop
0192
0193     // skip the event if the scattered electron is not found and we need it
0194     if(recEleFound < 1) {
0195       numNoEle++;
0196       continue; // TODO: only need to skip if we are using a recon method that needs it (`if reconMethod==...`)
0197     }
0198     else if(recEleFound>1) ErrorPrint("WARNING: found more than 1 reconstructed scattered electron in an event");
0199     else numEle++;
0200
0201     // subtract electron from hadronic final state variables
0202     kin->SubtractElectronFromHFS();
0203     kinTrue->SubtractElectronFromHFS();
0204
0205     // skip the event if there are no reconstructed particles (other than the
0206     // electron), otherwise hadronic recon methods will fail
0207     if(kin->countHadrons == 0) {
0208       numNoHadrons++;
0209       continue;
0210     };
0211
0212     // calculate DIS kinematics
0213     if(!(kin->CalculateDIS(reconMethod))) continue; // reconstructed
0214     if(!(kinTrue->CalculateDIS(reconMethod))) continue; // generated (truth)
0215
0216     // Get the weight for this event's Q2
0217     auto Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]);
0218
0219     // fill inclusive histograms, if only `inclusive` is included in output
0220     // (otherwise they will be filled in track and jet loops)
0221     if(includeOutputSet["inclusive_only"]) {
0222       auto wInclusive = Q2weightFactor * weightInclusive->GetWeight(*kinTrue);
0223       wInclusiveTotal += wInclusive;
0224       FillHistosInclusive(wInclusive);
0225     }
0226
0227     // loop over reconstructed particles again
0228     /* - calculate hadron kinematics
0229      * - fill output data structures (Histos, SidisTree, etc.)
0230      */
0231     for(auto part : recopart) {
0232       int pid_ = part.pid;
0233       int mcid_ = part.mcID;
0234
0235       // final state cut
0236       // - check PID, to see if it's a final state we're interested in for
0237       //   histograms; if not, proceed to next track
0238       auto kv = PIDtoFinalState.find(pid_);
0239       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
0240       if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
0241
0242       // calculate reconstructed hadron kinematics
0243       kin->vecHadron = part.vecPart;
0244       kin->CalculateHadronKinematics();
0245
0246       // add selected single hadron FS to HFS tree
0247       if( writeHFSTree ){
0248         kin->AddTrackToHFSTree(part.vecPart, part.pid);
0249       }
0250
0251       // find the matching truth hadron using mcID, and calculate its kinematics
0252       if(mcid_ > 0) {
0253         for(auto imc : mcpart) {
0254           if(mcid_ == imc.mcID) {
0255             kinTrue->vecHadron = imc.vecPart;
0256         // add tracks of interest for kinematic studies to HFSTree
0257             if( writeHFSTree ){
0258               kinTrue->AddTrackToHFSTree(imc.vecPart, imc.pid);
0259             }
0260             break;
0261           }
0262         }
0263       }
0264       /* // deprecated, since existence of truth match is checked earlier; in practice prox matching was never called
0265       else {
0266         // give it another shot: proximity matching
0267         double mineta = 4.0;
0268         numProxMatched++;
0269         for(int imc=0; imc<(int)mcpart.size(); imc++) {
0270           if(pid_ == mcpart[imc].pid) {
0271             double deta = abs(kin->vecHadron.Eta() - mcpart[imc].vecPart.Eta());
0272             if(deta < mineta) {
0273               mineta = deta;
0274               kinTrue->vecHadron = mcpart[imc].vecPart;
0275             }
0276           }
0277         }
0278       }
0279       */
0280       kinTrue->CalculateHadronKinematics();
0281
0282       // asymmetry injection
0283       // kin->InjectFakeAsymmetry(); // sets tSpin, based on reconstructed kinematics
0284       // kinTrue->InjectFakeAsymmetry(); // sets tSpin, based on generated kinematics
0285       // kin->tSpin = kinTrue->tSpin; // copy to "reconstructed" tSpin
0286
0287       // weighting
0288       auto wTrack = Q2weightFactor * weightTrack->GetWeight(*kinTrue);
0289       wTrackTotal += wTrack;
0290
0291       if(includeOutputSet["1h"]) {
0292         // fill single-hadron histograms in activated bins
0293         FillHistos1h(wTrack);
0294         FillHistosInclusive(wTrack);
0295
0296         // fill simple tree
0297         // - not binned
0298         // - `IsActiveEvent()` is only true if at least one bin gets filled for this track
0299         if( writeSidisTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
0300       }
0301
0302     }//hadron loop
0303
0304     // fill HFSTree for each event
0305     if( writeHFSTree && kin->nHFS > 0) HFST->FillTree(Q2weightFactor);
0306
0307   } while(tr.Next()); // tree reader loop
0308   cout << "end event loop" << endl;
0309   // event loop end =========================================================
0310
0311   // finish execution
0312   Finish();
0313
0314   // final printout
0315   cout << "Total number of scattered electrons found: " << numEle << endl;
0316   if(numNoEle>0)
0317     cerr << "WARNING: skipped " << numNoEle << " events which had no reconstructed scattered electron" << endl;
0318   if(numNoHadrons>0)
0319     cerr << "WARNING: skipped " << numNoHadrons << " events which had no reconstructed hadrons" << endl;
0320   if(numNoBeam>0)
0321     cerr << "WARNING: skipped " << numNoBeam << " events which had no beam particles" << endl;
0322   if(numProxMatched>0)
0323     cerr << "WARNING: " << numProxMatched << " recon. particles were proximity matched to truth (when mcID match failed)" << endl;
0324
0325 }