File indexing completed on 2024-09-27 07:03:26
0001
0002
0003 #ifdef INCLUDE_PODIO
0004
0005 #include "AnalysisEpicPodio.h"
0006
0007 AnalysisEpicPodio::AnalysisEpicPodio(TString infileName_, TString outfilePrefix_)
0008 : Analysis(infileName_, outfilePrefix_)
0009 , crossCheckKinematics(false)
0010 {};
0011
0012 AnalysisEpicPodio::~AnalysisEpicPodio() {};
0013
0014 void AnalysisEpicPodio::Execute()
0015 {
0016
0017 Prepare();
0018
0019
0020 std::vector<std::string> infilesFlat;
0021 for(const auto fileList : infiles)
0022 for(const auto fileName : fileList)
0023 infilesFlat.push_back(fileName);
0024
0025
0026 podioReader.openFiles(infilesFlat);
0027 evStore.setReader(&podioReader);
0028 ENT = podioReader.getEntries();
0029 if(maxEvents>0) ENT = std::min(maxEvents,ENT);
0030
0031
0032 CalculateEventQ2Weights();
0033
0034
0035
0036 const std::vector<std::string> upstreamReconMethodList = {
0037 "Truth",
0038 "Electron",
0039 "DA",
0040 "JB",
0041 "Sigma"
0042 };
0043
0044
0045 const std::map<TString,std::string> associatedUpstreamMethodMap = {
0046 { "Ele", "Electron" },
0047 { "DA", "DA" },
0048 { "JB", "JB" },
0049 { "Sigma", "Sigma" },
0050 { "Mixed", "NONE" },
0051 { "eSigma", "NONE" }
0052 };
0053
0054 const auto& associatedUpstreamMethod = associatedUpstreamMethodMap.at(reconMethod);
0055
0056
0057 fmt::print("begin event loop...\n");
0058 for(unsigned e=0; e<ENT; e++) {
0059 if(e%10000==0) fmt::print("{} events...\n",e);
0060 if(verbose) fmt::print("\n\n{:=<70}\n",fmt::format("EVENT {} ",e));
0061
0062
0063
0064 if(e>0) {
0065 evStore.clear();
0066 podioReader.endOfEvent();
0067 }
0068
0069
0070 kin->ResetHFS();
0071 kinTrue->ResetHFS();
0072 double mcPartElectronP = 0.0;
0073 bool double_counted_beam = false;
0074 int num_ele_beams = 0;
0075 int num_ion_beams = 0;
0076 int num_sim_electrons = 0;
0077 int num_rec_electrons = 0;
0078
0079
0080 const auto& simParts = evStore.get<edm4hep::MCParticleCollection>("MCParticles");
0081 const auto& recParts = evStore.get<edm4eic::ReconstructedParticleCollection>("ReconstructedParticles");
0082 const auto& mcRecAssocs = evStore.get<edm4eic::MCRecoParticleAssociationCollection>("ReconstructedParticlesAssoc");
0083
0084
0085 edm4hep::MCParticle mcPartEleBeam;
0086 edm4hep::MCParticle mcPartIonBeam;
0087 edm4hep::MCParticle mcPartElectron;
0088
0089
0090 if(verbose) fmt::print("\n{:-<60}\n","MCParticles ");
0091 for(auto simPart : simParts) {
0092
0093
0094
0095
0096
0097 auto simPDG = simPart.getPDG();
0098
0099
0100 kinTrue->AddToHFS(GetP4(simPart));
0101
0102
0103 if(simPart.getGeneratorStatus() == constants::statusBeam) {
0104 switch(simPDG) {
0105 case constants::pdgElectron:
0106 if(num_ele_beams>0) double_counted_beam = true;
0107 mcPartEleBeam = simPart;
0108 num_ele_beams++;
0109 break;
0110 case constants::pdgProton:
0111 if(num_ion_beams>0) double_counted_beam = true;
0112 mcPartIonBeam = simPart;
0113 num_ion_beams++;
0114 break;
0115 default:
0116 ErrorPrint(fmt::format("WARNING: Unknown beam particle with PDG={}",simPDG));
0117 }
0118 }
0119
0120
0121 if(simPart.getGeneratorStatus() == constants::statusFinal) {
0122 if(simPDG == constants::pdgElectron) {
0123 auto eleP = edm4hep::utils::p(simPart);
0124 if(eleP>mcPartElectronP) {
0125 mcPartElectron = simPart;
0126 mcPartElectronP = eleP;
0127 num_sim_electrons++;
0128 }
0129 }
0130 }
0131
0132 }
0133
0134
0135 if(num_ele_beams==0) { ErrorPrint("WARNING: missing MC electron beam"); continue; };
0136 if(num_ion_beams==0) { ErrorPrint("WARNING: missing MC ion beam"); continue; };
0137 if(num_sim_electrons==0) { ErrorPrint("WARNING: missing scattered electron"); continue; };
0138 if(double_counted_beam) { ErrorPrint("WARNING: found multiple beam particles"); continue; };
0139
0140
0141 kinTrue->vecEleBeam = GetP4(mcPartEleBeam);
0142 kinTrue->vecIonBeam = GetP4(mcPartIonBeam);
0143 kinTrue->vecElectron = GetP4(mcPartElectron);
0144
0145
0146 if(verbose) {
0147 if(verbose) fmt::print("\n{:-<60}\n","GENERATED BEAMS ");
0148 PrintParticle(mcPartEleBeam);
0149 PrintParticle(mcPartIonBeam);
0150 if(verbose) fmt::print("\n{:-<60}\n","GENERATED SCATTERED ELECTRON ");
0151 PrintParticle(mcPartElectron);
0152 }
0153
0154
0155
0156
0157
0158
0159
0160
0161 auto AllRecPartsToHFS = [&] (auto& simPart, auto& recPart, auto recPDG) {
0162 kin->AddToHFS(GetP4(recPart));
0163 };
0164 if(verbose) fmt::print("\n{:-<60}\n","MC<->Reco ASSOCIATIONS ");
0165 LoopMCRecoAssocs(mcRecAssocs, AllRecPartsToHFS, verbose);
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192 const auto& disCalcs = evStore.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsElectron");
0193 if(disCalcs.size() != 1) ErrorPrint(fmt::format("WARNING: disCalcs.size = {} != 1 for this event",disCalcs.size()));
0194 for(const auto& calc : disCalcs) {
0195 auto ele = calc.getScat();
0196 if( ! ele.isAvailable()) {
0197 ErrorPrint("WARNING: `disCalcs` scattered electron unavailable");
0198 continue;
0199 }
0200 num_rec_electrons++;
0201 kin->vecElectron = GetP4(ele);
0202 }
0203
0204
0205 if(num_rec_electrons == 0) { ErrorPrint("WARNING: reconstructed scattered electron not found"); continue; };
0206 if(num_rec_electrons > 1) { ErrorPrint("WARNING: found more than 1 reconstructed scattered electron"); };
0207
0208
0209 kin->SubtractElectronFromHFS();
0210 kinTrue->SubtractElectronFromHFS();
0211
0212
0213
0214 if(kin->countHadrons == 0) { ErrorPrint("WARNING: no hadrons"); };
0215
0216
0217 if( ! kin->CalculateDIS(reconMethod) ) continue;
0218 if( ! kinTrue->CalculateDIS(reconMethod) ) continue;
0219
0220
0221
0222
0223
0224
0225 auto Q2weightFactor = 1.0;
0226
0227
0228
0229 if(includeOutputSet["inclusive_only"]) {
0230 auto wInclusive = Q2weightFactor * weightInclusive->GetWeight(*kinTrue);
0231 wInclusiveTotal += wInclusive;
0232 FillHistosInclusive(wInclusive);
0233 }
0234
0235
0236
0237
0238
0239 auto SidisOutput = [&] (auto& simPart, auto& recPart, auto recPDG) {
0240
0241
0242
0243 auto kv = PIDtoFinalState.find(recPDG);
0244 if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else return;
0245 if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) return;
0246
0247
0248 kinTrue->vecHadron = GetP4(simPart);
0249 kinTrue->CalculateHadronKinematics();
0250 kin->vecHadron = GetP4(recPart);
0251 kin->CalculateHadronKinematics();
0252
0253
0254 auto wTrack = Q2weightFactor * weightTrack->GetWeight(*kinTrue);
0255 wTrackTotal += wTrack;
0256
0257
0258 FillHistos1h(wTrack);
0259 FillHistosInclusive(wTrack);
0260
0261
0262
0263
0264 if( writeSidisTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
0265 };
0266 LoopMCRecoAssocs(mcRecAssocs, SidisOutput);
0267
0268
0269
0270
0271 if(crossCheckKinematics) {
0272 auto PrintRow = [] <typename T> (std::string name, std::vector<T> vals, bool header=false) {
0273 fmt::print(" {:>16}",name);
0274 if(header) { for(auto val : vals) fmt::print(" {:>8}", val); }
0275 else { for(auto val : vals) fmt::print(" {:8.4f}", val); }
0276 fmt::print("\n");
0277 };
0278
0279 fmt::print("\n{:-<75}\n","KINEMATICS, calculated from upstream: ");
0280 PrintRow("", std::vector<std::string>({ "x", "Q2", "W", "y", "nu" }), true);
0281 for(const auto upstreamReconMethod : upstreamReconMethodList)
0282 for(const auto& calc : evStore.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematics"+upstreamReconMethod) )
0283 PrintRow( upstreamReconMethod, std::vector<float>({
0284 calc.getX(),
0285 calc.getQ2(),
0286 calc.getW(),
0287 calc.getY(),
0288 calc.getNu()
0289 }));
0290
0291 fmt::print("{:-<75}\n",fmt::format("KINEMATICS, calculated locally in EPIC-ANALYSIS, with method \"{}\": ",reconMethod));
0292 auto PrintKinematics = [&PrintRow] (std::string name, std::shared_ptr<Kinematics> K) {
0293 PrintRow( name, std::vector<Double_t>({
0294 K->x,
0295 K->Q2,
0296 K->W,
0297 K->y,
0298 K->Nu
0299 }));
0300 };
0301 PrintKinematics("Truth",kinTrue);
0302 PrintKinematics("Reconstructed",kin);
0303
0304 if(associatedUpstreamMethod != "NONE") {
0305 fmt::print("{:-<75}\n",fmt::format("DIFFERENCE: upstream({}) - local({}): ",associatedUpstreamMethod,reconMethod));
0306 for(const auto upstreamMethod : std::vector<std::string>({"Truth",associatedUpstreamMethod})) {
0307 const auto& upstreamCalcs = evStore.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematics"+upstreamMethod);
0308 for(const auto& upstreamCalc : upstreamCalcs) {
0309 auto K = upstreamMethod=="Truth" ? kinTrue : kin;
0310 auto name = upstreamMethod=="Truth" ? "Truth" : "Reconstructed";
0311 PrintRow( name, std::vector<Double_t>({
0312 upstreamCalc.getX() - K->x,
0313 upstreamCalc.getQ2() - K->Q2,
0314 upstreamCalc.getW() - K->W,
0315 upstreamCalc.getY() - K->y,
0316 upstreamCalc.getNu() - K->Nu
0317 }));
0318 }
0319 }
0320 }
0321 else fmt::print("{:-<75}\n method \"{}\" is not available upstream\n","DIFFERENCE: ",reconMethod);
0322 }
0323
0324 }
0325 fmt::print("end event loop\n");
0326
0327
0328 evStore.clear();
0329 podioReader.endOfEvent();
0330 podioReader.closeFile();
0331 Finish();
0332 }
0333
0334
0335
0336
0337 void AnalysisEpicPodio::PrintParticle(const edm4hep::MCParticle& P) {
0338 fmt::print("\n");
0339 fmt::print(" {:>20}: {}\n", "PDG", P.getPDG() );
0340 fmt::print(" {:>20}: {}\n", "Status", P.getGeneratorStatus() );
0341 fmt::print(" {:>20}: {}\n", "Energy", P.getEnergy() );
0342 fmt::print(" {:>20}: {}\n", "p=|Momentum|", edm4hep::utils::p(P) );
0343 fmt::print(" {:>20}: {}\n", "pT_lab", edm4hep::utils::pT(P) );
0344 fmt::print(" {:>20}: ({}, {}, {})\n",
0345 "3-Momentum",
0346 P.getMomentum().x,
0347 P.getMomentum().y,
0348 P.getMomentum().z
0349 );
0350 fmt::print(" {:>20}: ({}, {}, {})\n",
0351 "Vertex",
0352 P.getVertex().x,
0353 P.getVertex().y,
0354 P.getVertex().z
0355 );
0356 fmt::print(" {:>20}:\n", "Parents");
0357 for(const auto& parent : P.getParents())
0358 fmt::print(" {:>20}: {}\n", "PDG", parent.getPDG());
0359 fmt::print(" {:>20}:\n", "Daughters");
0360 for(const auto& daughter : P.getDaughters())
0361 fmt::print(" {:>20}: {}\n", "PDG", daughter.getPDG());
0362 }
0363
0364 void AnalysisEpicPodio::PrintParticle(const edm4eic::ReconstructedParticle& P) {
0365 fmt::print("\n");
0366 fmt::print(" {:>20}: ", "PDG");
0367 if(P.getParticleIDUsed().isAvailable()) fmt::print("{}\n", P.getParticleIDUsed().getPDG());
0368 else fmt::print("???\n");
0369 fmt::print(" {:>20}: {}\n", "Mass", P.getMass() );
0370 fmt::print(" {:>20}: {}\n", "Charge", P.getCharge() );
0371 fmt::print(" {:>20}: {}\n", "Energy", P.getEnergy() );
0372 fmt::print(" {:>20}: {}\n", "p=|Momentum|", edm4hep::utils::p(P) );
0373 fmt::print(" {:>20}: {}\n", "pT_lab", edm4hep::utils::pT(P) );
0374 fmt::print(" {:>20}: ({}, {}, {})\n",
0375 "3-Momentum",
0376 P.getMomentum().x,
0377 P.getMomentum().y,
0378 P.getMomentum().z
0379 );
0380 fmt::print(" {:>20}: {}\n", "# of clusters", P.clusters_size() );
0381 fmt::print(" {:>20}: {}\n", "# of tracks", P.tracks_size() );
0382 fmt::print(" {:>20}: {}\n", "# of PIDs", P.particleIDs_size() );
0383 fmt::print(" {:>20}: {}\n", "# of recParts", P.particles_size() );
0384
0385
0386
0387
0388
0389
0390 }
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401 void AnalysisEpicPodio::LoopMCRecoAssocs(
0402 const edm4eic::MCRecoParticleAssociationCollection& mcRecAssocs,
0403 std::function<void(const edm4hep::MCParticle&, const edm4eic::ReconstructedParticle&, int)> payload,
0404 bool printParticles
0405 )
0406 {
0407 for(const auto& assoc : mcRecAssocs ) {
0408
0409
0410 auto recPart = assoc.getRec();
0411 auto simPart = assoc.getSim();
0412
0413
0414
0415 if(printParticles) {
0416 fmt::print("\n {:->35}\n"," reconstructed particle:");
0417 PrintParticle(recPart);
0418 fmt::print("\n {:.>35}\n"," truth match:");
0419 if(simPart.isAvailable())
0420 PrintParticle(simPart);
0421 else
0422 fmt::print(" {:>35}\n","NO MATCH");
0423 fmt::print("\n");
0424 }
0425
0426
0427 bool usedTruthPID = false;
0428 auto recPDG = GetReconstructedPDG(simPart, recPart, usedTruthPID);
0429 if(verbose) fmt::print(" GetReconstructedPDG = {}\n",recPDG);
0430
0431
0432
0433 payload(simPart, recPart, recPDG);
0434
0435 }
0436 }
0437
0438
0439
0440
0441 int AnalysisEpicPodio::GetReconstructedPDG(
0442 const edm4hep::MCParticle& simPart,
0443 const edm4eic::ReconstructedParticle& recPart,
0444 bool& usedTruth
0445 )
0446 {
0447 int pdg = 0;
0448 usedTruth = false;
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460 if(pdg==0) {
0461 usedTruth = true;
0462 if(simPart.isAvailable())
0463 pdg = simPart.getPDG();
0464 }
0465
0466 return pdg;
0467 }
0468
0469 #endif