File indexing completed on 2024-09-27 07:03:25
0001
0002
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
0015 Prepare();
0016
0017
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
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
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
0050 TTreeReaderArray<Int_t> recparts_type(tr, "ReconstructedChargedParticles.type");
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
0059 std::string assoc_branch_name = "ReconstructedChargedParticleAssociations";
0060 if(listOfBranches->FindObject(assoc_branch_name.c_str()) == nullptr)
0061 assoc_branch_name = "ReconstructedChargedParticlesAssociations";
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
0067 CalculateEventQ2Weights();
0068
0069
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]++;
0080
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
0091 std::map <double,int> genidmap;
0092 std::map <int,int> mcidmap;
0093 std::map <int,int> recidmap;
0094 std::map <int,int> trackstatmap;
0095
0096
0097
0098 std::vector<Particles> genpart;
0099 std::vector<Particles> mcpart;
0100 std::vector<Particles> recpart;
0101
0102
0103
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
0119 Particles part;
0120
0121 part.pid=pid_;
0122
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
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
0145 Particles part;
0146
0147 part.pid=pid_;
0148
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;
0156
0157 mcidmap.insert({imc,igen});
0158
0159 }
0160
0161
0162
0163
0164
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
0178 Particles part;
0179
0180 part.pid=pid_;
0181
0182 part.vecPart.SetPxPyPzE(px_,py_,pz_,e_);
0183
0184 double m_ = part.vecPart.M();
0185
0186
0187
0188
0189
0190
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){
0199 part.mcID=idx_simID;
0200 break;
0201 }
0202 }
0203
0204 recpart.push_back(part);
0205 recidmap.insert({irec,part.mcID});
0206
0207 }
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221 for(auto mcpart_: mcpart){
0222
0223 int imc = mcpart_.mcID;
0224
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
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
0248
0249
0250 else if(genStat_ == 1 && mcidmap[mcpart_.mcID]>-1){
0251 kinTrue->AddToHFS(mcpart_.vecPart);
0252 }
0253
0254 }
0255
0256
0257 if(!foundBeamElectron || !foundBeamIon) { numNoBeam++; continue;};
0258
0259
0260
0261
0262
0263 int irec = 0;
0264 bool recEleFound=false;
0265 for(auto recpart_ : recpart){
0266
0267 if(recidmap[irec]!=-1){
0268
0269 if(recidmap[irec]==genEleID){
0270 recEleFound=true;
0271 kin->vecElectron= recpart_.vecPart;
0272 }
0273
0274 kin->AddToHFS(recpart_.vecPart);
0275
0276
0277 if( writeHFSTree ){
0278 kin->AddToHFSTree(recpart_.vecPart, pid);
0279 }
0280 }
0281 irec++;
0282 }
0283
0284
0285 if(recEleFound==false){
0286 numNoEle++;
0287 continue;
0288 }
0289 else
0290 numEle++;
0291
0292
0293 kin->SubtractElectronFromHFS();
0294 kinTrue->SubtractElectronFromHFS();
0295
0296
0297
0298 if(kin->countHadrons == 0) {
0299 numNoHadrons++;
0300 continue;
0301 };
0302
0303
0304 if(!(kin->CalculateDIS(reconMethod))) continue;
0305 if(!(kinTrue->CalculateDIS(reconMethod))) continue;
0306
0307
0308 auto Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]);
0309
0310
0311
0312 if(includeOutputSet["inclusive_only"]) {
0313 auto wInclusive = Q2weightFactor * weightInclusive->GetWeight(*kinTrue);
0314 wInclusiveTotal += wInclusive;
0315 FillHistosInclusive(wInclusive);
0316 }
0317
0318
0319
0320
0321
0322
0323
0324 int ipart = 0;
0325 for(auto part : recpart){
0326
0327 int pid_ = part.pid;
0328 int mcid_ = part.mcID;
0329
0330
0331
0332
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
0338 kin->vecHadron = part.vecPart;
0339 kin->CalculateHadronKinematics();
0340
0341
0342 if( writeHFSTree ){
0343 kin->AddTrackToHFSTree(part.vecPart, part.pid);
0344 }
0345
0346
0347
0348 if(mcid_ >= 0) {
0349 for(auto imc : mcpart) {
0350 if(mcid_ == imc.mcID) {
0351 kinTrue->vecHadron = imc.vecPart;
0352
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
0365 auto wTrack = Q2weightFactor * weightTrack->GetWeight(*kinTrue);
0366 wTrackTotal += wTrack;
0367 FillHistos1h(wTrack);
0368 FillHistosInclusive(wTrack);
0369
0370
0371
0372
0373 if( writeSidisTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
0374 }
0375 }
0376
0377 if( writeHFSTree && kin->nHFS > 0) HFST->FillTree(Q2weightFactor);
0378
0379
0380
0381
0382
0383
0384
0385 if ( writeParticleTree ) {
0386 int ipart = 0;
0387 for( auto part: recpart ){
0388 Particles mcpart_;
0389 int mcpart_idx=recidmap[ipart];
0390 int genStat_ = -1;
0391 if(mcpart_idx>-1){
0392 mcpart_ = mcpart.at(mcpart_idx);
0393 int imc = mcpart_.mcID;
0394 genStat_ = mcpart_genStat[imc];
0395 if(imc==genEleID)
0396 genStat_=2;
0397 }
0398 PT->FillTree(part.vecPart,
0399 mcpart_.vecPart,
0400 part.pid,
0401 genStat_);
0402 }
0403 ipart++;
0404 }
0405
0406 } while(tr.Next());
0407
0408
0409
0410 Finish();
0411
0412
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 }