File indexing completed on 2024-09-27 07:03:25
0001
0002
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
0013 AnalysisAthena::AnalysisAthena(TString configFileName_, TString outfilePrefix_) :
0014 Analysis(configFileName_, outfilePrefix_)
0015 { };
0016
0017
0018 AnalysisAthena::~AnalysisAthena() { };
0019
0020
0021
0022
0023
0024 void AnalysisAthena::Execute()
0025 {
0026
0027 Prepare();
0028
0029
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
0034 chain->Add(infiles[idx][idxF].c_str());
0035 }
0036 }
0037 chain->CanDeleteRefs();
0038
0039 TTreeReader tr(chain.get());
0040
0041
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
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
0065 CalculateEventQ2Weights();
0066
0067
0068 Long64_t numNoBeam, numEle, numNoEle, numNoHadrons, numProxMatched;
0069 numNoBeam = numEle = numNoEle = numNoHadrons = numProxMatched = 0;
0070
0071
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
0078 kin->ResetHFS();
0079 kinTrue->ResetHFS();
0080
0081
0082
0083
0084
0085
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];
0096 double px_ = mcparticles_psx[imc];
0097 double py_ = mcparticles_psy[imc];
0098 double pz_ = mcparticles_psz[imc];
0099 double mass_ = mcparticles_mass[imc];
0100 double p_ = sqrt(pow(mcparticles_psx[imc],2) + pow(mcparticles_psy[imc],2) + pow(mcparticles_psz[imc],2));
0101
0102 if(genStatus_ == 1) {
0103
0104
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
0112 kinTrue->AddToHFS(part.vecPart);
0113
0114
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) {
0125 if(pid_ == 11) {
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 {
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 }
0141
0142
0143 if(!foundBeamElectron || !foundBeamIon) { numNoBeam++; continue; };
0144
0145
0146
0147
0148
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;
0157
0158
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
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
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 }
0192
0193
0194 if(recEleFound < 1) {
0195 numNoEle++;
0196 continue;
0197 }
0198 else if(recEleFound>1) ErrorPrint("WARNING: found more than 1 reconstructed scattered electron in an event");
0199 else numEle++;
0200
0201
0202 kin->SubtractElectronFromHFS();
0203 kinTrue->SubtractElectronFromHFS();
0204
0205
0206
0207 if(kin->countHadrons == 0) {
0208 numNoHadrons++;
0209 continue;
0210 };
0211
0212
0213 if(!(kin->CalculateDIS(reconMethod))) continue;
0214 if(!(kinTrue->CalculateDIS(reconMethod))) continue;
0215
0216
0217 auto Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]);
0218
0219
0220
0221 if(includeOutputSet["inclusive_only"]) {
0222 auto wInclusive = Q2weightFactor * weightInclusive->GetWeight(*kinTrue);
0223 wInclusiveTotal += wInclusive;
0224 FillHistosInclusive(wInclusive);
0225 }
0226
0227
0228
0229
0230
0231 for(auto part : recopart) {
0232 int pid_ = part.pid;
0233 int mcid_ = part.mcID;
0234
0235
0236
0237
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
0243 kin->vecHadron = part.vecPart;
0244 kin->CalculateHadronKinematics();
0245
0246
0247 if( writeHFSTree ){
0248 kin->AddTrackToHFSTree(part.vecPart, part.pid);
0249 }
0250
0251
0252 if(mcid_ > 0) {
0253 for(auto imc : mcpart) {
0254 if(mcid_ == imc.mcID) {
0255 kinTrue->vecHadron = imc.vecPart;
0256
0257 if( writeHFSTree ){
0258 kinTrue->AddTrackToHFSTree(imc.vecPart, imc.pid);
0259 }
0260 break;
0261 }
0262 }
0263 }
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280 kinTrue->CalculateHadronKinematics();
0281
0282
0283
0284
0285
0286
0287
0288 auto wTrack = Q2weightFactor * weightTrack->GetWeight(*kinTrue);
0289 wTrackTotal += wTrack;
0290
0291 if(includeOutputSet["1h"]) {
0292
0293 FillHistos1h(wTrack);
0294 FillHistosInclusive(wTrack);
0295
0296
0297
0298
0299 if( writeSidisTree && HD->IsActiveEvent() ) ST->FillTree(wTrack);
0300 }
0301
0302 }
0303
0304
0305 if( writeHFSTree && kin->nHFS > 0) HFST->FillTree(Q2weightFactor);
0306
0307 } while(tr.Next());
0308 cout << "end event loop" << endl;
0309
0310
0311
0312 Finish();
0313
0314
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 }