File indexing completed on 2025-01-18 10:17:17
0001
0002 void reader(const char *dfname, const char *cfname, const char *dtname = 0)
0003 {
0004 #define _AEROGEL_
0005 #define _DRICH_
0006
0007
0008 #define _WAVE_LENGTH_CUTOFF_MIN_ (350.0)
0009 #define _WAVE_LENGTH_CUTOFF_MAX_ (650.0)
0010
0011
0012
0013 auto fdata = new TFile(dfname);
0014 TTree *t = (TTree*)fdata->Get("events");
0015
0016
0017
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
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 }
0031
0032 detector = (*detectors.begin()).second;
0033 }
0034
0035
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 }
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
0045
0046
0047
0048 #ifdef _AEROGEL_
0049
0050 auto ep = new TH1D("ep", "Emission Point", 100, -30.0, 30.0);
0051
0052 #else
0053
0054
0055
0056 #endif
0057 auto wl = new TH1D("wl", "Wave Length; #{Lambda} (nm)", 100, 200.0, 1000.0);
0058
0059
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
0065
0066
0067 gas ->m_AverageRefractiveIndex = gas ->n();
0068 aerogel->m_AverageRefractiveIndex = 1.020;
0069
0070
0071
0072 aerogel->SetGaussianSmearing(0.001);
0073
0074
0075
0076
0077 aerogel->SetTrajectoryBinCount(2);
0078
0079
0080
0081
0082 int pdg;
0083 int q;
0084
0085 auto event = new CherenkovEvent();
0086
0087
0088
0089
0090
0091
0092
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
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116 int evtcounter =0;
0117 while(myReader.Next()){
0118 printf("#################\n");
0119
0120
0121
0122 for(auto && mcpart : *mcparts){
0123
0124 if(mcpart.parents_begin!=mcpart.daughters_begin) continue;
0125 {
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135 }
0136 int phcounter =0;
0137 std::vector<OpticalPhoton*> photons;
0138 for(auto &&hit:*hits){
0139
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
0143 auto xx = hit.position.x; auto yy = hit.position.y;
0144
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
0155 photon->SetPhotonDetector(detector->m_PhotonDetectors[0]);
0156 photon->SetDetected(true); phcounter+=1;
0157
0158 photon->SetVolumeCopy(hit.cellID & detector->GetReadoutCellMask());
0159
0160 photons.push_back(photon);
0161 }
0162 printf("Set True Phot: %d\n",phcounter);
0163 auto particle = new ChargedParticle(mcpart.PDG);
0164 event->AddChargedParticle(particle);
0165 gas->ResetLocations();
0166
0167 particle->StartRadiatorHistory(std::make_pair(gas, new RadiatorHistory()));
0168 {
0169
0170
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
0174 TVector3 from, to;
0175 gas->GetFrontSide(0)->GetCrossing(x0, n0, &from);
0176 gas->GetRearSide (0)->GetCrossing(x0, n0, &to);
0177
0178
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());
0183 }
0184 {
0185 CherenkovPID pid;
0186
0187
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
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
0204 }
0205 }
0206 printf("&&&&&&&&&&&&&&\n");
0207 }
0208 cout<<"#################### &&&&&&&&&&&&& "<< evtcounter <<" &&&&&&&&&&&& ######################"<<endl;
0209 evtcounter++;
0210 }
0211
0212 np->Draw();
0213 }