irt/scripts/rejection.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.05
0010 #define _LSIZE_ 0.04
0011
0012 void rejection(const char *ifname)//, const char *ofname = 0)
0013 {
0014   const unsigned gdim = 11, qdim = 40, offset = 16;
0015   double thmin = 184.2, thstep = 0.2;
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", "",                        50,     180,      200);
0036   auto tq = new TH1D("tq", "",                        50,     180,      200);
0037   auto ri = new TH1D("ri", "Refractive Index - 1.0",  50,   0.018,    0.020);
0038   auto dt = new TH1D("dt", "Cherenkov theta diff",    50,     -10,       10);
0039   //auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -1,        1);
0040 #else
0041   unsigned id = 1;
0042   auto th = new TH1D("th", "",         50,      35,       41);
0043   auto tq = new TH1D("tq", "",         50,      35,       41);
0044   auto ri = new TH1D("ri", "Refractive Index - 1.0",  50, 0.00075,  0.00077);
0045   auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -2,        3);
0046   //auto dt = new TH1D("dt", "Cherenkov theta diff",    50,      -5,        5);
0047 #endif
0048
0049   // Use MC truth particles for a "main" loop;
0050   auto mctracks   = new std::vector<dd4pod::Geant4ParticleData>();
0051   auto rctracks   = new std::vector<edm4eic::ReconstructedParticleData>();
0052   auto cherenkov  = new std::vector<edm4eic::CherenkovParticleIDData>();
0053   it->SetBranchAddress("mcparticles", &mctracks);
0054
0055   // FIXME: or whatever the branches are called;
0056 #ifdef _USE_RECONSTRUCTED_TRACKS_
0057   it->SetBranchAddress("rcparticles", &rctracks);
0058 #endif
0059   it->SetBranchAddress((TString(_DETECTOR_) + "PID").Data(),   &cherenkov);
0060   auto options = new std::vector<edm4eic::CherenkovPdgHypothesis>();
0061   it->SetBranchAddress((TString(_DETECTOR_) + "PID_0").Data(), &options);
0062   auto angles  = new std::vector<edm4eic::CherenkovThetaAngleMeasurement>();
0063   it->SetBranchAddress((TString(_DETECTOR_) + "PID_1").Data(), &angles);
0064
0065   // Loop through all events;
0066   unsigned false_assignment_stat[2] = {0};
0067   for(int ev=0; ev<it->GetEntries(); ev++) {
0068     it->GetEntry(ev);
0069
0070 #ifdef _USE_RECONSTRUCTED_TRACKS_
0071     // First populate the reconstructed-to-simulated particle mapping table;
0072     std::map<eic::Index, const edm4eic::ReconstructedParticleData*> mc2rc;
0073     for(const auto &rctrack: *rctracks)
0074       mc2rc[rctrack.mcID] = &rctrack;
0075 #endif
0076
0077     // Then the Cherenkov-to-reconstructed mapping; FIXME: may want to use Cherenkov-to-simulated
0078     // mapping to start with, for the debugging purposes;
0079     std::map<eic::Index, const edm4eic::CherenkovParticleIDData*> rc2cherenkov;
0080     for(const auto &pid: *cherenkov)
0081       rc2cherenkov[pid.recID] = &pid;
0082
0083     //printf("Here! %d\n", mctracks->size());
0084     // Loop through all MC tracks;
0085     for(auto mctrack: *mctracks) {
0086       // FIXME: consider only primaries for now?;
0087       if (mctrack.g4Parent) continue;
0088
0089 #ifdef _USE_RECONSTRUCTED_TRACKS_
0090       // Find a matching reconstructed track;
0091       auto rctrack = mc2rc.find(mctrack.ID) == mc2rc.end() ? 0 : mc2rc[mctrack.ID];
0092       if (!rctrack) continue;
0093
0094       // Find a matching Cherenkov PID record;
0095       auto cherenkov = rc2cherenkov.find(rctrack.ID) == rc2cherenkov.end() ? 0 : rc2cherenkov[rctrack.ID];
0096 #else
0097       auto cherenkov = rc2cherenkov.find(mctrack.ID) == rc2cherenkov.end() ? 0 : rc2cherenkov[mctrack.ID];
0098 #endif
0099       if (!cherenkov) continue;
0100
0101
0102
0103       // BACK !!!;
0104       double pp = mctrack.ps.mag(), m = 0.493677;//mctrack.mass;
0105
0106
0107
0108       //continue;
0109       printf("%f %f (%4d) \n", mctrack.mass, mctrack.ps.mag(), mctrack.pdgID);
0110
0111       // Loop through all of the mass hypotheses available for this reconstructed track;
0112       {
0113     const edm4eic::CherenkovPdgHypothesis *best = 0;
0114
0115     for(unsigned iq=cherenkov->options_begin; iq<cherenkov->options_end; iq++) {
0116       const auto &option = (*options)[iq];
0117
0118       if (option.radiator != id) continue;
0119
0120       // Skip electron hypothesis; of no interest here;
0121       if (abs(option.pdg) == 11) continue;
0122
0123       if (abs(option.pdg) == _NPE_REFERENCE_) {
0124         np->Fill(option.npe);
0125
0126         //if (ofname) npvector.push_back(option.npe);
0127       } //if
0128
0129       if (!best || option.weight > best->weight) best = &option;
0130       printf("radiator %3d (pdg %5d): weight %7.2f, npe %7.2f\n",
0131          option.radiator, option.pdg, option.weight, option.npe);
0132     } //for
0133     printf("\n");
0134
0135     //if (!best) printf("@J@\n");
0136     // Check whether the true PDG got a highest score;
0137     if (!best || best->pdg != mctrack.pdgID) {
0138       //printf("@J@ %7.2f\n", best->npe);
0139       false_assignment_stat[best->npe >= 0 ? 0 : 1]++;
0140     }
0141
0142     // This assumes of course that at least one radiator was requested in juggler;
0143     double rindex = (*angles)[id].rindex, theta = (*angles)[id].theta;
0144     double argument = sqrt(pp*pp + m*m)/(rindex*pp);
0145     double thp = fabs(argument) <= 1.0 ? acos(argument) : theta;
0146
0147     (mctrack.pdgID == 321 ? th : tq)->Fill(1000*  theta);
0148     ri->Fill(rindex - 1.0);
0149     printf("<n> ~ %8.6f, <th> = %7.2f [mrad]\n", rindex - 1.0, 1000*thp);
0150
0151     {
0152       int thbin = (int)floor((1000*theta - thmin)/thstep);
0153       if (thbin >= 0 && thbin < qdim)
0154         thstat[mctrack.pdgID == 321 ? 0 : 1][thbin]++;
0155     }
0156       }
0157     } //for track
0158   } //for ev
0159
0160   printf("%3d (%3d) false out of %lld\n", false_assignment_stat[0],
0161      false_assignment_stat[1], it->GetEntries());
0162
0163   {
0164     gROOT->Reset();
0165
0166     gStyle->SetTextSize(0.02);
0167     gStyle->SetLabelSize(0.04,"xy");
0168     gStyle->SetFrameFillColor(0);
0169     gStyle->SetOptStat(0);
0170     gStyle->SetOptFit(0);
0171     gStyle->SetPadBottomMargin(0.10);
0172     gStyle->SetPadTopMargin(0.04);
0173     gStyle->SetPadRightMargin(0.05);
0174     gStyle->SetPadLeftMargin(0.12);
0175
0176     gStyle->SetStatBorderSize(0);
0177     gStyle->SetStatColor(kWhite);
0178     gStyle->SetStatFontSize(0.03);
0179     gStyle->SetStatFont(52);
0180     gStyle->SetStatW(.13);
0181     gStyle->SetFitFormat("5.2f");
0182
0183     TCanvas *cv = new TCanvas("cv", "", 0, 0, 1200, 600);
0184     cv->UseCurrentStyle();
0185     cv->SetBorderMode(0);
0186     cv->SetFrameBorderMode(0);
0187     cv->SetFrameLineColor(kWhite);
0188     cv->SetFillColor(kWhite);
0189
0190     cv->Divide(2, 1);
0191     th->SetLineWidth(2); tq->SetLineWidth(2);
0192
0193     tq->GetXaxis()->SetLabelFont(52);
0194     tq->GetYaxis()->SetLabelFont(52);
0195     tq->GetXaxis()->SetTitleFont(52);
0196     tq->GetYaxis()->SetTitleFont(52);
0197
0198     tq->GetXaxis()->SetTitle("Reconstructed Cherenkov angle, [mrad]");
0199     tq->GetXaxis()->SetTitleSize(_TSIZE_);
0200     tq->GetXaxis()->SetLabelSize(_LSIZE_);
0201     tq->GetYaxis()->SetTitle("Number of Events");
0202     tq->GetYaxis()->SetTitleSize(_TSIZE_);
0203     tq->GetYaxis()->SetLabelSize(_LSIZE_);
0204     tq->GetXaxis()->SetTitleOffset(0.90);
0205     tq->GetYaxis()->SetTitleOffset(1.10);
0206
0207     tq->GetXaxis()->SetNdivisions(408);
0208     tq->GetYaxis()->SetNdivisions(808);
0209
0210     cv->cd(1);
0211     gPad->SetGrid();
0212     tq->Draw();       tq->Fit("gaus");
0213     th->Draw("SAME"); th->Fit("gaus");
0214
0215     {
0216       // First calculate total kaon sample stat;
0217       unsigned kstat = 0, pistat = 0;//, offset = 16;
0218       for(unsigned iq=0; iq<qdim; iq++) {
0219     kstat  += thstat[0][iq];
0220     pistat += thstat[1][iq];
0221       } //for iq
0222       // This is not dramatically efficient;
0223       unsigned kaccu [qdim]; memset( kaccu, 0x00, sizeof( kaccu));
0224       unsigned piaccu[qdim]; memset(piaccu, 0x00, sizeof(piaccu));
0225       for(unsigned iq=0; iq<qdim; iq++) {
0226     kaccu [iq] = (iq ? kaccu [iq-1] : 0) + thstat[0][iq];
0227     piaccu[iq] = (iq ? piaccu[iq-1] : 0) + thstat[1][iq];
0228       } //for iq
0229       for(unsigned iq=0; iq<qdim; iq++)
0230     printf("%2d: %4d -> %5.3f   %4d -> %5.3f\n",
0231            iq, kaccu[iq], 1.*kaccu[iq]/kstat, piaccu[iq], 1.*piaccu[iq]/pistat);
0232
0233       float   eff[gdim];
0234       float   prf[gdim];
0235       for(unsigned ig=0; ig<gdim; ig++) {
0236     unsigned ibin = ig + offset;
0237     eff[ig] = 1.* kaccu[ibin]/kstat;
0238     prf[ig] = 1.*pistat/piaccu[ibin];
0239       } //for ig
0240       float  eeff[gdim] = { .001, .001, .001, .001, .001, .001, .001, .001, .001, .001, .001};
0241       float  eprf[gdim] = { .001, .001, .001, .001, .001, .001, .001, .001, .001, .001, .001};
0242
0243       auto greff = new TGraphErrors(gdim, eff, prf, eeff, eprf);
0244       greff->SetMarkerSize(1.3);
0245       greff->SetMarkerStyle(21);
0246       greff->SetMarkerColor(kBlue);
0247
0248       auto hdum = new TH1D("hdum", "", 10, 0.50, 1.00);
0249
0250       hdum->GetXaxis()->SetLabelFont(52);
0251       hdum->GetYaxis()->SetLabelFont(52);
0252       hdum->GetXaxis()->SetTitleFont(52);
0253       hdum->GetYaxis()->SetTitleFont(52);
0254
0255       hdum->GetXaxis()->SetTitle("Kaon detection efficiency, [0..1]");
0256       hdum->GetXaxis()->SetTitleSize(_TSIZE_);
0257       hdum->GetXaxis()->SetLabelSize(_LSIZE_);
0258       hdum->GetYaxis()->SetTitle("Pion rejection factor");
0259       hdum->GetYaxis()->SetTitleSize(_TSIZE_);
0260       hdum->GetYaxis()->SetLabelSize(_LSIZE_);
0261       hdum->GetXaxis()->SetTitleOffset(0.90);
0262       hdum->GetYaxis()->SetTitleOffset(1.10);
0263       hdum->SetMinimum(    0.00);
0264       hdum->SetMaximum(  200.);
0265
0266       hdum->GetXaxis()->SetNdivisions(408);
0267       hdum->GetYaxis()->SetNdivisions(808);
0268
0269       cv->cd(2);
0270       gPad->SetGrid();
0271       hdum->Draw();
0272       greff->Draw("PC");
0273     }
0274   } //if
0275 } // evaluation()