irt/scripts/evaluation.fig12.C

0001
0002 #define _DETECTOR_ "DRICH"
0003 //#define _DETECTOR_ "PFRICH"
0004 //#define _AEROGEL_
0005 //#define _NPE_REFERENCE_ 211
0006 #define _NPE_REFERENCE_ 321
0007 //#define _USE_RECONSTRUCTED_TRACKS_
0008
0009 #define _TSIZE_ 0.07
0010 #define _LSIZE_ 0.05
0011
0012 void evaluation(const char *ifname, const char *ofname = 0)
0013 {
0014   const unsigned gdim = 11, qdim = 40;
0015   double thmin = 36.0, thstep = 0.1;
0016   unsigned thstat[2][qdim]; memset(thstat, 0x00, sizeof(thstat));
0017
0018   // .root file with event tree;
0019   auto ifdata = new TFile(ifname);
0020   if (!ifdata) {
0021     printf("input file '%s' does not exist\n", ifname);
0022     exit(0);
0023   } //if
0024   TTree *it = dynamic_cast<TTree*>(ifdata->Get("events"));
0025   if (!it) {
0026     printf("input file '%s' does not have \"events\" tree\n", ifname);
0027     exit(0);
0028   } //if
0029
0030   std::vector<double> thvector, npvector;
0031
0032   auto np = new TH1D("np", "Photon count",            50,       0,       50);
0033 #ifdef _AEROGEL_
0034   unsigned id = 0;
0035   auto th = new TH1D("th", "Cherenkov theta",         50,     180,      200);
0036   auto ri = new TH1D("ri", "Refractive Index - 1.0",  50,   0.018,    0.020);
0037   auto dt = new TH1D("dt", "Cherenkov theta diff",    50,     -10,       10);
0038   //auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -1,        1);
0039 #else
0040   unsigned id = 1;
0041   auto th = new TH1D("th", "",         50,      35,       41);
0042   auto tq = new TH1D("tq", "",         50,      35,       41);
0043   auto ri = new TH1D("ri", "Refractive Index - 1.0",  50, 0.00075,  0.00077);
0044   auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -2,        3);
0045   //auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -5,        5);
0046 #endif
0047
0048   // Use MC truth particles for a "main" loop;
0049   auto mctracks   = new std::vector<dd4pod::Geant4ParticleData>();
0050   auto rctracks   = new std::vector<edm4eic::ReconstructedParticleData>();
0051   auto cherenkov  = new std::vector<edm4eic::CherenkovParticleIDData>();
0052   it->SetBranchAddress("mcparticles", &mctracks);
0053
0054   // FIXME: or whatever the branches are called;
0055 #ifdef _USE_RECONSTRUCTED_TRACKS_
0056   it->SetBranchAddress("rcparticles", &rctracks);
0057 #endif
0058   it->SetBranchAddress((TString(_DETECTOR_) + "PID").Data(),   &cherenkov);
0059   auto options = new std::vector<edm4eic::CherenkovPdgHypothesis>();
0060   it->SetBranchAddress((TString(_DETECTOR_) + "PID_0").Data(), &options);
0061   auto angles  = new std::vector<edm4eic::CherenkovThetaAngleMeasurement>();
0062   it->SetBranchAddress((TString(_DETECTOR_) + "PID_1").Data(), &angles);
0063
0064   // Loop through all events;
0065   unsigned false_assignment_stat[2] = {0};
0066   for(int ev=0; ev<it->GetEntries(); ev++) {
0067     it->GetEntry(ev);
0068
0069 #ifdef _USE_RECONSTRUCTED_TRACKS_
0070     // First populate the reconstructed-to-simulated particle mapping table;
0071     std::map<eic::Index, const edm4eic::ReconstructedParticleData*> mc2rc;
0072     for(const auto &rctrack: *rctracks)
0073       mc2rc[rctrack.mcID] = &rctrack;
0074 #endif
0075
0076     // Then the Cherenkov-to-reconstructed mapping; FIXME: may want to use Cherenkov-to-simulated
0077     // mapping to start with, for the debugging purposes;
0078     std::map<eic::Index, const edm4eic::CherenkovParticleIDData*> rc2cherenkov;
0079     for(const auto &pid: *cherenkov)
0080       rc2cherenkov[pid.recID] = &pid;
0081
0082     //printf("Here! %d\n", mctracks->size());
0083     // Loop through all MC tracks;
0084     for(auto mctrack: *mctracks) {
0085       // FIXME: consider only primaries for now?;
0086       if (mctrack.g4Parent) continue;
0087
0088 #ifdef _USE_RECONSTRUCTED_TRACKS_
0089       // Find a matching reconstructed track;
0090       auto rctrack = mc2rc.find(mctrack.ID) == mc2rc.end() ? 0 : mc2rc[mctrack.ID];
0091       if (!rctrack) continue;
0092
0093       // Find a matching Cherenkov PID record;
0094       auto cherenkov = rc2cherenkov.find(rctrack.ID) == rc2cherenkov.end() ? 0 : rc2cherenkov[rctrack.ID];
0095 #else
0096       auto cherenkov = rc2cherenkov.find(mctrack.ID) == rc2cherenkov.end() ? 0 : rc2cherenkov[mctrack.ID];
0097 #endif
0098       if (!cherenkov) continue;
0099
0100
0101
0102       // BACK !!!;
0103       double pp = mctrack.ps.mag(), m = 0.493677;//mctrack.mass;
0104
0105
0106
0107
0108       printf("%f %f (%4d) \n", mctrack.mass, mctrack.ps.mag(), mctrack.pdgID);
0109
0110       // Loop through all of the mass hypotheses available for this reconstructed track;
0111       {
0112     const edm4eic::CherenkovPdgHypothesis *best = 0;
0113
0114     for(unsigned iq=cherenkov->options_begin; iq<cherenkov->options_end; iq++) {
0115       const auto &option = (*options)[iq];
0116
0117       if (option.radiator != id) continue;
0118
0119       // Skip electron hypothesis; of no interest here;
0120       if (abs(option.pdg) == 11) continue;
0121
0122       if (abs(option.pdg) == _NPE_REFERENCE_) {
0123         np->Fill(option.npe);
0124
0125         if (ofname) npvector.push_back(option.npe);
0126       } //if
0127
0128       if (!best || option.weight > best->weight) best = &option;
0129       printf("radiator %3d (pdg %5d): weight %7.2f, npe %7.2f\n",
0130          option.radiator, option.pdg, option.weight, option.npe);
0131     } //for
0132     printf("\n");
0133
0134     //if (!best) printf("@J@\n");
0135     // Check whether the true PDG got a highest score;
0136     if (!best || best->pdg != mctrack.pdgID) {
0137       //printf("@J@ %7.2f\n", best->npe);
0138       false_assignment_stat[best->npe >= 0 ? 0 : 1]++;
0139     }
0140     //}
0141
0142     // This assumes of course that at least one radiator was requested in juggler;
0143     //{
0144     double rindex = (*angles)[id].rindex, theta = (*angles)[id].theta;
0145     double argument = sqrt(pp*pp + m*m)/(rindex*pp);
0146     double thp = fabs(argument) <= 1.0 ? acos(argument) : theta;
0147
0148     //th->Fill(1000*  theta);
0149     if (best->npe > 15) (mctrack.pdgID == 321 ? th : tq)->Fill(1000*  theta);
0150     //if (best->npe > 20)
0151     if (mctrack.pdgID == 321) dt->Fill(1000* (theta - thp));
0152     ri->Fill(rindex - 1.0);
0153     printf("<n> ~ %8.6f, <th> = %7.2f [mrad]\n", rindex - 1.0, 1000*thp);
0154
0155     {
0156       int thbin = (int)floor((1000*theta - thmin)/thstep);
0157       thstat[mctrack.pdgID == 321 ? 0 : 1][thbin]++;
0158     }
0159
0160     if (ofname) thvector.push_back(theta - thp);
0161       }
0162     } //for track
0163   } //for ev
0164
0165   printf("%3d (%3d) false out of %lld\n", false_assignment_stat[0],
0166      false_assignment_stat[1], it->GetEntries());
0167
0168   if (ofname) {
0169     ifdata->Close();
0170
0171     auto *ofdata = new TFile(ofname, "RECREATE");
0172
0173     if (!ofdata) {
0174       printf("was not able to create output file '%s'\n", ofname);
0175       exit(0);
0176     } //if
0177     auto *ot = new TTree("t", "My tree");
0178
0179     double thbff, npbff;
0180     ot->Branch("th", &thbff, "th/D");
0181     ot->Branch("np", &npbff, "np/D");
0182
0183     for(unsigned iq=0; iq<thvector.size(); iq++) {
0184       //for(auto thvar: thvector) {
0185       thbff = thvector[iq];
0186       npbff = npvector[iq];
0187
0188       ot->Fill();
0189     } //for iq
0190
0191     ot->Write();
0192     ofdata->Close();
0193     exit(0);
0194   } else {
0195     gROOT->Reset();
0196
0197     gStyle->SetTextSize(0.02);
0198     gStyle->SetLabelSize(0.04,"xy");
0199     gStyle->SetFrameFillColor(0);
0200     gStyle->SetOptStat(0);
0201     gStyle->SetOptFit(0);
0202     //gStyle->SetTitleSize(0.05,"x");
0203     gStyle->SetPadBottomMargin(0.15);
0204     gStyle->SetPadTopMargin(0.04);
0205     gStyle->SetPadRightMargin(0.05);
0206     gStyle->SetPadLeftMargin(0.12);
0207     //+  gStyle->SetMarkerColor(kBlack);
0208     //+  gStyle->SetMarkerStyle(25);
0209     //+  gStyle->SetMarkerSize(1.0);
0210
0211     gStyle->SetStatBorderSize(0);
0212     gStyle->SetStatColor(kWhite);
0213     gStyle->SetStatFontSize(0.03);
0214     gStyle->SetStatFont(52);
0215     gStyle->SetStatW(.13);
0216     gStyle->SetFitFormat("5.2f");
0217
0218     //gStyle->SetTitleFontSize(0.05);
0219     //gStyle->SetTitleFillColor(kWhite);
0220
0221     TCanvas *cv = new TCanvas("cv", "", 0, 0, 800, 1000);
0222     cv->UseCurrentStyle();
0223     cv->SetBorderMode(0);
0224     cv->SetFrameBorderMode(0);
0225     cv->SetFrameLineColor(kWhite);
0226     cv->SetFillColor(kWhite);
0227
0228     //auto cv = new TCanvas("cv", "", 800, 1000);
0229     cv->Divide(1, 2);
0230     th->SetLineWidth(2); tq->SetLineWidth(2);
0231
0232     tq->GetXaxis()->SetLabelFont(52);
0233     tq->GetYaxis()->SetLabelFont(52);
0234     tq->GetXaxis()->SetTitleFont(52);
0235     tq->GetYaxis()->SetTitleFont(52);
0236
0237     tq->GetXaxis()->SetTitle("Reconstructed Cherenkov angle, [mrad]");
0238     tq->GetXaxis()->SetTitleSize(_TSIZE_);
0239     tq->GetXaxis()->SetLabelSize(_LSIZE_);
0240     tq->GetYaxis()->SetTitle("Number of Events");
0241     tq->GetYaxis()->SetTitleSize(_TSIZE_);
0242     tq->GetYaxis()->SetLabelSize(_LSIZE_);
0243     tq->GetXaxis()->SetTitleOffset(0.90);
0244     tq->GetYaxis()->SetTitleOffset(0.90);
0245
0246     tq->GetXaxis()->SetNdivisions(408);
0247     tq->GetYaxis()->SetNdivisions(808);
0248
0249     cv->cd(1);
0250     gPad->SetGrid();
0251     tq->Draw();       tq->Fit("gaus");
0252     th->Draw("SAME"); th->Fit("gaus");
0253
0254     {
0255       //for(unsigned iq=0; iq<qdim; iq++)
0256       //    printf("%3d %3d\n", thstat[0][iq], thstat[1][iq]);
0257
0258       // First calculate total kaon sample stat;
0259       unsigned kstat = 0, pistat = 0, offset = 16;
0260       for(unsigned iq=0; iq<qdim; iq++) {
0261     kstat  += thstat[0][iq];
0262     pistat += thstat[1][iq];
0263       } //for iq
0264       // This is not dramatically efficient;
0265       unsigned kaccu [qdim]; memset( kaccu, 0x00, sizeof( kaccu));
0266       unsigned piaccu[qdim]; memset(piaccu, 0x00, sizeof(piaccu));
0267       for(unsigned iq=0; iq<qdim; iq++) {
0268     kaccu [iq] = (iq ? kaccu [iq-1] : 0) + thstat[0][iq];
0269     piaccu[iq] = (iq ? piaccu[iq-1] : 0) + thstat[1][iq];
0270       } //for iq
0271       for(unsigned iq=0; iq<qdim; iq++)
0272     printf("%2d: %4d -> %5.3f   %4d -> %5.3f\n",
0273            iq, kaccu[iq], 1.*kaccu[iq]/kstat, piaccu[iq], 1.*piaccu[iq]/pistat);
0274
0275       //float   eff[gdim] = { 0.70, 0.72, 0.75, 0.78, 0.80, 0.83, 0.85, 0.87, 0.90, 0.95, 0.99};
0276       //float   ctm[gdim] = { 0.00, 0.01, 0.03, 0.05, 0.07, 0.09, 0.12, 0.15, 0.22, 0.30, 0.33};
0277       float   eff[gdim];
0278       float   ctm[gdim];
0279       for(unsigned ig=0; ig<gdim; ig++) {
0280     unsigned ibin = ig + offset;
0281     eff[ig] = 1.* kaccu[ibin]/kstat;
0282     ctm[ig] = 1.*piaccu[ibin]/(kaccu[ibin] + piaccu[ibin]);
0283       } //for ig
0284       float  eeff[gdim] = { .001, .001, .001, .001, .001, .001, .001, .001, .001, .001, .001};
0285       float  ectm[gdim] = { .001, .001, .001, .001, .001, .001, .001, .001, .001, .001, .001};
0286
0287       auto greff = new TGraphErrors(gdim, eff, ctm, eeff, ectm);
0288       greff->SetMarkerSize(1.3);
0289       greff->SetMarkerStyle(21);
0290       greff->SetMarkerColor(kBlue);
0291
0292       auto hdum = new TH1D("hdum", "", 10, 0.50, 1.00);
0293
0294       hdum->GetXaxis()->SetLabelFont(52);
0295       hdum->GetYaxis()->SetLabelFont(52);
0296       hdum->GetXaxis()->SetTitleFont(52);
0297       hdum->GetYaxis()->SetTitleFont(52);
0298
0299       hdum->GetXaxis()->SetTitle("Kaon detection efficiency, [0..1]");
0300       hdum->GetXaxis()->SetTitleSize(_TSIZE_);
0301       hdum->GetXaxis()->SetLabelSize(_LSIZE_);
0302       hdum->GetYaxis()->SetTitle("Pion contamination, [0..1]");
0303       hdum->GetYaxis()->SetTitleSize(_TSIZE_);
0304       hdum->GetYaxis()->SetLabelSize(_LSIZE_);
0305       hdum->GetXaxis()->SetTitleOffset(0.90);
0306       hdum->GetYaxis()->SetTitleOffset(0.90);
0307       hdum->SetMinimum(    0.00);
0308       hdum->SetMaximum(    0.25);
0309
0310       hdum->GetXaxis()->SetNdivisions(408);
0311       hdum->GetYaxis()->SetNdivisions(808);
0312
0313       cv->cd(2);
0314       gPad->SetGrid();
0315       hdum->Draw();
0316       greff->Draw("PC");
0317       //grctm->Draw("P");
0318     }
0319
0320 #if 0
0321     auto cv = new TCanvas("cv", "", 1500, 500);
0322     cv->Divide(4, 1);
0323     cv->cd(1); np->Draw();       np->Fit("gaus");
0324     cv->cd(2); tq->Draw();       tq->Fit("gaus");
0325                th->Draw("SAME"); th->Fit("gaus");
0326     cv->cd(3); ri->Draw();       ri->Fit("gaus");
0327     cv->cd(4); dt->Draw();       //dt->Fit("gaus");
0328 #endif
0329   } //if
0330 } // evaluation()