irt/scripts/reader.C

0001
0002 void reader(const char *dfname, const char *cfname, const char *dtname = 0)
0003 {
0004 #define _AEROGEL_
0005 #define _DRICH_
0006
0007 // Optionally: mimic low wave length cutoff and average QE x Geometric sensor efficiency;
0008 #define _WAVE_LENGTH_CUTOFF_MIN_ (350.0)
0009 #define _WAVE_LENGTH_CUTOFF_MAX_ (650.0)
0010   //#define _AVERAGE_PDE_            ( 0.30)
0011
0012   // .root file with event tree;
0013   auto fdata = new TFile(dfname);
0014   TTree *t = (TTree*)fdata->Get("events");
0015   //t->SetMakeClass(1);
0016   //t->SetBranchStatus("*", 0);
0017   // .root file with the IRT configuration;
0018   auto fcfg  = new TFile(cfname);
0019   CherenkovDetector *detector = 0;
0020   auto geometry = dynamic_cast<CherenkovDetectorCollection*>(fcfg->Get("CherenkovDetectorCollection"));
0021
0022   if (dtname)
0023     detector = geometry->GetDetector(dtname);
0024   else {
0025     // Assume a single detector (PFRICH or DRICH);
0026     auto &detectors = geometry->GetDetectors();
0027     if (detectors.size() != 1) {
0028       printf("More than one detector in the provided IRT geometry config .root file!\n");
0029       exit(0);
0030     } //if
0031
0032     detector = (*detectors.begin()).second;
0033   } //if
0034
0035   // Either this or that way, the detector should be available;
0036   if (!detector) {
0037     printf("Was not able to find a valid Cherenkov detector in the provided IRT geometry config .root file!\n");
0038     exit(0);
0039   } //if
0040
0041   auto ed = new TH2D("ed","ed;x[cm];y[cm]",1000,-200,200, 1000,-200,200);
0042   auto nq = new TH1D("nq", "Photon count",            50,      0,    100);
0043   auto np = new TH1D("np", "Photon count",            50,      0,    100);
0044   //auto fi = new TH1D("fi", "Cherenkov phi angle",      30,-180.0,  180.0);
0045   //auto xi = new TH1D("xi", "Chi square CDFC (photon)",100,   0.0,    1.0);
0046   //auto wi = new TH1D("wi", "Chi square CDFC (track)",  20,   0.0,    1.0);
0047   //fi->SetMinimum(0); xi->SetMinimum(0); wi->SetMinimum(0);
0048 #ifdef _AEROGEL_
0049   //auto pl = new TH1D("pl", "Radiator path Length",    100,   0.0,   50.0);
0050   auto ep = new TH1D("ep", "Emission Point",          100, -30.0,   30.0);
0051   //auto ri = new TH1D("ri", "1.0 - Refractive Index",  100, 0.015,  0.025);
0052 #else
0053   //auto pl = new TH1D("pl", "Radiator path Length",    100,   0.0, 2000.0);
0054   //auto ep = new TH1D("ep", "Emission Point",          100, -70.0,   70.0);
0055   //auto ri = new TH1D("ri", "1.0 - Refractive Index",  100, 0.015,  0.025);
0056 #endif
0057   auto wl = new TH1D("wl", "Wave Length; #{Lambda} (nm)",             100, 200.0, 1000.0);
0058   //auto th = new TH1D("th", "Cherenkov angle",         100, -10,  10);
0059   //auto tq = new TH1D("tq", "Average Cherenkov angle", 100, -10,  10);
0060   auto p = new TH1D("p","Mom; p(GeV/c)",50,49.,51.);
0061
0062   auto gas      = detector->GetRadiator("GasVolume");
0063   auto aerogel  = detector->GetRadiator("Aerogel");
0064   //auto acrylic  = detector->GetRadiator("Filter");
0065   // Assume the reference value was close enough in PFRICH_geo.cpp; since QE was not accounted,
0066   // this may not be true;
0067   gas    ->m_AverageRefractiveIndex = gas    ->n();
0068   aerogel->m_AverageRefractiveIndex = 1.020;//aerogel->n();
0069   //acrylic->m_AverageRefractiveIndex = acrylic->n();
0070
0071   //#ifdef _DRICH_
0072   aerogel->SetGaussianSmearing(0.001);
0073   //#else
0074   //aerogel->SetUniformSmearing(0.005);
0075   //#endif
0076   // Be aware, that AddLocations() part should take this into account;
0077   aerogel->SetTrajectoryBinCount(2);
0078   // This may be bogus for a blob-like operation mode;
0079   //gas    ->SetUniformSmearing(0.003);
0080
0081
0082   int pdg;
0083   int q;
0084   // TTree interface variable;
0085   auto event = new CherenkovEvent();
0086
0087   // Use MC truth particles, and deal with just pfRICH hits here; however the interface
0088   // should work for combinations like pfRICH+DIRC, eventually;
0089   //edm4hep::MCParticleData  *MCParticle;     //= new std::vector<edm4hep::MCParticleData> ();
0090   //std::vector<edm4hep::SimTrackerHitData> *hits   = new std::vector<edm4hep::SimTrackerHitData>();
0091   //t->SetBranchAddress("MCParticles",&MCParticle);
0092   //t->SetBranchStatus("MCParticles.PDG",1);
0093   TTreeReader myReader("events",fdata);
0094   TTreeReaderValue <std::vector<edm4hep::MCParticleData>> mcparts(myReader,"MCParticles");
0095   TTreeReaderValue <std::vector<edm4hep::SimTrackerHitData>> hits(myReader,"DRICHHits");
0096 /*
0097   t->SetBranchAddress("MCParticles.PDG", &pdg);
0098   t->SetBranchAddress("MCParticles.charge", &q);
0099  // {
0100    // TString hname; hname.Form("%sHits", detector->GetName());
0101    // t->SetBranchAddress(hname,   &hits);
0102  // }
0103   int myentries = t->GetEntries();
0104   printf("Entries: %d\n",myentries);
0105   printf("Here!\n");
0106   // Loop through all events;
0107   //unsigned false_assignment_stat = 0;
0108   for(int ev=0; ev<myentries; ev++) {
0109     t->GetEntry(ev);
0110     //int Size = MCParticle->size();
0111     cout<<"EV: "<<ev<<" Charge:  "<<q<<endl;
0112     //printf("%d\n", tsize);
0113
0114   }//ev
0115 */
0116   int evtcounter =0;
0117   while(myReader.Next()){
0118     printf("#################\n");
0119     //evtcounter++;
0120     //cout<<"MC Size: "<<mcparts->size()<<endl;
0121     //cout<<"Hit Size: "<<hits->size()<<endl;
0122     for(auto  && mcpart : *mcparts){
0123        //cout<<val.size()<<endl;
0124       if(mcpart.parents_begin!=mcpart.daughters_begin) continue;
0125       {
0126          /*if(mcpart.PDG == 22){
0127            double pmag = TMath::Sqrt(TMath::Power(mcpart.momentum.x,2)+ TMath::Power(mcpart.momentum.y,2) + TMath::Power(mcpart.momentum.z,2));
0128            cout<<"EV: "<< evtcounter<< " PDG: "<<mcpart.PDG<<" ep x: "<<pmag<<endl;
0129            p->Fill(pmag);
0130            double wave_length = 1239.84/(1E9*pmag);
0131            wl->Fill(wave_length);
0132            cout<<mcpart.vertex.z<<" "<<mcpart.endpoint.z<<endl;
0133          }*/
0134
0135       }
0136       int phcounter =0;
0137       std::vector<OpticalPhoton*> photons;
0138       for(auto &&hit:*hits){
0139         //cout<<"hitsx   "<<hit.momentum.x<<endl;
0140         double pmag = TMath::Sqrt(TMath::Power(hit.momentum.x,2)+ TMath::Power(hit.momentum.y,2) + TMath::Power(hit.momentum.z,2));
0141         double wave_length = 1239.84/(1E9*pmag);
0142         //cout<<"Lambda: "<<wave_length<<endl;
0143         auto xx = hit.position.x; auto yy = hit.position.y;
0144         //ed->Fill(xx/10,yy/10);
0145         wl->Fill(wave_length);
0146         auto photon = new OpticalPhoton();
0147         {
0148           auto x = hit.position.x;
0149           auto y = hit.position.y;
0150           auto z = hit.position.z;
0151           printf("Recorded Hit Poistion ------> %f %f %f\n",x,y,z);
0152           photon->SetDetectionPosition(TVector3(x, y, z));
0153         }
0154         // A single photodetector type is used;
0155         photon->SetPhotonDetector(detector->m_PhotonDetectors[0]);  //?
0156         photon->SetDetected(true); phcounter+=1;
0157         // Get cell index; mask out everything apart from {module,sector};
0158         photon->SetVolumeCopy(hit.cellID & detector->GetReadoutCellMask());
0159         //photon->SetVolumeCopy(hit.cellID & detector->GetReadoutCellMask());
0160         photons.push_back(photon);
0161       }//hit
0162       printf("Set True Phot: %d\n",phcounter);
0163       auto particle = new ChargedParticle(mcpart.PDG);
0164       event->AddChargedParticle(particle);
0165       gas->ResetLocations();
0166       // Create a fake (empty) history; then track locations at the gas boundaries;
0167       particle->StartRadiatorHistory(std::make_pair(gas, new RadiatorHistory()));
0168       {
0169         // FIXME: need it not at vertex, but in the radiator; as coded here, this can
0170         // hardly work once the magnetic field is turned on;
0171         auto x0 = TVector3(mcpart.vertex.x, mcpart.vertex.y, mcpart.vertex.z), p0 = TVector3(mcpart.momentum.x, mcpart.momentum.y, mcpart.momentum.z), n0 = p0.Unit();
0172         printf("Momentum: %0.2f\n",p0.Mag());
0173         // So, give the algorithm gas surface boundaries as encoded in PFRICH_geo.cpp;
0174         TVector3 from, to;
0175         gas->GetFrontSide(0)->GetCrossing(x0, n0, &from);
0176         gas->GetRearSide (0)->GetCrossing(x0, n0, &to);
0177
0178         // Move the points a bit inwards;
0179         TVector3 nn = (to - from).Unit(); from += (0.010)*nn; to -= (0.010)*nn;
0180         gas->AddLocation(from, p0);
0181         gas->AddLocation(  to, p0);
0182         printf("@@@ %f %f\n", from.z(), to.z());// - from.z());
0183       }
0184       {
0185         CherenkovPID pid;
0186
0187         // Consider just pi/K case for now;
0188         pid.AddMassHypothesis(0.140);
0189         pid.AddMassHypothesis(0.494);
0190
0191         printf("Entering PID Rec:%zu\n",photons.size());
0192         for(auto photon : photons) particle->FindRadiatorHistory(gas)->AddOpticalPhoton(photon);
0193         particle->PIDReconstruction(pid);
0194         {
0195           auto pion = pid.GetHypothesis(0), kaon = pid.GetHypothesis(1);
0196           double wt0 = pion->GetWeight(gas), wt1 = kaon->GetWeight(gas);
0197
0198           //th->Fill(pion->GetTheta(gas));
0199
0200           printf("%10.3f (%10.3f) vs %10.3f (%10.3f) ...  %3d %d\n",
0201                  wt0, pion->GetNpe(gas), wt1, kaon->GetNpe(gas), particle->GetPDG(), wt0 > wt1);
0202
0203           //if (wt0 <= wt1) false_assignment_stat++;
0204         }
0205       }
0206       printf("&&&&&&&&&&&&&&\n");
0207     }//mctrack
0208     cout<<"#################### &&&&&&&&&&&&&   "<< evtcounter   <<"   &&&&&&&&&&&& ######################"<<endl;
0209     evtcounter++;
0210   }//ev
0211   //evtcounter++;
0212   np->Draw();
0213 } // reader()