pfRICH/scripts/reco-epic-LUT.C

0001 void reco_epic_LUT(const char *dfname, const char *output, const char *cfname = 0)
0002 {
0003   auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH");
0004
0005   // const char* output="test.hist.root";
0006   TFile *ofile = TFile::Open(output, "recreate");
0007
0008   TTree *tree = new TTree("T", "Tree of Quantities");
0009
0010   Int_t partPDG;
0011   Double_t partMom;
0012   Double_t partEta;
0013   Double_t partTheta;
0014   Double_t partPhi;
0015   Double_t partVertX;
0016   Double_t partVertY;
0017   Double_t partVertZ;
0018
0019   Int_t recoPDG;
0020   Int_t recoPhotons;
0021   Double_t recoTrackCherenkovAngle;
0022
0023   tree->Branch("partPDG", &partPDG, "partPDG/I");
0024   tree->Branch("partMom", &partMom, "partMom/D");
0025   tree->Branch("partEta", &partEta, "partEta/D");
0026   tree->Branch("partTheta", &partTheta, "partTheta/D");
0027   tree->Branch("partPhi", &partPhi, "partPhi/D");
0028   tree->Branch("partVertX", &partVertX, "partVertX/D");
0029   tree->Branch("partVertY", &partVertY, "partVertY/D");
0030   tree->Branch("partVertZ", &partVertZ, "partVertZ/D");
0031
0032   tree->Branch("recoPDG", &recoPDG, "recoPDG/I");
0033   tree->Branch("recoPhotons", &recoPhotons, "recoPhotons/I");
0034   tree->Branch("recoTrackCherenkovAngle", &recoTrackCherenkovAngle, "recoTrackCherenkovAngle/D");
0035
0036   // Factory configuration part;
0037   //
0038   // reco->IgnoreTimingInChiSquare();
0039   // reco->IgnorePoissonTermInChiSquare();
0040   // reco->SetSingleHitCCDFcut(0.005);
0041   // reco->RemoveAmbiguousHits();
0042   // This only affects the calibration stage;
0043   // reco->SetDefaultSinglePhotonThetaResolution(0.0040);
0044   // Sensor active area pixelated will be pixellated NxN in digitization;
0045   // reco->SetSensorActiveAreaPixellation(24);
0046   // [rad] (should match SPE sigma) & [ns];
0047   // auto *a1 = reco->UseRadiator("Aerogel225",      0.0040);
0048
0049   auto *a1 = reco->UseRadiator("BelleIIAerogel3");
0050   // auto *a1 = reco->UseRadiator("BelleIIAerogel3");
0051
0052   // reco->SetSinglePhotonTimingResolution(0.030);
0053   // reco->SetQuietMode();
0054   // reco->AddHypothesis(-11);
0055   reco->AddHypothesis(11);
0056   reco->AddHypothesis("pi+");
0057   reco->AddHypothesis(321);
0058   reco->AddHypothesis(2212);
0059   // reco->IgnoreMcTruthPhotonDirectionSeed();
0060
0061   // Mark all pads hit by "calibration" (above the QE curve) photons as "useless";
0062   reco->AddBlackoutRadiator("QuartzWindow");
0063   reco->AddBlackoutRadiator("Acrylic");
0064   reco->AddBlackoutRadiator("GasVolume");
0065   // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits;
0066   reco->SetBlackoutBlowupValue(3);
0067   reco->ImportTrackingSmearing("./database/dtheta_seed_param.reformatted.dat", "./database/dphi_seed_param.reformatted.dat");
0068
0069   reco->UseAutomaticCalibration();
0070   // This call is mandatory; second argument: statistics (default: all events);
0071   reco->PerformCalibration(200);
0072
0073   // int numElec = 0;
0074   // int numPion = 0;
0075   // int numKaon = 0;
0076   // int numProt = 0;
0077   {
0078     CherenkovEvent *event;
0079
0080     // Event loop;
0081     while ((event = reco->GetNextEvent()))
0082     {
0083       for (auto mcparticle : event->ChargedParticles())
0084       {
0085         if (!mcparticle->IsPrimary())
0086           continue;
0087
0088         partPDG = mcparticle->GetPDG();
0089         partMom = mcparticle->GetVertexMomentum().Mag();
0090         partEta = mcparticle->GetVertexMomentum().PseudoRapidity();
0091         partPhi = mcparticle->GetVertexMomentum().Phi();
0092         partTheta = mcparticle->GetVertexMomentum().Theta();
0093         partVertX = mcparticle->GetVertexPosition().X();
0094         partVertY = mcparticle->GetVertexPosition().Y();
0095         partVertZ = mcparticle->GetVertexPosition().Z();
0096
0097         recoPDG = mcparticle->GetRecoPdgCode();
0098         recoPhotons = mcparticle->GetRecoCherenkovPhotonCount();
0099         recoTrackCherenkovAngle = 1000.0 * mcparticle->GetRecoCherenkovAverageTheta(a1);
0100
0101         tree->Fill();
0102
0103       } // for mcparticle
0104     }   // while
0105   }
0106
0107   tree->Write();
0108   ofile->Close();
0109 }