epic-analysis/src/AnalysisEcce.cxx

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Ralf Seidl, Christopher Dilks, Sanghwa Park
0003
0004 #include "AnalysisEcce.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 AnalysisEcce::AnalysisEcce(TString configFileName_, TString outfilePrefix_) :
0014   Analysis(configFileName_, outfilePrefix_),
0015   trackSource(0) // default track source is "all tracks"
0016 { };
0017
0018 // destructor
0019 AnalysisEcce::~AnalysisEcce() { };
0020
0021
0022 //=============================================
0023 // perform the analysis
0024 //=============================================
0025 void AnalysisEcce::Execute()
0026 {
0027   // setup
0028   Prepare();
0029
0030   // read EventEvaluator tree
0031   auto chain = std::make_unique<TChain>("event_tree");
0032   for(Int_t idx=0; idx<infiles.size(); ++idx) {
0033     for(std::size_t idxF=0; idxF<infiles[idx].size(); ++idxF) {
0034       // std::cout << "Adding " << infiles[idx][idxF] << std::endl;
0035       chain->Add(infiles[idx][idxF].c_str());
0036     }
0037   }
0038   chain->CanDeleteRefs();
0039
0040   TTreeReader tr(chain.get());
0041
0042   TBranch        *b_clusters;   //!
0043   TBranch        *b_cluster_E;   //!
0044   TBranch        *b_cluster_Eta;   //!
0045   TBranch        *b_cluster_Phi;   //!
0046   TBranch        *b_cluster_NTower;   //!
0047   TBranch        *b_cluster_trueID;   //!
0048
0049   Int_t           cluster_N;
0050
0051
0052
0053   // Truth
0054
0055   TTreeReaderArray<Int_t> hepmcp_status(tr, "hepmcp_status");
0056   TTreeReaderArray<Int_t> hepmcp_PDG(tr,    "hepmcp_PDG");
0057   TTreeReaderArray<Float_t> hepmcp_E(tr,      "hepmcp_E");
0058   TTreeReaderArray<Float_t> hepmcp_psx(tr,    "hepmcp_px");
0059   TTreeReaderArray<Float_t> hepmcp_psy(tr,    "hepmcp_py");
0060   TTreeReaderArray<Float_t> hepmcp_psz(tr,    "hepmcp_pz");
0061
0062   TTreeReaderArray<Int_t> hepmcp_BCID(tr, "hepmcp_BCID");
0063   TTreeReaderArray<Int_t> hepmcp_m1(tr, "hepmcp_m1");
0064   TTreeReaderArray<Int_t> hepmcp_m2(tr, "hepmcp_m2");
0065
0066
0067   // All true particles (including secondaries, etc)
0068   TTreeReaderArray<Int_t> mcpart_ID(tr,        "mcpart_ID");
0069   TTreeReaderArray<Int_t> mcpart_ID_parent(tr, "mcpart_ID_parent");
0070   TTreeReaderArray<Int_t> mcpart_PDG(tr,       "mcpart_PDG");
0071   TTreeReaderArray<Float_t> mcpart_E(tr,         "mcpart_E");
0072   TTreeReaderArray<Float_t> mcpart_psx(tr,       "mcpart_px");
0073   TTreeReaderArray<Float_t> mcpart_psy(tr,       "mcpart_py");
0074   TTreeReaderArray<Float_t> mcpart_psz(tr,       "mcpart_pz");
0075   TTreeReaderArray<Int_t> mcpart_BCID(tr,      "mcpart_BCID");
0076
0077
0078   // Reco tracks
0079   TTreeReaderArray<Float_t> tracks_id(tr,  "tracks_ID"); // needs to be made an int eventually in actual EE code
0080   TTreeReaderArray<Float_t> tracks_p_x(tr, "tracks_px");
0081   TTreeReaderArray<Float_t> tracks_p_y(tr, "tracks_py");
0082   TTreeReaderArray<Float_t> tracks_p_z(tr,  "tracks_pz");
0083   TTreeReaderArray<Float_t> tracks_trueID(tr,  "tracks_trueID");
0084   TTreeReaderArray<UShort_t> tracks_source(tr,  "tracks_source");
0085    // TTreeReaderArray<Short_t> tracks_charge(tr,  "tracks_charge");
0086
0087
0088   //define the different calorimeter clusters
0089   std::vector<std::string> calos = {"FEMC","CEMC","EEMC","FHCAL","EHCAL","HCALOUT","HCALIN"};
0090
0091
0092
0093   // calculate Q2 weights
0094   CalculateEventQ2Weights();
0095
0096   // counters
0097   Long64_t numNoBeam, numEle, numNoEle, numNoHadrons, numProxMatched;
0098   numNoBeam = numEle = numNoEle = numNoHadrons = numProxMatched = 0;
0099
0100   // event loop =========================================================
0101   cout << "begin event loop..." << endl;
0102   tr.SetEntriesRange(1,maxEvents);
0103   do {
0104     if(tr.GetCurrentEntry()%10000==0) cout << tr.GetCurrentEntry() << " events..." << endl;
0105
0106     // resets
0107     kin->ResetHFS();
0108     kinTrue->ResetHFS();
0109
0110     // a few maps needed to get the associated info between tracks, true particles, etec
0111
0112     std::map <int,int> mcidmap;
0113     std::map <int,int> mcbcidmap;
0114     std::map <int,int> mcbcidmap2;
0115     std::map <int,int> trackmap; // mapping of all the tracks
0116
0117     for (int imc =0; imc < mcpart_ID.GetSize(); imc++){
0118       if (mcpart_E[imc]<0.1) continue;
0119       int id = (int)mcpart_ID[imc];
0120       int bcid = (int)mcpart_BCID[imc];
0121       mcidmap.insert({id,imc});
0122       mcbcidmap.insert({bcid,imc});
0123     }
0124
0125     for(int itrack=0; itrack<tracks_id.GetSize(); itrack++) {
0126       if(tracks_source[itrack]!=trackSource) continue;
0127       trackmap.insert({tracks_trueID[itrack],itrack});
0128     }
0129
0130
0131
0132
0133
0134     // generated truth loop
0135     /* - add truth particle to `mcpart`
0136      * - add to hadronic final state sums (momentum, sigma, etc.)
0137      * - find scattered electron
0138      * - find beam particles
0139      */
0140     std::vector<Particles> mcpart;
0141     double maxPt = 0;
0142     int genEleID = -1;
0143     bool foundBeamElectron = false;
0144     bool foundBeamIon = false;
0145     int genEleBCID = -1;
0146
0147     for(int imc=0; imc<hepmcp_PDG.GetSize(); imc++) {
0148
0149       int pid_ = hepmcp_PDG[imc];
0150
0151       int genStatus_ = hepmcp_status[imc]; // genStatus 4: beam particle,  1: final state
0152
0153       double px_ = hepmcp_psx[imc];
0154       double py_ = hepmcp_psy[imc];
0155       double pz_ = hepmcp_psz[imc];
0156       double e_  = hepmcp_E[imc];
0157
0158       double p_ = sqrt(pow(hepmcp_psx[imc],2) + pow(hepmcp_psy[imc],2) + pow(hepmcp_psz[imc],2));
0159       double pt_ = sqrt(pow(hepmcp_psx[imc],2) + pow(hepmcp_psy[imc],2));
0160       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.;
0161
0162       // add to `mcpart`
0163       Particles part;
0164
0165       //cout << genStatus_ << " " << pid_ << " " << part.mcID << " " << hepmcp_BCID[imc] << " " << hepmcp_m1[imc] << " " << hepmcp_m2[imc] << " "
0166       //      << px_ << " " << py_ << " " << pz_ << " " << e_ << endl;
0167
0168       if(genStatus_ == 1) { // final state
0169
0170     int imcpart = -1;// matched truthtrack
0171     auto search = mcbcidmap.find((int)(hepmcp_BCID[imc]));
0172     if (search != mcbcidmap.end()) {
0173       imcpart = search->second;
0174     }
0175
0176
0177     if (imcpart >-1){
0178       px_ = mcpart_psx[imcpart];
0179       py_ = mcpart_psy[imcpart];
0180       pz_ = mcpart_psz[imcpart];
0181       e_  = mcpart_E[imcpart];
0182       p_ = sqrt(pow(mcpart_psx[imcpart],2) + pow(mcpart_psy[imcpart],2) + pow(mcpart_psz[imcpart],2));
0183       pt_ = sqrt(pow(mcpart_psx[imcpart],2) + pow(mcpart_psy[imcpart],2));
0184       part.mcID = mcpart_ID[imcpart];
0185     }
0186       else
0187         part.mcID = -1;
0188     //cout << genStatus_ << " " << pid_ << " " << part.mcID << " " << hepmcp_BCID[imc] << " " << hepmcp_m1[imc] << " " << hepmcp_m2[imc] << " "
0189     //    << px_ << " " << py_ << " " << pz_ << " " << e_ << endl;
0190
0191       part.pid = pid_;
0192       part.vecPart.SetPxPyPzE(px_, py_, pz_, e_);
0193
0194       mcpart.push_back(part);
0195
0196       // add to hadronic final state sums
0197       kinTrue->AddToHFS(part.vecPart);
0198
0199       // identify scattered electron by max momentum
0200       if(pid_ == 11) {
0201         if(pt_ > maxPt) {
0202           maxPt = pt_;
0203           kinTrue->vecElectron.SetPxPyPzE(px_, py_, pz_, e_);
0204           genEleID = part.mcID; //mcpart_ID[imc];
0205           genEleBCID = hepmcp_BCID[imc];
0206           //          cout  << "\t\t\t found scattered electron  " << Form(" %6.2f %6.2f %6.2f %6.2f %6.2f  %5.3f %6.2f %6.2f id %3d\n",px_,py_,pz_, sqrt(p_*p_ + mass_*mass_),p_,mass_,hepmcp_E[imc],mcpart_E[imcpart],genEleID);
0207         }
0208       }
0209       }
0210
0211       else if(genStatus_ == 4) { // beam particles
0212         if(pid_ == 11) { // electron beam
0213           if(!foundBeamElectron) {
0214             foundBeamElectron = true;
0215             kinTrue->vecEleBeam.SetPxPyPzE(px_, py_, pz_, sqrt(p_*p_ + mass_*mass_));
0216         //      cout  << "\t\t\t found beam electron  " << Form(" %4.2f %4.2f %4.2f \n",px_,py_,pz_);
0217           }
0218           else { ErrorPrint("ERROR: Found two beam electrons in one event"); }
0219         }
0220         else { // ion beam
0221           if(!foundBeamIon) {
0222             foundBeamIon = true;
0223             kinTrue->vecIonBeam.SetPxPyPzE(px_, py_, pz_, sqrt(p_*p_ + mass_*mass_));
0224         //      cout  << "\t\t\t found beam ion  " << Form(" %4.2f %4.2f %4.2f \n",px_,py_,pz_);
0225           }
0226           else { ErrorPrint("ERROR: Found two beam ions in one event"); }
0227         }
0228       }
0229     } // end truth loop
0230
0231     // looking for radiated photons (mother particle == scattered electron)
0232     // requires that we already know scattered electron index
0233     for(int imc=0; imc<hepmcp_PDG.GetSize(); imc++) {
0234       int pid_ = hepmcp_PDG[imc];
0235       int genStatus_ = hepmcp_status[imc]; // genStatus 4: beam particle,  1: final state
0236
0237       double px_ = hepmcp_psx[imc];
0238       double py_ = hepmcp_psy[imc];
0239       double pz_ = hepmcp_psz[imc];
0240       double e_  = hepmcp_E[imc];
0241
0242       double p_ = sqrt(pow(hepmcp_psx[imc],2) + pow(hepmcp_psy[imc],2) + pow(hepmcp_psz[imc],2));
0243
0244       int mother = hepmcp_m1[imc];
0245       if(pid_ == 22 || pid_ == 23){
0246     if(genStatus_ == 1){
0247       if( mother == genEleBCID ){
0248         TLorentzVector FSRmom(px_, py_, pz_, e_);
0249         TLorentzVector eleCorr = kinTrue->vecElectron + FSRmom;
0250         // correcting true scattered electron with FSR
0251         kinTrue->vecElectron.SetPxPyPzE(eleCorr.Px(), eleCorr.Py(), eleCorr.Pz(), eleCorr.E());
0252       }
0253     }
0254       }
0255
0256
0257     }
0258     // check beam finding
0259     if(!foundBeamElectron || !foundBeamIon) { numNoBeam++; continue; };
0260
0261     // reconstructed particles loop
0262     /* - add reconstructed particle to `recopart`
0263      * - find the scattered electron
0264      *
0265      */
0266     std::vector<Particles> recopart;
0267     int recEleFound = 0;
0268     for(int ireco=0; ireco<tracks_id.GetSize(); ireco++) {
0269
0270       // skip track if not from the user-specified source `trackSource`
0271       if(tracks_source[ireco]!=trackSource) continue;
0272
0273       int pid_ = 0; //tracks_pid[ireco];
0274       double reco_mass = 0.;
0275       int imc = -1;// matched truthtrack
0276       auto search = mcidmap.find((int)(tracks_trueID[ireco]));
0277       if (search != mcidmap.end()) {
0278     imc = search->second;
0279     pid_ = (int)(mcpart_PDG[imc]);
0280       }
0281       // later also use the likelihoods instead for pid
0282       //      cout  << "\t\t track " << ireco << " PDG  " << pid_ << " mcpart imc " << imc << endl;
0283       if(pid_ == 0) continue; // pid==0: reconstructed tracks with no matching truth pid
0284
0285       // add reconstructed particle `part` to `recopart`
0286       Particles part;
0287       part.pid = pid_;
0288       part.mcID = tracks_trueID[ireco];
0289       //      part.charge = tracks_charge[ireco];
0290       part.charge = (pid_ == 211 || pid_ == 321 || pid_ == 2212 || pid_ == -11 || pid_ == -13)?1:(pid_ == -211 || pid_ == -321 || pid_ == -2212 || pid_ == 11 || pid_ == 13)?-1:0;
0291
0292       double reco_px = tracks_p_x[ireco];
0293       double reco_py = tracks_p_y[ireco];
0294       double reco_pz = tracks_p_z[ireco];
0295       reco_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.;
0296
0297       double reco_p = sqrt(reco_px*reco_px + reco_py*reco_py + reco_pz*reco_pz);
0298       double reco_E = sqrt(reco_p*reco_p + reco_mass * reco_mass);
0299
0300       part.vecPart.SetPxPyPzE(reco_px, reco_py, reco_pz, sqrt(reco_p*reco_p + reco_mass*reco_mass));
0301
0302
0303       //cout  << "\t\t\t track  " << Form(" %4.2f %4.2f %4.2f true id %4d imc %3d mcid %3d \n",reco_px,reco_py,reco_pz,tracks_trueID[ireco],imc,part.mcID);
0304
0305       // add to `recopart` and hadronic final state sums only if there is a matching truth particle
0306       if(part.mcID > 0 && part.mcID != genEleID) {
0307     if(imc>-1) {
0308       //  cout  << "\t\t\t add  to hadfs  \n" ;
0309       recopart.push_back(part);
0310       kin->AddToHFS(part.vecPart);
0311       if( writeHFSTree ){
0312         kin->AddToHFSTree(part.vecPart, part.pid);
0313       }
0314     }
0315       }
0316
0317       // find scattered electron, by matching to truth // TODO: not realistic... is there an upstream electron finder?
0318       if(pid_ == 11 && part.mcID == genEleID) {
0319         recEleFound++;
0320     //  cout  << "\t\t\t reco electron " << Form(" %6.2f %6.2f %6.2f %6.2f  id %3d\n",reco_px,reco_py,reco_pz,sqrt(reco_p*reco_p + reco_mass*reco_mass),genEleID);
0321         kin->vecElectron.SetPxPyPzE(reco_px, reco_py, reco_pz, sqrt(reco_p*reco_p + reco_mass*reco_mass));
0322       }
0323
0324     } // end reco loop
0325
0326
0327     //add all clusters for hadronic final state
0328     int jentryt = tr.GetCurrentEntry();
0329     Long64_t localEntry = chain->LoadTree(jentryt);
0330
0331     //loop over all calorimeters with their different names
0332     for(auto calo : calos){
0333
0334       std::string dummy =  "cluster_"+calo+"_N";
0335       if ( chain->GetBranch(dummy.data())){
0336     chain->SetBranchAddress(dummy.data(), &cluster_N, &b_clusters);
0337
0338     b_clusters->GetEntry(localEntry);
0339     b_clusters->SetAutoDelete(kTRUE);
0340
0341     Float_t         cluster_E[cluster_N];
0342     Float_t         cluster_Eta[cluster_N];
0343     Float_t         cluster_Phi[cluster_N];
0344     Int_t           cluster_NTower[cluster_N];
0345     Int_t           cluster_trueID[cluster_N];
0346     dummy = "cluster_"+calo+"_E";
0347     chain->SetBranchAddress(dummy.data(), &cluster_E, &b_cluster_E);
0348     dummy = "cluster_"+calo+"_Eta";
0349     chain->SetBranchAddress(dummy.data(), &cluster_Eta, &b_cluster_Eta);
0350     dummy = "cluster_"+calo+"_Phi";
0351     chain->SetBranchAddress(dummy.data(), &cluster_Phi, &b_cluster_Phi);
0352     dummy = "cluster_"+calo+"_NTower";
0353     chain->SetBranchAddress(dummy.data(), &cluster_NTower, &b_cluster_NTower);
0354     dummy = "cluster_"+calo+"_trueID";
0355     chain->SetBranchAddress(dummy.data(), &cluster_trueID, &b_cluster_trueID);
0356     b_cluster_E->GetEntry(localEntry);
0357     b_cluster_E->SetAutoDelete(kTRUE);
0358     b_cluster_Eta->GetEntry(localEntry);
0359     b_cluster_Eta->SetAutoDelete(kTRUE);
0360     b_cluster_Phi->GetEntry(localEntry);
0361     b_cluster_Phi->SetAutoDelete(kTRUE);
0362     b_cluster_NTower->GetEntry(localEntry);
0363     b_cluster_NTower->SetAutoDelete(kTRUE);
0364     b_cluster_trueID->GetEntry(localEntry);
0365     b_cluster_trueID->SetAutoDelete(kTRUE);
0366
0367     //loop over clusters of this calorimeter
0368     for ( int iclus=0; iclus < cluster_N; iclus ++){
0369
0370       Particles part;
0371       part.mcID = cluster_trueID[iclus];
0372       part.charge = 0;
0373       TVector3 track3;
0374       track3.SetPtEtaPhi(cluster_E[iclus]/cosh(cluster_Eta[iclus]),cluster_Eta[iclus],cluster_Phi[iclus]);
0375
0376       part.vecPart.SetPxPyPzE(track3.Px(), track3.Py(), track3.Pz(), sqrt(track3.Mag2()));
0377
0378
0379       // currently use the true PID to select clusters that are not from charged pions, kaons, protons, electrons or muons for the hadronic final state
0380       auto search = mcidmap.find((int)(cluster_trueID[iclus]));
0381       if (search != mcidmap.end()) {
0382         int imc = search->second;
0383         part.pid = (int)(mcpart_PDG[imc]);
0384         if ( fabs(part.pid) != 11 &&  fabs(part.pid) != 13 && fabs(part.pid) != 211 && fabs(part.pid) != 321 && fabs(part.pid) != 2212 && fabs(part.pid) != 0  ){
0385           // only add to hadronic final state, not to the particle list (could be changed for pi0 in the future
0386           //  recopart.push_back(part);
0387           kin->AddToHFS(part.vecPart);
0388         }
0389       }
0390     }
0391
0392       }
0393
0394     }
0395
0396
0397     // skip the event if the scattered electron is not found and we need it
0398     if(recEleFound < 1) {
0399       numNoEle++;
0400       continue; // TODO: only need to skip if we are using a recon method that needs it (`if reconMethod==...`)
0401     }
0402     else if(recEleFound>1) ErrorPrint("WARNING: found more than 1 reconstructed scattered electron in an event");
0403     else numEle++;
0404
0405     // subtract electron from hadronic final state variables
0406     kin->SubtractElectronFromHFS();
0407     kinTrue->SubtractElectronFromHFS();
0408
0409     // skip the event if there are no reconstructed particles (other than the
0410     // electron), otherwise hadronic recon methods will fail
0411     if(kin->countHadrons == 0) {
0412       numNoHadrons++;
0413       continue;
0414     };
0415
0416
0417
0418     /*    printf("hadronic final state %5.2f %5.2f %5.2f %5.2f \n", kin->sigmah, kin->Pxh, kin->Pyh, kin->hadronSumVec.E() );
0419     getchar();
0420     */
0421
0422     // calculate DIS kinematics
0423     if(!(kin->CalculateDIS(reconMethod))) continue; // reconstructed
0424     if(!(kinTrue->CalculateDIS(reconMethod))) continue; // generated (truth)
0425
0426     // Get the weight for this event's Q2
0427     auto Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]);
0428
0429     //fill HFS tree for event
0430     if( writeHFSTree && kin->nHFS > 0) HFST->FillTree(Q2weightFactor);
0431
0432     // fill inclusive histograms, if only `inclusive` is included in output
0433     // (otherwise they will be filled in track and jet loops)
0434     if(includeOutputSet["inclusive_only"]) {
0435       auto wInclusive = Q2weightFactor * weightInclusive->GetWeight(*kinTrue);
0436       wInclusiveTotal += wInclusive;
0437       FillHistosInclusive(wInclusive);
0438     }
0439
0440     // loop over reconstructed particles again
0441     /* - calculate hadron kinematics
0442      * - fill output data structures (Histos, SidisTree, etc.)
0443      */
0444     for(auto part : recopart) {
0445       int pid_ = part.pid;
0446       int mcid_ = part.mcID;
0447
0448       // final state cut
0449       // - check PID, to see if it's a final state we're interested in for
0450       //   histograms; if not, proceed to next track
0451       auto kv = PIDtoFinalState.find(pid_);
0452       if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
0453       if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
0454
0455       // calculate reconstructed hadron kinematics
0456       kin->vecHadron = part.vecPart;
0457       kin->CalculateHadronKinematics();
0458
0459       // add selected single hadron FS to HFS tree
0460       if( writeHFSTree ){
0461         kin->AddTrackToHFSTree(part.vecPart, part.pid);
0462       }
0463
0464       // find the matching truth hadron using mcID, and calculate its kinematics
0465       if(mcid_ > 0) {
0466         for(auto imc : mcpart) {
0467           if(mcid_ == imc.mcID) {
0468             kinTrue->vecHadron = imc.vecPart;
0469         // add tracks of interest for kinematic studies to HFSTree
0470             if( writeHFSTree ){
0471               kinTrue->AddTrackToHFSTree(imc.vecPart, imc.pid);
0472             }
0473             break;
0474           }
0475         }
0476       }
0477       /* // deprecated, since existence of truth match is checked earlier; in practice prox matching was never called
0478      else {
0479      // give it another shot: proximity matching
0480      double mineta = 4.0;
0481      numProxMatched++;
0482      for(int imc=0; imc<(int)mcpart.size(); imc++) {
0483      if(pid_ == mcpart[imc].pid) {
0484      double deta = abs(kin->vecHadron.Eta() - mcpart[imc].vecPart.Eta());
0485             if(deta < mineta) {
0486               mineta = deta;
0487               kinTrue->vecHadron = mcpart[imc].vecPart;
0488             }
0489           }
0490         }
0491       }
0492       */
0493       kinTrue->CalculateHadronKinematics();
0494
0495       // asymmetry injection
0496       // kin->InjectFakeAsymmetry(); // sets tSpin, based on reconstructed kinematics
0497       // kinTrue->InjectFakeAsymmetry(); // sets tSpin, based on generated kinematics
0498       // kin->tSpin = kinTrue->tSpin; // copy to "reconstructed" tSpin
0499
0500       // weighting
0501       auto wTrack = Q2weightFactor * weightTrack->GetWeight(*kinTrue);
0502       wTrackTotal += wTrack;
0503
0504       if(includeOutputSet["1h"]) {
0505         // fill single-hadron histograms in activated bins
0506         FillHistos1h(wTrack);
0507         FillHistosInclusive(wTrack);
0508
0509         // fill simple tree
0510         // - not binned
0511         // - `IsActiveEvent()` is only true if at least one bin gets filled for this track
0512         if( writeSidisTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
0513       }
0514
0515     }//hadron loop
0516
0517     if( writeHFSTree && kin->nHFS > 0) HFST->FillTree(Q2weightFactor);
0518
0519   } while(tr.Next()); // tree reader loop
0520   cout << "end event loop" << endl;
0521   // event loop end =========================================================
0522
0523   // finish execution
0524   Finish();
0525
0526   // final printout
0527   cout << "Total number of scattered electrons found: " << numEle << endl;
0528   if(numNoEle>0)
0529     cerr << "WARNING: skipped " << numNoEle << " events which had no reconstructed scattered electron" << endl;
0530   if(numNoHadrons>0)
0531     cerr << "WARNING: skipped " << numNoHadrons << " events which had no reconstructed hadrons" << endl;
0532   if(numNoBeam>0)
0533     cerr << "WARNING: skipped " << numNoBeam << " events which had no beam particles" << endl;
0534   if(numProxMatched>0)
0535     cerr << "WARNING: " << numProxMatched << " recon. particles were proximity matched to truth (when mcID match failed)" << endl;
0536
0537 }