irt/scripts/evaluation.C

0001 R__LOAD_LIBRARY(podioDict)
0002 R__LOAD_LIBRARY(podioRootIO)
0003 R__LOAD_LIBRARY(libedm4hepDict)
0004 #include "podio/EventStore.h"
0005 #include "podio/ROOTReader.h"
0006 #include "edm4hep/utils/kinematics.h"
0007
0008 // #define _AEROGEL_
0009
0010 #define _NPE_REFERENCE_ 211
0011 //#define _NPE_REFERENCE_ (11)
0012 //#define _NPE_REFERENCE_ 321
0013
0014 void evaluation(const char *ifname, const char *ofname = 0)
0015 {
0016   // open recon output .root file with podio
0017   podio::EventStore podio_store;
0018   podio::ROOTReader podio_reader;
0019   podio_reader.openFile(ifname);
0020   podio_store.setReader(&podio_reader);
0021
0022   auto np = new TH1D("np", "Photon count",            50,       0,       50);
0023 #ifdef _AEROGEL_
0024   unsigned id = 0;
0025   auto th = new TH1D("th", "Cherenkov theta",         50,     180,      200);
0026   auto ri = new TH1D("ri", "Refractive Index - 1.0",  50,   0.018,    0.020);
0027   auto dt = new TH1D("dt", "Cherenkov theta diff",    50,     -10,       10);
0028   //auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -1,        1);
0029 #else
0030   unsigned id = 1;
0031   auto th = new TH1D("th", "",         50,      35,       41);
0032   auto tq = new TH1D("tq", "",         50,      35,       41);
0033   auto ri = new TH1D("ri", "Refractive Index - 1.0",  50, 0.00075,  0.00077);
0034   auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -3,        3);
0035   //auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -5,        5);
0036 #endif
0037   auto wl = new TH1D("wl", "Wave length",             50,     350,      900);
0038
0039   unsigned false_assignment_stat[2] = {0};
0040   std::vector<double> thvector, npvector;
0041
0042   // event loop
0043   for(unsigned ev=0; ev<podio_reader.getEntries(); ev++) {
0044     if(ev%100==0) printf("read event %d\n",ev);
0045
0046     // get collections
0047     auto& cherenkovs = podio_store.get<edm4eic::CherenkovParticleIDCollection>("DRICHPID");
0048     auto& mctracks   = podio_store.get<edm4hep::MCParticleCollection>("MCParticles");
0049
0050     // Then the Cherenkov-to-reconstructed mapping;
0051     // FIXME: may want to use Cherenkov-to-simulated mapping to start with, for the debugging purposes;
0052     // FIXME: if we loop over reconstructed tracks, rather than MC particles, then we
0053     // have the 1-1 relation edm4hep::ReconstructedParticle::getParticleIDUsed(),
0054     // which returns the type edm4hep::ParticleID
0055     std::map<edm4hep::MCParticle,edm4eic::CherenkovParticleID> rc2cherenkov;
0056     for(const auto &pid : cherenkovs)
0057       rc2cherenkov[pid.getAssociatedParticle()] = pid;
0058
0059     // Loop through all MC tracks;
0060     for(const auto &mctrack : mctracks) {
0061       // FIXME: consider only primaries for now?; equivalent to mctrack.getGeneratorStatus()==1?
0062       if (mctrack.parents_size()>0) continue;
0063
0064       edm4eic::CherenkovParticleID cherenkov;
0065       if(rc2cherenkov.find(mctrack) != rc2cherenkov.end()) cherenkov = rc2cherenkov[mctrack];
0066       else continue;
0067
0068       double pp = edm4hep::utils::p(mctrack);
0069       double m  = mctrack.getMass();
0070
0071       //printf("m=%5.3f p=%5.1f (%4d) \n", m, pp, mctrack.getPDG());
0072
0073       // Loop through all of the mass hypotheses available for this reconstructed track;
0074       {
0075         const edm4eic::CherenkovPdgHypothesis *best = 0;
0076
0077         for(const auto &option : cherenkov.getOptions()) {
0078
0079           if (option.radiator != id) continue;
0080
0081           // Skip electron hypothesis; of no interest here;
0082           //if (abs(option.pdg) == 11) continue;
0083
0084           if (abs(option.pdg) == _NPE_REFERENCE_) {
0085             np->Fill(option.npe);
0086
0087             if (ofname) npvector.push_back(option.npe);
0088           } //if
0089
0090           if (!best || option.weight > best->weight) best = &option;
0091           printf("radiator %3d (pdg %5d): weight %7.2f, npe %7.2f\n",
0092               option.radiator, option.pdg, option.weight, option.npe);
0093         } //for option
0094         printf("\n");
0095
0096         // Check whether the true PDG got a highest score;
0097         if (!best || best->pdg != mctrack.getPDG()) false_assignment_stat[best->npe >= 5 ? 0 : 1]++;
0098
0099         // This assumes of course that at least one radiator was requested in juggler;
0100         double rindex = cherenkov.getAngles()[id].rindex;
0101         double theta  = cherenkov.getAngles()[id].theta;
0102         double lambda = cherenkov.getAngles()[id].wavelength;
0103         double argument = sqrt(pp*pp + m*m)/(rindex*pp);
0104         double thp = fabs(argument) <= 1.0 ? acos(argument) : theta;
0105
0106         th->Fill(1000 * theta);
0107         //if (mctrack.getPDG() == 321)
0108         dt->Fill(1000* (theta - thp));
0109         ri->Fill(rindex - 1.0);
0110         wl->Fill(lambda);//rindex - 1.0);
0111         printf("<n> ~ %8.6f, <th> = %7.2f [mrad]\n", rindex - 1.0, 1000*thp);
0112
0113         if (ofname) thvector.push_back(theta - thp);
0114       }
0115     } //for track
0116
0117     // next event
0118     podio_store.clear();
0119     podio_reader.endOfEvent();
0120   } //for ev
0121
0122   // end of event loop
0123   printf("%3d (%3d) false out of %d\n", false_assignment_stat[0], false_assignment_stat[1], podio_reader.getEntries());
0124   podio_reader.closeFile();
0125
0126   // write
0127   if (ofname) {
0128
0129     auto *ofdata = new TFile(ofname, "RECREATE");
0130
0131     if (!ofdata) {
0132       printf("was not able to create output file '%s'\n", ofname);
0133       exit(0);
0134     } //if
0135     auto *ot = new TTree("t", "My tree");
0136
0137     double thbff, npbff;
0138     ot->Branch("th", &thbff, "th/D");
0139     ot->Branch("np", &npbff, "np/D");
0140
0141     for(unsigned iq=0; iq<thvector.size(); iq++) {
0142       thbff = thvector[iq];
0143       npbff = npvector[iq];
0144
0145       ot->Fill();
0146     } //for iq
0147
0148     ot->Write();
0149     ofdata->Close();
0150     exit(0);
0151   } else {
0152     auto cv = new TCanvas("cv", "", 1700, 500);
0153     cv->Divide(5, 1);
0154     cv->cd(1); np->Draw();       np->Fit("gaus");
0155     cv->cd(2); th->Draw("SAME"); th->Fit("gaus");
0156     cv->cd(3); ri->Draw();       ri->Fit("gaus");
0157     cv->cd(4); dt->Draw();       dt->Fit("gaus");
0158     cv->cd(5); wl->Draw();       //dt->Fit("gaus");
0159   } //if
0160 } // evaluation()