pfRICH/scripts/reco-cern.C

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