drich/script/dRICHAna.cpp

0001 #include "edm4hep/MCParticleCollection.h"
0002 #include "edm4eic/CherenkovParticleIDCollection.h"
0003 #include "podio/ROOTReader.h"
0004 #include "podio/Frame.h"
0005 #include "TH1F.h"
0006 #include "TString.h"
0007 using namespace std;
0008
0009
0010 void extractSPEres(const char* filename, const char* outname, const char* outdir, int radiator) {
0011
0012   double thlow, thhigh;
0013   int nbins;
0014   if (radiator == 0) {
0015     thlow = 150;
0016     thhigh = 220;
0017     nbins = 200;
0018   } else {
0019     thlow = 20;
0020     thhigh = 60;
0021     nbins = 100;
0022   }
0023
0024   TH1F* hSingleThetaError = new TH1F("hSingleThetaError", "", 100, -10, 10);
0025   TH1F* hSingleTheta = new TH1F("hSingleTheta", "", nbins, thlow, thhigh);
0026   TH1F* hnPhotons = new TH1F("hnPhotons", "", 60, -0.5, 60.5);
0027
0028   podio::ROOTReader reader;
0029   reader.openFile(filename);
0030
0031   int nev = reader.getEntries("events");
0032   double nThrown = 0;
0033   double ndRICHDet = 0;
0034   // event loop
0035   for (int i = 0; i < nev; i++) {
0036     const auto event = podio::Frame(reader.readNextEntry("events"));
0037
0038     std::string pidCollection;
0039     double n;
0040     if (radiator == 0) {
0041       pidCollection = "DRICHAerogelIrtCherenkovParticleID";
0042       n = 1.019;
0043     } else {
0044       pidCollection = "DRICHGasIrtCherenkovParticleID";
0045       n = 1.00076;
0046     }
0047
0048     auto& dRichCherenkov = event.get<edm4eic::CherenkovParticleIDCollection>(pidCollection);
0049     auto& MCParticles = event.get<edm4hep::MCParticleCollection>("MCParticles");
0050
0051     double px, py, pz, p, mass;
0052     double betaTrue;
0053     // get true momentum from thrown particle
0054     if (MCParticles.isValid()) {
0055       px = MCParticles[0].getMomentum().x;
0056       py = MCParticles[0].getMomentum().y;
0057       pz = MCParticles[0].getMomentum().z;
0058       p = sqrt(px*px+py*py+pz*pz);
0059       mass = MCParticles[0].getMass();
0060       betaTrue = p/(sqrt(p*p+mass*mass));
0061     } else {
0062       cout << "Error: no thrown particles" << endl;
0063       continue;
0064     }
0065
0066     nThrown += 1.;
0067     if (dRichCherenkov.isValid()) {
0068       double chExpected = acos(1/(n*betaTrue))*1000;
0069
0070       for (const auto& pid : dRichCherenkov) {
0071         auto thetaPhi = pid.getThetaPhiPhotons();
0072
0073         int nPhotons = pid.getNpe();
0074         if (nPhotons == 0) {
0075           // if no photons, consider this to be missed
0076           continue;
0077         }
0078         if (nPhotons > 5) {
0079           ndRICHDet += 1.; // if > 5 photons, consider this to be accepted
0080         }
0081
0082         hnPhotons->Fill(nPhotons);
0083         for (const auto& theta : thetaPhi) {
0084           hSingleThetaError->Fill(abs(theta[0]*1000 - chExpected));
0085           hSingleTheta->Fill(theta[0]*1000);
0086         }
0087       }
0088     }
0089   }
0090
0091   TString outname_wdir = TString(outdir) + TString(outname);
0092   FILE *outfile = fopen(outname_wdir.Data(),"w");
0093
0094
0095   fprintf(outfile, "%lf %lf %lf %lf \n",
0096       hnPhotons->GetMean(), hSingleTheta->GetMean(),
0097       hSingleThetaError->GetMean(), // MAE
0098       ndRICHDet/nThrown);
0099
0100   return;
0101 }
0102
0103
0104 int main(int argc, char* argv[]) {
0105   if (argc < 2) {
0106     cout << "usage: dRICHAna [filename] [outputname (txt)] [output dir] [radiator: 0 - aerogel, 1 - gas] \n";
0107     return 1;
0108   }
0109   int rad;
0110   stringstream s(argv[4]);
0111   s >> rad;
0112
0113   extractSPEres(argv[1], argv[2], argv[3], rad);
0114   return 0;
0115 }