Warning, file /snippets/JetsAndHF/JetValidation/MakeJetValidationHists.C was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include <edm4eic/EDM4eicVersion.h>
0017 #include <TChain.h>
0018 #include <TFile.h>
0019 #include <TH1D.h>
0020 #include <TH2D.h>
0021 #include <TTree.h>
0022 #include <TTreeReader.h>
0023 #include <TTreeReaderArray.h>
0024 #include <TLegend.h>
0025 #include <TVector3.h>
0026 #include <fstream>
0027 #include <string>
0028 #include <vector>
0029
0030
0031 const std::string DefaultOutFile = "hists.files26071.py8ncdis10x100q100t1000.root";
0032
0033
0034 const std::string DefaultInFileList = "filelists/files26071.py8ncdis10x100q100t1000.list";
0035
0036
0037
0038 const std::int32_t DefaultNFiles = -1;
0039
0040
0041
0042 const std::int32_t DefaultNEvents = -1;
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052 bool branchExists(TTree* tree, const std::string& branch) {
0053 bool exists = false;
0054 if(tree->GetBranch(branch.c_str())) {
0055 exists = true;
0056 }
0057 return exists;
0058 }
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071 int MakeJetValidationHists(
0072 const std::string& results_file = DefaultOutFile,
0073 const std::string& filelist = DefaultInFileList,
0074 const std::int32_t n_files = DefaultNFiles,
0075 const std::int32_t n_events = DefaultNEvents
0076 ) {
0077
0078 const bool PRINT = true;
0079
0080 std::ifstream mylist(filelist);
0081 if (!mylist.is_open()) {
0082 std::cerr << "PANIC: Couldn't open fileslist '" << filelist << "'!" << std::endl;
0083 return 1;
0084 }
0085
0086
0087 std::string file;
0088 std::int32_t i_file = 0;
0089 std::vector<std::string> rec_files;
0090 while (std::getline(mylist, file)) {
0091 rec_files.push_back(file);
0092 ++i_file;
0093 if ((n_files > -1) && (i_file == n_files)) {
0094 break;
0095 }
0096 }
0097
0098
0099 TChain *mychain = new TChain("events");
0100 for (const auto& rec_file : rec_files) {
0101 mychain->Add(rec_file.c_str());
0102 }
0103
0104
0105 TTreeReader tree_reader(mychain);
0106
0107
0108 float DELTARCUT = 0.05;
0109
0110
0111
0112 bool useNewEDM = false;
0113 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8,9,0)
0114 {
0115 auto file = TFile::Open(rec_files.front().c_str(), "READ");
0116 auto tree = file->Get<TTree>("events");
0117 bool hasRecoArea = branchExists(tree, "ReconstructedChargedJets.area");
0118 bool hasGenArea = branchExists(tree, "GeneratedChargedJets.area");
0119 useNewEDM = hasRecoArea && hasGenArea;
0120 }
0121 #endif
0122 if (PRINT) {
0123 std::cout << "INFO: Using new EDM? " << useNewEDM << std::endl;
0124 }
0125
0126
0127
0128
0129 std::string recoAreaBranch = "ReconstructedChargedJets.energy";
0130 std::string genAreaBranch = "GeneratedChargedJets.energy";
0131 if(useNewEDM) {
0132 recoAreaBranch = "ReconstructedChargedJets.area";
0133 genAreaBranch = "GeneratedChargedJets.area";
0134 }
0135
0136
0137 std::string recoCstsBeginBranch = "ReconstructedChargedJets.particles_begin";
0138 std::string recoCstsEndBranch = "ReconstructedChargedJets.particles_end";
0139 std::string recoCstIndexBranch = "_ReconstructedChargedJets_particles.index";
0140 if (useNewEDM) {
0141 recoCstsBeginBranch = "ReconstructedChargedJets.constituents_begin";
0142 recoCstsEndBranch = "ReconstructedChargedJets.constituents_end";
0143 recoCstIndexBranch = "_ReconstructedChargedJets_constituents.index";
0144 }
0145
0146 std::string genCstsBeginBranch = "GeneratedChargedJets.particles_begin";
0147 std::string genCstsEndBranch = "GeneratedChargedJets.particles_end";
0148 std::string genCstIndexBranch = "_GeneratedChargedJets_particles.index";
0149 if (useNewEDM) {
0150 genCstsBeginBranch = "GeneratedChargedJets.constituents_begin";
0151 genCstsEndBranch = "GeneratedChargedJets.constituents_end";
0152 genCstIndexBranch = "_GeneratedChargedJets_constituents.index";
0153 }
0154
0155
0156 TTreeReaderArray<int> recoType = {tree_reader, "ReconstructedChargedJets.type"};
0157 TTreeReaderArray<float> recoNRG = {tree_reader, "ReconstructedChargedJets.energy"};
0158 TTreeReaderArray<float> recoMomX = {tree_reader, "ReconstructedChargedJets.momentum.x"};
0159 TTreeReaderArray<float> recoMomY = {tree_reader, "ReconstructedChargedJets.momentum.y"};
0160 TTreeReaderArray<float> recoMomZ = {tree_reader, "ReconstructedChargedJets.momentum.z"};
0161 TTreeReaderArray<float> recoArea = {tree_reader, recoAreaBranch.c_str()};
0162 TTreeReaderArray<unsigned int> recoCstsBegin = {tree_reader, recoCstsBeginBranch.c_str()};
0163 TTreeReaderArray<unsigned int> recoCstsEnd = {tree_reader, recoCstsEndBranch.c_str()};
0164
0165 TTreeReaderArray<int> recoCstIndex = {tree_reader, recoCstIndexBranch.c_str()};
0166
0167
0168 TTreeReaderArray<float> recoPartMomX = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
0169 TTreeReaderArray<float> recoPartMomY = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
0170 TTreeReaderArray<float> recoPartMomZ = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0171 TTreeReaderArray<float> recoPartM = {tree_reader, "ReconstructedChargedParticles.mass"};
0172 TTreeReaderArray<int> recoPartPDG = {tree_reader, "ReconstructedChargedParticles.PDG"};
0173 TTreeReaderArray<float> recoPartNRG = {tree_reader, "ReconstructedChargedParticles.energy"};
0174
0175 TTreeReaderArray<int> recoPartAssocRec = {tree_reader, "_ReconstructedChargedParticleAssociations_rec.index"};
0176 TTreeReaderArray<int> recoPartAssocSim = {tree_reader, "_ReconstructedChargedParticleAssociations_sim.index"};
0177 TTreeReaderArray<float> recoPartAssocWeight = {tree_reader, "ReconstructedChargedParticleAssociations.weight"};
0178
0179
0180 TTreeReaderArray<int> genType = {tree_reader, "GeneratedChargedJets.type"};
0181 TTreeReaderArray<float> genNRG = {tree_reader, "GeneratedChargedJets.energy"};
0182 TTreeReaderArray<float> genMomX = {tree_reader, "GeneratedChargedJets.momentum.x"};
0183 TTreeReaderArray<float> genMomY = {tree_reader, "GeneratedChargedJets.momentum.y"};
0184 TTreeReaderArray<float> genMomZ = {tree_reader, "GeneratedChargedJets.momentum.z"};
0185 TTreeReaderArray<float> genArea = {tree_reader, genAreaBranch.c_str()};
0186 TTreeReaderArray<unsigned int> genCstsBegin = {tree_reader, genCstsBeginBranch.c_str()};
0187 TTreeReaderArray<unsigned int> genCstsEnd = {tree_reader, genCstsEndBranch.c_str()};
0188
0189 TTreeReaderArray<int> genPartIndex = {tree_reader, genCstIndexBranch.c_str()};
0190
0191
0192
0193
0194 TTreeReaderArray<float> mcMomX = {tree_reader, "GeneratedParticles.momentum.x"};
0195 TTreeReaderArray<float> mcMomY = {tree_reader, "GeneratedParticles.momentum.y"};
0196 TTreeReaderArray<float> mcMomZ = {tree_reader, "GeneratedParticles.momentum.z"};
0197 TTreeReaderArray<float> mcM = {tree_reader, "GeneratedParticles.mass"};
0198 TTreeReaderArray<int> pdg = {tree_reader, "GeneratedParticles.PDG"};
0199
0200 TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
0201 TTreeReaderArray<double> mcMomXPart = {tree_reader, "MCParticles.momentum.x"};
0202 TTreeReaderArray<double> mcMomYPart = {tree_reader, "MCParticles.momentum.y"};
0203 TTreeReaderArray<double> mcMomZPart = {tree_reader, "MCParticles.momentum.z"};
0204 TTreeReaderArray<double> mcMPart = {tree_reader, "MCParticles.mass"};
0205 TTreeReaderArray<int> pdgMCPart = {tree_reader, "MCParticles.PDG"};
0206
0207
0208 TH1D *counter = new TH1D("counter","",10,0.,10.);
0209
0210
0211 TFile* out_file = new TFile(results_file.c_str(), "recreate");
0212 std::cout << "INFO: writing histograms to " << out_file->GetName() << std::endl;
0213
0214
0215 TH1D *numRecoChargedJetsECutHist = new TH1D("numRecoChargedJetsECut","",20,0.,20.);
0216 TH1D *recoChargedJetEHist = new TH1D("recoChargedJetE","",300,0.,300.);
0217 TH1D *recoChargedJetEtaECutHist = new TH1D("recoChargedJetEtaECut","",60,-3.,3.);
0218 TH1D *recoChargedJetAreaECutHist = nullptr;
0219 TH2D *recoChargedJetEvsAreaHist = nullptr;
0220 if (useNewEDM) {
0221 recoChargedJetAreaECutHist = new TH1D("recoChargedJetAreaECut","",250,0.,5.);
0222 recoChargedJetEvsAreaHist = new TH2D("recoChargedJetEvsArea","",250,0.,5.,300,0.,300.);
0223 }
0224 TH2D *recoChargedJetEvsEtaHist = new TH2D("recoChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
0225 TH2D *recoChargedJetPhiVsEtaECutHist = new TH2D("recoChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0226
0227 TH1D *numRecoChargedJetsECutNoElecHist = new TH1D("numRecoChargedJetsECutNoElec","",20,0.,20.);
0228 TH1D *recoChargedJetENoElecHist = new TH1D("recoChargedJetENoElec","",300,0.,300.);
0229 TH1D *recoChargedJetEtaECutNoElecHist = new TH1D("recoChargedJetEtaECutNoElec","",60,-3.,3.);
0230 TH1D *recoChargedJetAreaECutNoElecHist = nullptr;
0231 TH2D *recoChargedJetEvsAreaNoElecHist = nullptr;
0232 if (useNewEDM) {
0233 recoChargedJetAreaECutNoElecHist = new TH1D("recoHargedJetAreaECutNoElec","",250,0.,5.);
0234 recoChargedJetEvsAreaNoElecHist = new TH2D("recoChargedJetEvsAreaNoElec","",250,0.,5.,300,0.,300.);
0235 }
0236 TH2D *recoChargedJetEvsEtaNoElecHist = new TH2D("recoChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
0237 TH2D *recoChargedJetPhiVsEtaECutNoElecHist = new TH2D("recoChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0238
0239 TH1D *numRecoChargedJetPartsHist = new TH1D("numRecoChargedJetParts","",20,0.,20.);
0240 TH1D *recoChargedJetPartPHist = new TH1D("recoChargedJetPartP","",500,0.,100.);
0241 TH1D *recoChargedJetPartEtaHist = new TH1D("recoChargedJetPartEta","",80,-4.,4.);
0242 TH2D *recoChargedJetPartPvsEtaHist = new TH2D("recoChargedJetPartPvsEta","",80,-4.,4.,500,0.,100.);
0243 TH2D *recoChargedJetPartPhiVsEtaHist = new TH2D("recoChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0244
0245 TH1D *numRecoChargedJetPartsNoElecHist = new TH1D("numRecoChargedJetPartsNoElec","",20,0.,20.);
0246 TH1D *recoChargedJetPartPNoElecHist = new TH1D("recoChargedJetPartPNoElec","",500,0.,100.);
0247 TH1D *recoChargedJetPartEtaNoElecHist = new TH1D("recoChargedJetPartEtaNoElec","",80,-4.,4.);
0248 TH2D *recoChargedJetPartPvsEtaNoElecHist = new TH2D("recoChargedJetPartPvsEtaNoElec","",80,-4.,4.,500,0.,100.);
0249 TH2D *recoChargedJetPartPhiVsEtaNoElecHist = new TH2D("recoChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0250
0251 TH1D *recoChargedJetPartPairwiseDeltaRHist = new TH1D("recoChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
0252
0253
0254 TH1D *numGenChargedJetsECutHist = new TH1D("numGenChargedJetsECut","",20,0.,20.);
0255 TH1D *genChargedJetEHist = new TH1D("genChargedJetE","",300,0.,300.);
0256 TH1D *genChargedJetEtaECutHist = new TH1D("genChargedJetEtaECut","",60,-3.,3.);
0257 TH1D *genChargedJetAreaECutHist = nullptr;
0258 TH2D *genChargedJetEvsAreaHist = nullptr;
0259 if (useNewEDM) {
0260 genChargedJetAreaECutHist = new TH1D("genChargedJetAreaECut","",250,0.,5.);
0261 genChargedJetEvsAreaHist = new TH2D("genChargedJetEvsAreaHist","",250,0.,5.,300,0.,300.);
0262 }
0263 TH2D *genChargedJetEvsEtaHist = new TH2D("genChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
0264 TH2D *genChargedJetPhiVsEtaECutHist = new TH2D("genChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0265
0266 TH1D *numGenChargedJetsECutNoElecHist = new TH1D("numGenChargedJetsECutNoElec","",20,0.,20.);
0267 TH1D *genChargedJetENoElecHist = new TH1D("genChargedJetENoElec","",300,0.,300.);
0268 TH1D *genChargedJetEtaECutNoElecHist = new TH1D("genChargedJetEtaECutNoElec","",60,-3.,3.);
0269 TH1D *genChargedJetAreaECutNoElecHist = nullptr;
0270 TH2D *genChargedJetEvsAreaNoElecHist = nullptr;
0271 if (useNewEDM) {
0272 genChargedJetAreaECutNoElecHist = new TH1D("genChargedJetAreaECutNoElec","",250,0.,5.);
0273 genChargedJetEvsAreaNoElecHist = new TH2D("genChargedJetEvsAreaNoElec","",250,0.,5.,300,0.,300.);
0274 }
0275 TH2D *genChargedJetEvsEtaNoElecHist = new TH2D("genChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
0276 TH2D *genChargedJetPhiVsEtaECutNoElecHist = new TH2D("genChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0277
0278 TH1D *numGenChargedJetPartsHist = new TH1D("numGenChargedJetParts","",20,0.,20.);
0279 TH1D *genChargedJetPartPHist = new TH1D("genChargedJetPartP","",500,0.,100.);
0280 TH1D *genChargedJetPartEtaHist = new TH1D("genChargedJetPartEta","",80,-4.,4.);
0281 TH2D *genChargedJetPartPvsEtaHist = new TH2D("genChargedJetPartPvsEta","",80,-4.,4.,500,0.,100.);
0282 TH2D *genChargedJetPartPhiVsEtaHist = new TH2D("genChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0283
0284 TH1D *numGenChargedJetPartsNoElecHist = new TH1D("numGenChargedJetPartsNoElec","",20,0.,20.);
0285 TH1D *genChargedJetPartPNoElecHist = new TH1D("genChargedJetPartPNoElec","",500,0.,100.);
0286 TH1D *genChargedJetPartEtaNoElecHist = new TH1D("genChargedJetPartEtaNoElec","",80,-4.,4.);
0287 TH2D *genChargedJetPartPvsEtaNoElecHist = new TH2D("genChargedJetPartPvsEtaNoElec","",80,-4.,4.,500,0.,100.);
0288 TH2D *genChargedJetPartPhiVsEtaNoElecHist = new TH2D("genChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0289
0290 TH1D *genChargedJetPartPairwiseDeltaRHist = new TH1D("genChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
0291
0292
0293 TH1D *matchJetDeltaRHist = new TH1D("matchJetDeltaR","",5000,0.,5.);
0294 TH1D *matchJetDeltaRBackHist = new TH1D("matchJetDeltaRBack","",5000,0.,5.);
0295 TH2D *recoVsGenChargedJetEtaHist = new TH2D("recoVsGenChargedJetEta","",80,-4.,4.,80,-4.,4.);
0296 TH2D *recoVsGenChargedJetPhiHist = new TH2D("recoVsGenChargedJetPhi","",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
0297 TH2D *recoVsGenChargedJetAreaHist = nullptr;
0298 if (useNewEDM) {
0299 recoVsGenChargedJetAreaHist = new TH2D("recoVsGenChargedJetArea","",250,0.,5.,250,0.,5.);
0300 }
0301 TH2D *recoVsGenChargedJetEHist = new TH2D("recoVsGenChargedJetE","",100,0.,100.,100,0.,100.);
0302 TH2D *recoVsGenChargedJetENoDRHist = new TH2D("recoVsGenChargedJetENoDRHist","",100,0.,100.,100,0.,100.);
0303 TH2D *recoVsGenChargedJetENoDupHist = new TH2D("recoVsGenChargedJetENoDup","",100,0.,100.,100,0.,100.);
0304
0305 TH2D *jetResVsEtaHist = new TH2D("jetResVsEta","",80,-4.,4.,10000,-10.,10.);
0306 TH2D *jetResVsEHist = new TH2D("jetResVsE","",100,0.,100.,10000,-10.,10.);
0307 TH2D *jetResVsENegEtaHist = new TH2D("jetResVsENegEta","",20,0.,100.,10000,-10.,10.);
0308 TH2D *jetResVsEMidEtaHist = new TH2D("jetResVsEMidEta","",20,0.,100.,10000,-10.,10.);
0309 TH2D *jetResVsEPosEtaHist = new TH2D("jetResVsEPosEta","",20,0.,100.,10000,-10.,10.);
0310
0311 TH2D *jetResVsENegEtaNoDupHist = new TH2D("jetResVsENegEtaNoDup","",20,0.,100.,10000,-10.,10.);
0312 TH2D *jetResVsEMidEtaNoDupHist = new TH2D("jetResVsEMidEtaNoDup","",20,0.,100.,10000,-10.,10.);
0313 TH2D *jetResVsEPosEtaNoDupHist = new TH2D("jetResVsEPosEtaNoDup","",20,0.,100.,10000,-10.,10.);
0314
0315
0316
0317 int NEVENTS = 0;
0318 while(tree_reader.Next()) {
0319
0320 if ((n_events > -1) && (NEVENTS == n_events)) {
0321 break;
0322 }
0323 if (NEVENTS%10000 == 0) cout << "Events Processed: " << NEVENTS << endl;
0324
0325 counter->Fill(0);
0326
0327
0328
0329
0330 int numRecoChargedJets = 0;
0331 int numRecoChargedJetsNoElec = 0;
0332 for(unsigned int i=0; i<recoType.GetSize(); i++)
0333 {
0334 TVector3 jetMom(recoMomX[i],recoMomY[i],recoMomZ[i]);
0335
0336 counter->Fill(3);
0337
0338
0339 if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0340
0341
0342 bool ECut = recoNRG[i] > 5.0;
0343
0344 if(ECut) numRecoChargedJets++;
0345
0346 recoChargedJetEHist->Fill(recoNRG[i]);
0347 if(ECut) recoChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
0348 recoChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0349 if(useNewEDM) {
0350 if(ECut) recoChargedJetAreaECutHist->Fill(recoArea[i]);
0351 recoChargedJetEvsAreaHist->Fill(recoArea[i],recoNRG[i]);
0352 }
0353 if(ECut) recoChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0354
0355
0356 bool noElectron = true;
0357 for(unsigned int m=recoCstsBegin[i]; m<recoCstsEnd[i]; m++)
0358 {
0359 int elecIndex = -1;
0360 double elecIndexWeight = -1.0;
0361 int chargePartIndex = recoCstIndex[m];
0362 for(unsigned int n=0; n<recoPartAssocRec.GetSize(); n++)
0363 {
0364 if(recoPartAssocRec[n] == chargePartIndex)
0365 {
0366 if(recoPartAssocWeight[n] > elecIndexWeight)
0367 {
0368 elecIndex = recoPartAssocSim[n];
0369 elecIndexWeight = recoPartAssocWeight[n];
0370 }
0371 }
0372 }
0373
0374 if(pdgMCPart[elecIndex] == 11)
0375 noElectron = false;
0376 }
0377
0378 if(ECut)
0379 {
0380 for(unsigned int j=recoCstsBegin[i]; j<recoCstsEnd[i]; j++)
0381 {
0382
0383
0384 double mX = recoPartMomX[recoCstIndex[j]];
0385 double mY = recoPartMomY[recoCstIndex[j]];
0386 double mZ = recoPartMomZ[recoCstIndex[j]];
0387 double mM = recoPartM[recoCstIndex[j]];
0388
0389
0390 TVector3 partMom(mX,mY,mZ);
0391
0392 recoChargedJetPartPHist->Fill(partMom.Mag());
0393 recoChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
0394 recoChargedJetPartPvsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0395 recoChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0396
0397 if(noElectron)
0398 {
0399 recoChargedJetPartPNoElecHist->Fill(partMom.Mag());
0400 recoChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
0401 recoChargedJetPartPvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0402 recoChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0403 }
0404
0405
0406 if(j<(recoCstsEnd[i]-1))
0407 {
0408 for(unsigned int k=j+1; k<recoCstsEnd[i]; k++)
0409 {
0410 double mXB = recoPartMomX[recoCstIndex[k]];
0411 double mYB = recoPartMomY[recoCstIndex[k]];
0412 double mZB = recoPartMomZ[recoCstIndex[k]];
0413
0414 TVector3 partMomB(mXB,mYB,mZB);
0415
0416 double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0417 double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0418 double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0419
0420 recoChargedJetPartPairwiseDeltaRHist->Fill(dR);
0421 }
0422 }
0423 }
0424 numRecoChargedJetPartsHist->Fill(recoCstsEnd[i] - recoCstsBegin[i]);
0425 if(noElectron) numRecoChargedJetPartsNoElecHist->Fill(recoCstsEnd[i] - recoCstsBegin[i]);
0426 }
0427
0428
0429 if(noElectron)
0430 {
0431 recoChargedJetENoElecHist->Fill(recoNRG[i]);
0432 if(ECut) recoChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
0433 recoChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0434 if(useNewEDM) {
0435 if(ECut) recoChargedJetAreaECutNoElecHist->Fill(recoArea[i]);
0436 recoChargedJetEvsAreaNoElecHist->Fill(recoArea[i],recoNRG[i]);
0437 }
0438 if(ECut) recoChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0439
0440 if(ECut) numRecoChargedJetsNoElec++;
0441 }
0442 }
0443 numRecoChargedJetsECutHist->Fill(numRecoChargedJets);
0444 numRecoChargedJetsECutNoElecHist->Fill(numRecoChargedJetsNoElec);
0445
0446
0447
0448
0449 int numGenChargedJets = 0;
0450 int numGenChargedJetsNoElec = 0;
0451 for(unsigned int i=0; i<genType.GetSize(); i++)
0452 {
0453 TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0454
0455 counter->Fill(4);
0456
0457
0458 if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0459
0460
0461 bool ECut = genNRG[i] > 5.0;
0462
0463 if(ECut) numGenChargedJets++;
0464
0465 genChargedJetEHist->Fill(genNRG[i]);
0466 if(ECut) genChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
0467 genChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0468 if(useNewEDM) {
0469 if(ECut) genChargedJetAreaECutHist->Fill(genArea[i]);
0470 genChargedJetEvsAreaHist->Fill(genArea[i],genNRG[i]);
0471 }
0472 if(ECut) genChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0473
0474
0475 bool noElectron = true;
0476 for(unsigned int m=genCstsBegin[i]; m<genCstsEnd[i]; m++)
0477 {
0478 if(pdg[genPartIndex[m]] == 11)
0479 noElectron = false;
0480 }
0481
0482 if(ECut)
0483 {
0484 for(unsigned int j=genCstsBegin[i]; j<genCstsEnd[i]; j++)
0485 {
0486
0487
0488 double mX = mcMomX[genPartIndex[j]];
0489 double mY = mcMomY[genPartIndex[j]];
0490 double mZ = mcMomZ[genPartIndex[j]];
0491 double mM = mcM[genPartIndex[j]];
0492
0493
0494 TVector3 partMom(mX,mY,mZ);
0495
0496 genChargedJetPartPHist->Fill(partMom.Mag());
0497 genChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
0498 genChargedJetPartPvsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0499 genChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0500
0501 if(noElectron)
0502 {
0503 genChargedJetPartPNoElecHist->Fill(partMom.Mag());
0504 genChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
0505 genChargedJetPartPvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0506 genChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0507 }
0508
0509
0510 if(j<(genCstsEnd[i]-1))
0511 {
0512 for(unsigned int k=j+1; k<genCstsEnd[i]; k++)
0513 {
0514 double mXB = mcMomX[genPartIndex[k]];
0515 double mYB = mcMomY[genPartIndex[k]];
0516 double mZB = mcMomZ[genPartIndex[k]];
0517
0518 TVector3 partMomB(mXB,mYB,mZB);
0519
0520 double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0521 double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0522 double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0523
0524 genChargedJetPartPairwiseDeltaRHist->Fill(dR);
0525 }
0526 }
0527 }
0528 numGenChargedJetPartsHist->Fill(genCstsEnd[i] - genCstsBegin[i]);
0529 if(noElectron) numGenChargedJetPartsNoElecHist->Fill(genCstsEnd[i] - genCstsBegin[i]);
0530 }
0531
0532
0533 if(noElectron)
0534 {
0535 genChargedJetENoElecHist->Fill(genNRG[i]);
0536 if(ECut) genChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
0537 genChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0538 if (useNewEDM) {
0539 if(ECut) genChargedJetAreaECutNoElecHist->Fill(genArea[i]);
0540 genChargedJetEvsAreaHist->Fill(genArea[i],genNRG[i]);
0541 }
0542 if(ECut) genChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0543
0544 if(ECut) numGenChargedJetsNoElec++;
0545 }
0546 }
0547 numGenChargedJetsECutHist->Fill(numGenChargedJets);
0548 numGenChargedJetsECutNoElecHist->Fill(numGenChargedJetsNoElec);
0549
0550
0551
0552
0553
0554 for(unsigned int i=0; i<genType.GetSize(); i++)
0555 {
0556 TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0557
0558
0559
0560
0561
0562
0563
0564
0565 bool hasElectron = false;
0566
0567 for(unsigned int m=genCstsBegin[i]; m<genCstsEnd[i]; m++)
0568 {
0569 if(pdg[genPartIndex[m]] == 11)
0570 hasElectron = true;
0571 }
0572
0573
0574
0575 double minDeltaR = 999.;
0576 int minIndex = -1;
0577 for(unsigned int j=0; j<recoType.GetSize(); j++)
0578 {
0579 TVector3 recoMom(recoMomX[j],recoMomY[j],recoMomZ[j]);
0580
0581 double dEta = jetMom.PseudoRapidity() - recoMom.PseudoRapidity();
0582 double dPhi = TVector2::Phi_mpi_pi(jetMom.Phi() - recoMom.Phi());
0583 double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0584
0585 if(dR < minDeltaR)
0586 {
0587 minDeltaR = dR;
0588 minIndex = j;
0589 }
0590 }
0591
0592
0593 double minDeltaRBack = 999.;
0594 double minIndexBack = -1;
0595 if(minIndex > -1)
0596 {
0597 TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
0598 for(unsigned int j=0; j<genType.GetSize(); j++)
0599 {
0600 TVector3 genMom(genMomX[j],genMomY[j],genMomZ[j]);
0601
0602 double dEta = recoMatchMom.PseudoRapidity() - genMom.PseudoRapidity();
0603 double dPhi = TVector2::Phi_mpi_pi(recoMatchMom.Phi() - genMom.Phi());
0604 double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0605
0606 if(dR < minDeltaRBack)
0607 {
0608 minDeltaRBack = dR;
0609 minIndexBack = j;
0610 }
0611 }
0612 }
0613
0614
0615 if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1 && !hasElectron) matchJetDeltaRHist->Fill(minDeltaR);
0616 if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1) matchJetDeltaRBackHist->Fill(minDeltaR);
0617 if(minIndex > -1 && genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && !hasElectron)
0618 {
0619 TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
0620
0621 recoVsGenChargedJetENoDRHist->Fill(genNRG[i],recoNRG[minIndex]);
0622
0623 if(minDeltaR < DELTARCUT)
0624 {
0625 recoVsGenChargedJetEtaHist->Fill(jetMom.PseudoRapidity(),recoMatchMom.PseudoRapidity());
0626 recoVsGenChargedJetPhiHist->Fill(jetMom.Phi(),recoMatchMom.Phi());
0627 if(useNewEDM) {
0628 recoVsGenChargedJetAreaHist->Fill(genArea[i],recoArea[minIndex]);
0629 }
0630 recoVsGenChargedJetEHist->Fill(genNRG[i],recoNRG[minIndex]);
0631
0632 double jetERes = (recoNRG[minIndex] - genNRG[i])/genNRG[i];
0633
0634 jetResVsEtaHist->Fill(jetMom.PseudoRapidity(),jetERes);
0635 jetResVsEHist->Fill(genNRG[i],jetERes);
0636 if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
0637 jetResVsENegEtaHist->Fill(genNRG[i],jetERes);
0638 if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
0639 jetResVsEMidEtaHist->Fill(genNRG[i],jetERes);
0640 if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
0641 jetResVsEPosEtaHist->Fill(genNRG[i],jetERes);
0642
0643
0644 bool noDuplicate = true;
0645 for(unsigned int j=recoCstsBegin[minIndex]; j<recoCstsEnd[minIndex]; j++)
0646 {
0647 double mX = recoPartMomX[recoCstIndex[j]];
0648 double mY = recoPartMomY[recoCstIndex[j]];
0649 double mZ = recoPartMomZ[recoCstIndex[j]];
0650 double mM = recoPartM[recoCstIndex[j]];
0651 double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0652
0653 TVector3 partMom(mX,mY,mZ);
0654
0655
0656 if(j<(recoCstsEnd[minIndex]-1))
0657 {
0658 for(unsigned int k=j+1; k<recoCstsEnd[minIndex]; k++)
0659 {
0660 double mXB = recoPartMomX[recoCstIndex[k]];
0661 double mYB = recoPartMomY[recoCstIndex[k]];
0662 double mZB = recoPartMomZ[recoCstIndex[k]];
0663
0664 TVector3 partMomB(mXB,mYB,mZB);
0665
0666 double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0667 double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0668 double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0669
0670 if(dR < 0.02) noDuplicate = false;
0671 }
0672 }
0673 }
0674
0675 if(noDuplicate)
0676 {
0677 recoVsGenChargedJetENoDupHist->Fill(genNRG[i],recoNRG[minIndex]);
0678
0679 if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
0680 jetResVsENegEtaNoDupHist->Fill(genNRG[i],jetERes);
0681 if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
0682 jetResVsEMidEtaNoDupHist->Fill(genNRG[i],jetERes);
0683 if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
0684 jetResVsEPosEtaNoDupHist->Fill(genNRG[i],jetERes);
0685 }
0686 }
0687 }
0688 }
0689
0690 NEVENTS++;
0691 }
0692
0693
0694 out_file->cd();
0695 out_file->Write();
0696 out_file->Close();
0697
0698 delete mychain;
0699 return 0;
0700 }
