obsolete/ak/e-eval.C

0001
0002 #define _MASS_      (0.140)
0003 #define _PDG_         211
0004 //#define _MASS_       (0.494)
0005 //#define _PDG_         321
0006
0007 //#define _AEROGEL_ ("Aerogel155")
0008 #define _AEROGEL_ ("Aerogel225")
0009 //#define _AEROGEL_ ("BelleIIAerogel1")
0010
0011 // Kind of digitization; of course only one of the two must be used;
0012 //#define _GAUSSIAN_SMEARING_            1.0
0013 #define _DISCRETE_PIXEL_PITCH_         4.0
0014
0015 // [rad] & [mm];
0016 //#define _PHOTON_DETECTION_WINDOW_    0.012
0017 #define _PHOTON_DETECTION_WINDOW_    0.004
0018 #define _AVERAGE_ATTENUATION_LENGTH_   8.4
0019
0020 void e_eval(const char *dfname, const char *cfname = 0)
0021 {
0022
0023   auto fcfg  = new TFile(cfname ? cfname : dfname);
0024   auto geometry = dynamic_cast<CherenkovDetectorCollection*>(fcfg->Get("CherenkovDetectorCollection"));
0025   auto fdata = new TFile(dfname);
0026   TTree *t = dynamic_cast<TTree*>(fdata->Get("t"));
0027   auto event = new CherenkovEvent();
0028   t->SetBranchAddress("e", &event);
0029   auto ri = new TH1D("ri", "1.0 - Refractive Index",  500, 0.015,  0.055);
0030   auto qh = new TH1D("qh", "Cherenkov angle (SPE)",   100,   -30,     30);
0031   auto th = new TH1D("th", "Cherenkov angle (track)", 100,   -10,     10);
0032   auto np = new TH1D("np", "Photon count (pion?)",     50,     0,     50);
0033   auto nq = new TH1D("nq", "Photon count (radiator)",  50,     0,     50);
0034   auto wl = new TH1D("wl", "Wave Length",             200, 150.0,  800.0);
0035   auto tt = new TH2D("tt", "Cherenkov angle vs time",  50,  -0.2, 0.6, 50, -100,  100);
0036   auto z0 = new TH1D("z0", "True emission point",      80,   -40,     40);
0037   auto al = new TH1D("al", "True attenuation length", 100,     0,    100);
0038
0039   int nEvents = t->GetEntries();
0040
0041   auto rich = geometry->GetDetector("pfRICH");
0042   auto aerogel = rich->GetRadiator(_AEROGEL_);
0043
0044   // First calculate average refractive index and theta for all detected photons
0045   // for every radiator;
0046   {
0047     for(unsigned ev=0; ev<nEvents; ev++) {
0048       t->GetEntry(ev);
0049
0050       for(auto particle: event->ChargedParticles()) {
0051     if (particle->GetPDG() != _PDG_) continue;
0052
0053     for(auto rhistory: particle->GetRadiatorHistory()) {
0054       auto radiator = particle->GetRadiator(rhistory);
0055       auto history  = particle->GetHistory (rhistory);
0056
0057       if (history->StepCount() >= 2) {
0058         auto from = history->GetStep(0), to = history->GetStep(history->StepCount()-1);
0059         TVector3 ntrack = (from->GetMomentum() + to->GetMomentum()).Unit();
0060
0061         for(auto photon: history->Photons()) {
0062           if (!photon->WasDetected() /*|| photon->m_ReflectionPoints.size() != 1*/) continue;
0063
0064           radiator->m_Stat++;
0065           radiator->m_AverageTheta += acos(ntrack.Dot(photon->GetVertexMomentum().Unit()));
0066           radiator->m_AverageRefractiveIndex += photon->GetVertexRefractiveIndex();
0067           radiator->m_AverageVertexPosition  += photon->GetVertexPosition();
0068         } //for photon
0069       } //if
0070     } //for radiator
0071       } //for particle
0072     } //for ev
0073
0074     for(auto rptr: rich->Radiators()) {
0075       printf("%s\n", rptr.first.Data());
0076       auto radiator = rptr.second;
0077
0078       if (radiator->m_Stat) {
0079     radiator->m_AverageTheta           /= radiator->m_Stat;
0080     radiator->m_AverageRefractiveIndex /= radiator->m_Stat;
0081     radiator->m_AverageVertexPosition  *= 1.0/radiator->m_Stat;
0082       } //if
0083
0084       printf("%5d photons, <theta> = %7.2f mrad; <n> = %8.5f\n", radiator->m_Stat,
0085          1000*radiator->m_AverageTheta, radiator->m_AverageRefractiveIndex);
0086     } //for dt..radiator
0087   }
0088
0089   aerogel->SetTrajectoryBinCount(10);
0090   //aerogel->SetUniformSmearing(_PHOTON_DETECTION_WINDOW_);
0091   aerogel->SetGaussianSmearing(_PHOTON_DETECTION_WINDOW_);
0092   //aerogel->SetReferenceRefractiveIndex(_AVERAGE_REFRACTIVE_INDEX_);
0093   aerogel->SetReferenceAttenuationLength(_AVERAGE_ATTENUATION_LENGTH_);
0094
0095   unsigned false_assignment_stat = 0;
0096
0097   printf("%d\n", nEvents);
0098   for(unsigned ev=0; ev<nEvents; ev++) {
0099     t->GetEntry(ev);
0100
0101     for(auto particle: event->ChargedParticles()) {
0102       TVector3 mid;
0103
0104       for(auto rhistory: particle->GetRadiatorHistory()) {
0105     auto history  = particle->GetHistory (rhistory);
0106     auto radiator = particle->GetRadiator(rhistory);
0107     const unsigned zdim = radiator->GetTrajectoryBinCount();
0108
0109     radiator->ResetLocations();
0110     if (history->StepCount() >= 2) {
0111       auto fstep = history->GetStep(0), lstep = history->GetStep(history->StepCount()-1);
0112
0113       auto s1 = radiator->GetFrontSide(0);//ifsec);
0114       auto s2 = radiator->GetRearSide (0);//ilsec);
0115
0116       // FIXME: check what's wrong with 's2';
0117       if (s1 && s2) {
0118         TVector3 from, to;//, p0;
0119
0120         auto x0 = fstep->GetPosition();
0121         auto n0 = fstep->GetDirection();
0122         // Go backwards and ignore surface orientation mismatch;
0123         bool b1 = s1->GetCrossing(x0, -1*n0, &from, false);
0124
0125         auto x1 = lstep->GetPosition();
0126         auto n1 = lstep->GetDirection();
0127         bool b2 = s2->GetCrossing(x1,    n1, &to);
0128
0129         assert(b1 && b2);
0130
0131         //auto fstep = history->GetStep(0), lstep = history->GetStep(stCount-1);
0132         //auto from = history->GetStep(0), to = history->GetStep(history->StepCount()-1);
0133         TVector3 p0 = 0.5*(fstep->GetMomentum() + lstep->GetMomentum());
0134
0135         const unsigned zbins = radiator->GetTrajectoryBinCount();
0136         TVector3 nn = (to - from).Unit();
0137         // FIXME: hardcoded;
0138         from += (0.010)*nn;
0139         to   -= (0.010)*nn;
0140
0141         if (radiator == aerogel) mid = 0.5*(from + to);
0142
0143         // FIXME: assume a straight line to the moment;
0144         auto span = to - from;
0145         double tlen = span.Mag();
0146         double step = tlen / zbins;
0147         for(unsigned iq=0; iq<zbins+1; iq++) {
0148           radiator->AddLocation(from + iq*step*nn, p0);
0149           //printf("(2) %f %f %f\n",
0150           //       (from + iq*step*nn).x(),
0151           //       (from + iq*step*nn).y(),
0152           //       (from + iq*step*nn).z());
0153         } //for iq
0154       } //if
0155       //TVector3 ptnx = (to->GetPosition() - from->GetPosition()).Unit();
0156
0157       //+auto vstart = from->GetPosition() + 0.001*ptnx, vend = to->GetPosition() - 0.001*ptnx;
0158       //auto vstart = from->GetPosition() + 0.001*ptnx, vend = vstart + 30.0*ptnx;
0159       //double vlen = (vend - vstart).Mag(), step = vlen/zdim;
0160       //if (radiator == aerogel) printf("%f %f -> %f\n", from->GetPosition().Z(), to->GetPosition().Z(), vlen);
0161
0162       //for(unsigned iq=0; iq<zdim+1; iq++)
0163       //radiator->AddLocation(vstart + iq*step*ptnx, p);
0164     } //if
0165       } //for rhistory
0166
0167       // Now that all internal track-level structures are populated, run IRT code;
0168       CherenkovPID pid;
0169
0170       // Consider just pi/K case for now;
0171       pid.AddMassHypothesis(0.140);
0172       pid.AddMassHypothesis(0.494);
0173
0174       for(auto rhistory: particle->GetRadiatorHistory())
0175     for(auto photon: particle->GetHistory(rhistory)->Photons()) {
0176       TVector3 phx = photon->GetDetectionPosition();
0177
0178 #ifdef _GAUSSIAN_SMEARING_
0179       phx += TVector3(gRandom->Gaus(0.0, _GAUSSIAN_SMEARING_), gRandom->Gaus(0.0, _GAUSSIAN_SMEARING_), 0.0);
0180 #endif
0181 #ifdef _DISCRETE_PIXEL_PITCH_
0182       double pitch = _DISCRETE_PIXEL_PITCH_;
0183       phx = TVector3(pitch*rint(phx.X()/pitch), pitch*rint(phx.Y()/pitch), phx.Z());
0184 #endif
0185       photon->SetDetectionPosition(phx);
0186     } //for rhistory..photon
0187
0188       // FIXME: use true photon 3D direction at birth; otherwise conical mirror case
0189       // is problematic; do it better later;
0190       particle->PIDReconstruction(pid, true);
0191       {
0192     auto pion = pid.GetHypothesis(0), kaon = pid.GetHypothesis(1);
0193     double wt0 = pion->GetWeight(aerogel), wt1 = kaon->GetWeight(aerogel);
0194
0195     printf("%4d -> %10.3f (%10.3f) vs %10.3f (%10.3f) ...  %3d %d\n",
0196            ev, wt0, pion->GetNpe(aerogel), wt1, kaon->GetNpe(aerogel),
0197            particle->GetPDG(), wt0 > wt1);
0198     np->Fill(pion->GetNpe(aerogel));
0199
0200     if (wt0 <= wt1) false_assignment_stat++;
0201       }
0202
0203       {
0204     unsigned npe = 0;
0205     double dtheta = 0.0;
0206     //double ridx = 0.0;
0207     double alsum = 0.0;
0208     //double wlen = 0.0;
0209     double wtsum = 0.0;
0210     double dth = aerogel->GetSmearing();
0211
0212     for(auto rhistory: particle->GetRadiatorHistory()) {
0213       auto radiator = particle->GetRadiator(rhistory);
0214       if (radiator != aerogel) continue;
0215
0216       auto history  = particle->GetHistory (rhistory);
0217
0218       // This loop goes only over the photons which belong to a given radiator;
0219       for(auto photon: history->Photons()) {
0220         bool selected = false;
0221         // Check whether this photon was selected by at least one mass hypothesis;
0222         for(auto entry: photon->_m_Selected)
0223           if (entry.second == aerogel) {
0224         selected = true;
0225         break;
0226           } //if
0227
0228         if (selected) {
0229           npe++;
0230           double wt = 1.0;//photon->_m_PDF[aerogel].GetWeight();//1.0;//fabs(sin(photon->m_Phi[radiator]));
0231           //+wtsum += wt;
0232           wtsum += 1.0;
0233           {
0234         double m = _MASS_, pp = photon->GetVertexParentMomentum().Mag();
0235         double thp = acos(sqrt(pp*pp + m*m)/(radiator->m_AverageRefractiveIndex*pp));
0236         //+dtheta += wt*photon->_m_PDF[aerogel].GetAverage(thp - 3*dth, thp + 3*dth) - thp;
0237         {
0238           double qwtsum = 0.0, qthavg = 0.0;
0239
0240           for(auto member: photon->_m_PDF[aerogel].GetMembers()) {
0241             qwtsum += member->GetWeight();
0242             qthavg += member->GetWeight()*member->GetAverage();
0243           } // for member
0244
0245           qthavg /= qwtsum;
0246           dtheta += qthavg - thp;
0247           qh->Fill(1000*(qthavg - thp));//(member->GetAverage() - thp), member->GetWeight());
0248           //qh->Fill(1000*(member->GetAverage() - thp), member->GetWeight());
0249         }
0250
0251         double len = (photon->GetDetectionPosition() - mid).Mag();
0252         double beta = 1/sqrt(1.0 + pow(m/pp, 2)), lspeed = 299.792458, velocity = beta*lspeed;
0253         double te = len / velocity;
0254         //printf("%f\n", photon->GetDetectionTime() - te);
0255         tt->Fill(photon->GetDetectionTime() - te - 4.25, 1000*(photon->_m_PDF[aerogel].GetAverage() - thp));
0256           }
0257           //printf("%f\n", photon->GetVertexRefractiveIndex());//ridx);
0258           //ridx  +=    photon->GetVertexRefractiveIndex();
0259           ri->Fill(photon->GetVertexRefractiveIndex() - 1.0);
0260           // FIXME: hardcoded;
0261           //wlen  +=    1239.8/(photon->GetVertexMomentum().Mag());
0262           wl->Fill(1239.8/(photon->GetVertexMomentum().Mag()));
0263
0264           //printf("%f\n", photon->GetVertexPosition().Z());
0265           z0->Fill(photon->GetVertexPosition().Z() - 1185.5);
0266           //if (1239.8/(photon->GetVertexMomentum().Mag()) > 400)
0267           alsum += photon->GetVertexAttenuationLength();
0268         } //if
0269       } //for photon
0270     } //for rhistory
0271
0272     nq->Fill(npe);
0273     if (wtsum) th->Fill(1000.*(dtheta / wtsum));
0274     if (npe) {
0275       //ri->Fill((ridx  / npe) - 1.0);
0276       //wl->Fill( wlen  / npe);
0277       al->Fill( alsum / npe);
0278     } //if
0279       }
0280     } //for particle
0281   } //for ev
0282
0283   printf("%d false out of %lld\n", false_assignment_stat, t->GetEntries());
0284
0285   auto cv = new TCanvas("cv", "", 1700, 800);
0286   cv->Divide(4, 2);
0287   cv->cd(1); np->Draw();
0288   cv->cd(2); nq->Draw();
0289   z0->GetXaxis()->SetTitle("Local coordinate in the aerogel block, [mm]");
0290   z0->GetYaxis()->SetTitle("Detected photon count");
0291   cv->cd(3); z0->Draw();//wl->Draw();
0292   cv->cd(4); th->Fit("gaus");
0293   cv->cd(5); ri->Draw();
0294   cv->cd(6); wl->Draw();//tt->Draw("COLZ");
0295   cv->cd(7); qh->Fit("gaus");//Draw();
0296   cv->cd(8); al->Draw();
0297 } // e_eval()