pfRICH/scripts/reco-simple.C

0001
0002
0003 void reco_simple(const char *dfname, const char *cfname = 0)
0004 {
0005   auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH-Simple");
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
0017   // Sensor active area pixelated will be pixellated NxN in digitization; 3.25mm * 400 ~ 1300mm;
0018   reco->SetSensorActiveAreaPixellation(400);
0019
0020   // [rad] (should match SPE sigma) & [ns];
0021   auto *a1 = reco->UseRadiator("Aerogel225");//,      0.0040);
0022   //auto *a1 = reco->UseRadiator("BelleIIAerogel3");
0023
0024   //reco->SetSinglePhotonTimingResolution(0.030);
0025   //reco->SetQuietMode();
0026   reco->AddHypothesis("pi+");
0027   reco->AddHypothesis(321);
0028   //reco->IgnoreMcTruthPhotonDirectionSeed();
0029
0030   // Mark all pads hit by "calibration" (above the QE curve) photons as "useless";
0031   reco->AddBlackoutRadiator("QuartzWindow");
0032   reco->AddBlackoutRadiator("Acrylic");
0033   reco->AddBlackoutRadiator("GasVolume");
0034   // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits;
0035   reco->SetBlackoutBlowupValue(3);
0036
0037   auto hmatch = new TH1D("hmatch", "PID evaluation correctness",       3,    0,      3);
0038   auto hthtr1 = new TH1D("thtr1",  "Cherenkov angle (track)",        200,  220,    320);
0039
0040   reco->UseAutomaticCalibration();
0041   // This call is mandatory; second argument: statistics (default: all events);
0042   reco->PerformCalibration(200);
0043   {
0044     CherenkovEvent *event;
0045
0046     // Event loop;
0047     while((event = reco->GetNextEvent())) {
0048       for(auto mcparticle: event->ChargedParticles()) {
0049     if (!mcparticle->IsPrimary()) continue;
0050
0051     if (mcparticle->GetRecoCherenkovHitCount() <= 3) {
0052       hmatch->Fill(0.5);
0053     } else if (mcparticle->GetPDG() == mcparticle->GetRecoPdgCode()) {
0054       hmatch->Fill(1.5);
0055     } else {
0056       hmatch->Fill(2.5);
0057     } //if
0058
0059     hthtr1->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a1));
0060       } //for mcparticle
0061     } //while
0062   }
0063
0064   auto cv = new TCanvas("cv", "", 1600, 1000);
0065   cv->Divide(4, 3);
0066   cv->cd(1); reco->hthph()->Fit("gaus");
0067   cv->cd(2); reco->hccdfph()->SetMinimum(0); reco->hccdfph()->Draw();
0068   cv->cd(3); reco->hccdftr()->SetMinimum(0); reco->hccdftr()->Draw();
0069   cv->cd(4); reco->hccdfev()->SetMinimum(0); reco->hccdfev()->Draw();
0070   cv->cd(5); reco->hnpetr()->Draw();
0071   cv->cd(6); reco->hmatch()->SetMinimum(0); reco->hmatch()->Draw();
0072   cv->cd(7);       hmatch  ->SetMinimum(0);       hmatch  ->Draw();
0073   cv->cd(8); reco->hdtph()->Fit("gaus");
0074   cv->cd(9); hthtr1->Fit("gaus");
0075   //cv->cd(10); hthtr2->Fit("gaus");
0076   cv->cd(10); reco->hwl()->Draw();
0077   cv->cd(11); reco->hvtx()->Draw();
0078   cv->cd(12); reco->hri()->Draw();
0079 } // reco_simple()