pfRICH/scripts/reco-epic.C

0001
0002
0003 void reco_epic(const char *dfname, const char *cfname = 0)
0004 {
0005   auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH");
0006
0007   //
0008   // Factory configuration part;
0009   //
0010   //reco->IgnoreTimingInChiSquare();
0011   //reco->IgnorePoissonTermInChiSquare();
0012   //reco->SetSingleHitCCDFcut(0.005);
0013   //reco->RemoveAmbiguousHits();
0014   // This only affects the calibration stage;
0015   //reco->SetDefaultSinglePhotonThetaResolution(0.0040);
0016   // Sensor active area pixelated will be pixellated NxN in digitization;
0017   //reco->SetSensorActiveAreaPixellation(24);
0018   // [rad] (should match SPE sigma) & [ns];
0019   //auto *a1 = reco->UseRadiator("Aerogel225",      0.0040);
0020
0021   auto *a1 = reco->UseRadiator("BelleIIAerogel3");
0022   //auto *a1 = reco->UseRadiator("BelleIIAerogel3");
0023
0024
0025   //reco->SetSinglePhotonTimingResolution(0.030);
0026   //reco->SetQuietMode();
0027   reco->AddHypothesis("pi+");
0028   reco->AddHypothesis(321);
0029   //reco->IgnoreMcTruthPhotonDirectionSeed();
0030
0031   // Mark all pads hit by "calibration" (above the QE curve) photons as "useless";
0032   reco->AddBlackoutRadiator("QuartzWindow");
0033   reco->AddBlackoutRadiator("Acrylic");
0034   reco->AddBlackoutRadiator("GasVolume");
0035   // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits;
0036   reco->SetBlackoutBlowupValue(3);
0037
0038   reco->ImportTrackingSmearing("./database/dtheta_seed_param.reformatted.dat", "./database/dphi_seed_param.reformatted.dat");
0039
0040   auto hmatch = new TH1D("hmatch", "PID evaluation correctness",       3,    0,      3);
0041   //auto hthtr1 = new TH1D("thtr1",  "Cherenkov angle (track)",        200,  220,    320);
0042   auto hthtr1 = new TH1D("thtr1",  "Cherenkov angle (track)",        40,  270, 290);
0043   // For a dual aerogel configuration;
0044   //auto hthtr2  = new TH1D("thtr2",   "Cherenkov angle (track)",        200,  220,    320);
0045
0046   reco->UseAutomaticCalibration();
0047   // This call is mandatory; second argument: statistics (default: all events);
0048   reco->PerformCalibration(200);
0049   {
0050     CherenkovEvent *event;
0051
0052     // Event loop;
0053     while((event = reco->GetNextEvent())) {
0054       for(auto mcparticle: event->ChargedParticles()) {
0055     if (!mcparticle->IsPrimary()) continue;
0056
0057     if (mcparticle->GetRecoCherenkovHitCount() <= 3) {
0058       hmatch->Fill(0.5);
0059     } else if (mcparticle->GetPDG() == mcparticle->GetRecoPdgCode()) {
0060       hmatch->Fill(1.5);
0061     } else {
0062       hmatch->Fill(2.5);
0063     } //if
0064
0065     hthtr1->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a1));
0066     //hthtr2->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a2));
0067       } //for mcparticle
0068     } //while
0069   }
0070
0071   auto cv = new TCanvas("cv", "", 1600, 1000);
0072   cv->Divide(4, 3);
0073   cv->cd(1); reco->hthph()->Fit("gaus");
0074   cv->cd(2); reco->hccdfph()->SetMinimum(0); reco->hccdfph()->Draw();
0075   cv->cd(3); reco->hccdftr()->SetMinimum(0); reco->hccdftr()->Draw();
0076   cv->cd(4); reco->hccdfev()->SetMinimum(0); reco->hccdfev()->Draw();
0077   cv->cd(5); reco->hnpetr()->Draw();
0078   cv->cd(6); reco->hmatch()->SetMinimum(0); reco->hmatch()->Draw();
0079   cv->cd(7);       hmatch  ->SetMinimum(0);       hmatch  ->Draw();
0080   cv->cd(8); reco->hdtph()->Fit("gaus");
0081   cv->cd(9); hthtr1->Fit("gaus");
0082   //cv->cd(10); hthtr2->Fit("gaus");
0083   cv->cd(10); reco->hwl()->Draw();
0084   cv->cd(11); reco->hvtx()->Draw();
0085   cv->cd(12); reco->hri()->Draw();
0086 } // reco_epic()