delphes_EIC/SimpleAnalysis/TaggingStudyModule.cc

0001 #include "TaggingStudyModule.h"
0002
0003 #include "TClonesArray.h"
0004
0005 #include "AnalysisFunctions.cc"
0006
0007 #include <iostream>
0008 #include <iomanip>
0009 #include <fstream>
0010
0011 #include "TreeHandler.h"
0012
0013 TaggingStudyModule::TaggingStudyModule(ExRootTreeReader* data)
0014   : Module(data)
0015 {
0016
0017 }
0018
0019 TaggingStudyModule::~TaggingStudyModule()
0020 {
0021
0022 }
0023
0024 void TaggingStudyModule::initialize()
0025 {
0026   TreeHandler *tree_handler = tree_handler->getInstance();
0027
0028   if (tree_handler->getTree() != nullptr) {
0029     tree_handler->getTree()->Branch("jet_pt", &_jet_pt, "jet_pt/F");
0030     tree_handler->getTree()->Branch("jet_eta", &_jet_eta, "jet_eta/F");
0031     tree_handler->getTree()->Branch("jet_flavor", &_jet_flavor, "jet_flavor/F");
0032     // tree_handler->getTree()->Branch("jet_tagged", &_jet_tagged, "jet_tagged/F");
0033     tree_handler->getTree()->Branch("jet_btag", &_jet_btag, "jet_btag/F");
0034   }
0035
0036   // create the study variations map
0037   study_variations = std::map<std::string, std::map<std::string, int>>();
0038
0039   study_variations["all"]["charm"] = 0;
0040   study_variations["all"]["light"] = 0;
0041   for (auto minTrack : minTrackVar) {
0042     for (auto minTrkPt : minTrkPTVar) {
0043       for (auto minSignif : minSignifVar) {
0044     std::string varName = std::string(Form("MinTrk: %d;TrkPT: %.2f;MinSig: %.2f",minTrack,minTrkPt,minSignif));
0045     //if (study_variations.find(varName) == study_variations.end()
0046     study_variations[varName]["charm"] = 0;
0047     study_variations[varName]["light"] = 0;
0048       }
0049     }
0050   }
0051 }
0052
0053 void TaggingStudyModule::finalize()
0054 {
0055   std::vector<std::string> jets{"light", "charm"};
0056   // Report the yield for each variation as well as the Punzi Figure of Merit
0057
0058   float allLight = float(study_variations["all"]["light"]);
0059   float allCharm = float(study_variations["all"]["charm"]);
0060
0061   float maxFOM = 0.0;
0062   std::string bestVariation = "";
0063
0064   std::cout << std::setw(40) << "VARIATION" << std::setw(6) << "LIGHT"
0065         << std::setw(6) << "CHARM" << std::setw(8) << "PUNZI" << std::endl;
0066
0067
0068   ofstream csvfile;
0069   csvfile.open("tagging_study.csv");
0070
0071   csvfile << "Variation,Light,Charm" << std::endl;
0072
0073   for (auto& [variation, jet] : study_variations) {
0074     float nLight = float(study_variations[variation]["light"]);
0075     float nCharm = float(study_variations[variation]["charm"]);
0076     float PunziFOM = -1;
0077     if (allCharm>0)
0078       PunziFOM = (nCharm/allCharm)/( 1.5 + TMath::Sqrt(nLight));
0079
0080     std::cout << std::setw(40) << variation << std::setw(6) << int(nLight)
0081           << std::setw(6) << int(nCharm) << std::setw(8) << Form("%.4f", PunziFOM) << std::endl;
0082
0083     csvfile << "\"" << variation << "\"," << int(nLight) << "," << nCharm << std::endl;
0084
0085     if (PunziFOM > maxFOM) {
0086       maxFOM = PunziFOM;
0087       bestVariation = variation;
0088     }
0089   }
0090
0091   std::cout << "======================================================================" << std::endl;
0092   std::cout << "Best Variation: " << bestVariation << " (FOM: " << maxFOM << ")" << std::endl;
0093
0094   csvfile.close();
0095
0096
0097 }
0098 bool TaggingStudyModule::execute(std::map<std::string, std::any>* DataStore)
0099 {
0100   auto data = getData();
0101   TreeHandler *tree_handler = tree_handler->getInstance();
0102
0103   // If event contains at least 1 jet
0104   std::vector<Jet*> all_jets;
0105   for (int ijet = 0; ijet < getJets()->GetEntries(); ijet++)
0106     {
0107       // Take first jet
0108       Jet *jet = (Jet*) getJets()->At(ijet);
0109       all_jets.push_back(jet);
0110     }
0111
0112   // MET cut first
0113   bool passed = true;
0114
0115
0116   // Get the MET object and use it
0117   MissingET* MET = nullptr;
0118   for (int imet = 0; imet < getMET()->GetEntries(); imet++) {
0119     MET = static_cast<MissingET*>(getMET()->At(imet));
0120   }
0121
0122   if (MET == nullptr) {
0123     passed = false;
0124   }
0125
0126   if (passed == true && MET->MET > 10.0) {
0127     passed = true;
0128   } else {
0129     passed = false;
0130   }
0131
0132
0133   if (passed == false)
0134     return false;
0135
0136   std::vector<Jet*> fiducial_jets = SelectorFcn<Jet>(all_jets, [](Jet* j){ return (TMath::Abs(j->Eta) < 3.0 && j->PT > 5.0); });
0137
0138   std::vector<Jet*> charmJets = SelectorFcn<Jet>(fiducial_jets, [](Jet* j){ return (j->Flavor == 4); });
0139
0140   std::vector<Jet*> lightJets = SelectorFcn<Jet>(fiducial_jets, [](Jet* j){ return (j->Flavor < 4 || j->Flavor == 21); });
0141
0142   // Resolution settings
0143   // float d0err = 0.020; // mm
0144   // float z0err = 0.020; // mm
0145
0146
0147   // Clone the tracks and adjust their errors to suit this study
0148   auto ModifiedTracks = static_cast<TClonesArray*>(getEFlowTracks()->Clone());
0149   for (int i = 0; i < ModifiedTracks->GetEntries(); i++) {
0150     auto track = static_cast<Track*>(ModifiedTracks->At(i));
0151     // track->ErrorD0 = d0err;
0152     // track->ErrorDZ = z0err;
0153   }
0154
0155
0156   study_variations["all"]["charm"] += charmJets.size();
0157   for (auto charmjet : charmJets) {
0158     _jet_pt = charmjet->PT;
0159     _jet_eta = charmjet->Eta;
0160     _jet_flavor = charmjet->Flavor;
0161     // _jet_tagged = Tagged_sIP3D(charmjet);
0162     _jet_btag = charmjet->BTag;
0163     // tree_handler->getTree()->Fill();
0164
0165     for (auto minTrack : minTrackVar) {
0166       for (auto minTrkPt : minTrkPTVar) {
0167     for (auto minSignif : minSignifVar) {
0168       std::string varName = std::string(Form("MinTrk: %d;TrkPT: %.2f;MinSig: %.2f",minTrack,minTrkPt,minSignif));
0169       if (Tagged_sIP3D(charmjet, *ModifiedTracks, minSignif, minTrkPt, minTrack))
0170         study_variations[varName]["charm"] += 1;
0171     }
0172       }
0173     }
0174
0175   }
0176
0177   study_variations["all"]["light"] += lightJets.size();
0178   for (auto lightjet : lightJets) {
0179     _jet_pt = lightjet->PT;
0180     _jet_eta = lightjet->Eta;
0181     _jet_flavor = lightjet->Flavor;
0182     // _jet_tagged = Tagged_sIP3D(lightjet);
0183     _jet_btag = lightjet->BTag;
0184     // tree_handler->getTree()->Fill();
0185
0186     for (auto minTrack : minTrackVar) {
0187       for (auto minTrkPt : minTrkPTVar) {
0188     for (auto minSignif : minSignifVar) {
0189       std::string varName = std::string(Form("MinTrk: %d;TrkPT: %.2f;MinSig: %.2f",minTrack,minTrkPt,minSignif));
0190       if (Tagged_sIP3D(lightjet, *ModifiedTracks, minSignif, minTrkPt, minTrack))
0191         study_variations[varName]["light"] += 1;
0192     }
0193       }
0194     }
0195   }
0196
0197   if (ModifiedTracks)
0198     delete ModifiedTracks;
0199
0200   return true;
0201 }
0202
0203
0204 // bool TaggingStudyModule::Tagged(Jet* jet,
0205 //              float minSignif, float minPT, int minTracks, float errd0, float errz0)
0206 // {
0207 //   bool tagged = false;
0208
0209 //   const TLorentzVector &jetMomentum = jet->P4();
0210 //   float jpx = jetMomentum.Px();
0211 //   float jpy = jetMomentum.Py();
0212 //   float jpz = jetMomentum.Pz();
0213
0214 //   auto jet_constituents = *(getEFlowTracks());
0215
0216 //   int N_sIPtrack = 0;
0217
0218 //   for (int iconst = 0; iconst < jet_constituents.GetEntries(); iconst++) {
0219
0220
0221 //     if (N_sIPtrack >= minTracks)
0222 //       break;
0223
0224 //     auto constituent = jet_constituents.At(iconst);
0225
0226 //     if(constituent == 0) continue;
0227
0228 //     if (constituent->IsA() == Track::Class()) {
0229 //       auto track = static_cast<Track*>(constituent);
0230
0231 //       const TLorentzVector &trkMomentum = track->P4();
0232 //       float tpt = trkMomentum.Pt();
0233 //       if(tpt < minPT) continue;
0234
0235 //       if (trkMomentum.DeltaR(jetMomentum) > 0.5)
0236 //  continue;
0237
0238 //       float d0 = TMath::Abs(track->D0);
0239 //       if (d0 > 3.0)
0240 //  continue;
0241
0242 //       float xd = track->Xd;
0243 //       float yd = track->Yd;
0244 //       float zd = track->Zd;
0245 //       float dd0 = errd0;
0246 //       if (dd0 < 0)
0247 //  TMath::Abs(track->ErrorD0);
0248 //       float dz = TMath::Abs(track->DZ);
0249 //       float ddz = errz0;
0250 //       if (ddz < 0)
0251 //  TMath::Abs(track->ErrorDZ);
0252
0253 //       int sign = (jpx * xd + jpy * yd + jpz * zd > 0.0) ? 1 : -1;
0254 //       //add transverse and longitudinal significances in quadrature
0255 //       float sip = sign * TMath::Sqrt(TMath::Power(d0 / dd0, 2) + TMath::Power(dz / ddz, 2));
0256
0257 //       if(sip > minSignif) N_sIPtrack++;
0258 //     }
0259
0260
0261
0262
0263 //   }
0264
0265 //   tagged = (N_sIPtrack >= minTracks);
0266
0267 //   return tagged;
0268 // }