Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-17 09:57:14 

0001 void readTree(const char *input, const char *output, const char *output_txt)
0002 {
0003 
0004       cout << "running readTree.C" << endl;
0005 
0006       TFile *fa = new TFile(input);
0007 
0008       TTree *ta = (TTree *)fa->Get("T");
0009 
0010       TFile *ofile = TFile::Open(output, "recreate");
0011 
0012       Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount;
0013       Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle;
0014       Double_t zero = pow(10., -50.);
0015 
0016       ta->SetBranchAddress("partPDG", &partPDG);
0017       ta->SetBranchAddress("partMom", &partMom);
0018       ta->SetBranchAddress("partEta", &partEta);
0019       ta->SetBranchAddress("partTheta", &partTheta);
0020       ta->SetBranchAddress("partPhi", &partPhi);
0021       ta->SetBranchAddress("partVertX", &partVertX);
0022       ta->SetBranchAddress("partVertY", &partVertY);
0023       ta->SetBranchAddress("partVertZ", &partVertZ);
0024       ta->SetBranchAddress("recoPDG", &recoPDG);
0025       ta->SetBranchAddress("recoPhotons", &recoPhotons);
0026       ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle);
0027 
0028       TH1D *h_p = new TH1D("p", "", 37, 0.4, 15.2);
0029       TH1D *h_t = new TH1D("theta", "", 20, 2.65, 3.1);
0030       TH1D *h_phi = new TH1D("phi", "", 120, -1.0 * TMath::Pi(), TMath::Pi());
0031 
0032       static Float_t fee[37][20][120], fepi[37][20][120], fek[37][20][120], fep[37][20][120], feb[37][20][120];
0033       static Float_t fpie[37][20][120], fpipi[37][20][120], fpik[37][20][120], fpip[37][20][120], fpib[37][20][120];
0034       static Float_t fke[37][20][120], fkpi[37][20][120], fkk[37][20][120], fkp[37][20][120], fkb[37][20][120];
0035       static Float_t fpe[37][20][120], fppi[37][20][120], fpk[37][20][120], fpp[37][20][120], fpb[37][20][120];
0036 
0037       for (int pbin = 0; pbin < 37; pbin++)
0038       {
0039             for (int tbin = 0; tbin < 20; tbin++)
0040             {
0041                   for (int phibin = 0; phibin < 120; phibin++)
0042                   {
0043                         fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0.;
0044                         fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0.;
0045                         fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0.;
0046                         fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0.;
0047                   }
0048             }
0049       }
0050 
0051       Int_t nEntries = (Int_t)ta->GetEntries();
0052       for (int i = 0; i < nEntries; i++)
0053       {
0054             if (i % 1000000 == 0)
0055                   cout << "Processing Event " << i << endl;
0056             ta->GetEntry(i);
0057 
0058             Int_t pbin = h_p->FindBin(partMom) - 1;
0059             Int_t tbin = h_t->FindBin(partTheta) - 1;
0060             Int_t phibin = h_phi->FindBin(partPhi) - 1;
0061 
0062             // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1
0063 
0064             if (pbin < 0)
0065                   continue;
0066             if (tbin < 0)
0067                 continue;
0068             if (phibin < 0)
0069                   continue;
0070                   
0071             if (partPDG == 11) // Look at Thrown Electrons
0072             {
0073                   if (recoPhotons < 3)
0074                   {
0075                         feb[pbin][tbin][phibin]++;
0076                   }
0077                   else
0078                   {
0079                         if (recoPDG == 11)
0080                               fee[pbin][tbin][phibin]++;
0081                         if (recoPDG == 211)
0082                               fepi[pbin][tbin][phibin]++;
0083                         if (recoPDG == 321)
0084                               fek[pbin][tbin][phibin]++;
0085                         if (recoPDG == 2212)
0086                               fep[pbin][tbin][phibin]++;
0087                   }
0088             }
0089 
0090             if (partPDG == 211) // Look at Thrown Pions
0091             {
0092                   if (recoPhotons < 3)
0093                   {
0094                         fpib[pbin][tbin][phibin]++;
0095                   }
0096                   else
0097                   {
0098                         if (recoPDG == 11)
0099                               fpie[pbin][tbin][phibin]++;
0100                         if (recoPDG == 211)
0101                               fpipi[pbin][tbin][phibin]++;
0102                         if (recoPDG == 321)
0103                               fpik[pbin][tbin][phibin]++;
0104                         if (recoPDG == 2212)
0105                               fpip[pbin][tbin][phibin]++;
0106                   }
0107             }
0108 
0109             if (partPDG == 321) // Look at Thrown Kaons
0110             {
0111                   if (recoPhotons < 3)
0112                   {
0113                         fkb[pbin][tbin][phibin]++;
0114                   }
0115                   else
0116                   {
0117                         if (recoPDG == 11)
0118                               fke[pbin][tbin][phibin]++;
0119                         if (recoPDG == 211)
0120                               fkpi[pbin][tbin][phibin]++;
0121                         if (recoPDG == 321)
0122                               fkk[pbin][tbin][phibin]++;
0123                         if (recoPDG == 2212)
0124                               fkp[pbin][tbin][phibin]++;
0125                   }
0126             }
0127 
0128             if (partPDG == 2212) // Look at Thrown Protons
0129             {
0130                   if (recoPhotons < 3)
0131                   {
0132                         fpb[pbin][tbin][phibin]++;
0133                   }
0134                   else
0135                   {
0136                         if (recoPDG == 11)
0137                               fpe[pbin][tbin][phibin]++;
0138                         if (recoPDG == 211)
0139                               fppi[pbin][tbin][phibin]++;
0140                         if (recoPDG == 321)
0141                               fpk[pbin][tbin][phibin]++;
0142                         if (recoPDG == 2212)
0143                               fpp[pbin][tbin][phibin]++;
0144                   }
0145             }
0146       }
0147 
0148       TH2D *h[37][20][120];
0149       for (int pbin = 0; pbin < 37; pbin++)
0150       {
0151             for (int tbin = 0; tbin < 20; tbin++)
0152             {
0153                   for (int phibin = 0; phibin < 120; phibin++)
0154                   {
0155                         h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.);
0156                   }
0157             }
0158       }
0159 
0160       for (int pbin = 0; pbin < 37; pbin++)
0161       {
0162             for (int tbin = 0; tbin < 20; tbin++)
0163             {
0164                   for (int phibin = 0; phibin < 120; phibin++)
0165                   {
0166 
0167                         Double_t fe_tot = fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero;
0168                         fee[pbin][tbin][phibin] /= fe_tot;
0169                         fepi[pbin][tbin][phibin] /= fe_tot;
0170                         fek[pbin][tbin][phibin] /= fe_tot;
0171                         fep[pbin][tbin][phibin] /= fe_tot;
0172                         feb[pbin][tbin][phibin] /= fe_tot;
0173 
0174                         Double_t fpi_tot = fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero;
0175                         fpie[pbin][tbin][phibin] /= fpi_tot;
0176                         fpipi[pbin][tbin][phibin] /= fpi_tot;
0177                         fpik[pbin][tbin][phibin] /= fpi_tot;
0178                         fpip[pbin][tbin][phibin] /= fpi_tot;
0179                         fpib[pbin][tbin][phibin] /= fpi_tot;
0180 
0181                         Double_t fk_tot = fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero;
0182                         fke[pbin][tbin][phibin] /= fk_tot;
0183                         fkpi[pbin][tbin][phibin] /= fk_tot;
0184                         fkk[pbin][tbin][phibin] /= fk_tot;
0185                         fkp[pbin][tbin][phibin] /= fk_tot;
0186                         fkb[pbin][tbin][phibin] /= fk_tot;
0187 
0188                         Double_t fp_tot = fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero;
0189                         fpe[pbin][tbin][phibin] /= fp_tot;
0190                         fppi[pbin][tbin][phibin] /= fp_tot;
0191                         fpk[pbin][tbin][phibin] /= fp_tot;
0192                         fpp[pbin][tbin][phibin] /= fp_tot;
0193                         fpb[pbin][tbin][phibin] /= fp_tot;
0194 
0195                         h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]);
0196                         h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]);
0197                         h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]);
0198                         h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]);
0199                         h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]);
0200 
0201                         h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]);
0202                         h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]);
0203                         h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]);
0204                         h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]);
0205                         h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]);
0206 
0207                         h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]);
0208                         h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]);
0209                         h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]);
0210                         h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]);
0211                         h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]);
0212 
0213                         h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]);
0214                         h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]);
0215                         h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]);
0216                         h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]);
0217                         h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]);
0218 
0219                         Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1);
0220                         Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1);
0221                         Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1);
0222 
0223                         if (pbin == 0 && tbin == 0 && phibin == 0)
0224                         {
0225                               cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl;
0226                               cout << "11   -1 " << p << " " << t << " " << phi << " " << fee[pbin][tbin][phibin] << " " << fepi[pbin][tbin][phibin] << " " << fek[pbin][tbin][phibin] << " " << fep[pbin][tbin][phibin] << " " << feb[pbin][tbin][phibin] << endl;
0227                               cout << "211  1  " << p << " " << t << " " << phi << " " << fpie[pbin][tbin][phibin] << " " << fpipi[pbin][tbin][phibin] << " " << fpik[pbin][tbin][phibin] << " " << fpip[pbin][tbin][phibin] << " " << fpib[pbin][tbin][phibin] << endl;
0228                               cout << "321  1  " << p << " " << t << " " << phi << " " << fke[pbin][tbin][phibin] << " " << fkpi[pbin][tbin][phibin] << " " << fkk[pbin][tbin][phibin] << " " << fkp[pbin][tbin][phibin] << " " << fkb[pbin][tbin][phibin] << endl;
0229                               cout << "2212 1  " << p << " " << t << " " << phi << " " << fpe[pbin][tbin][phibin] << " " << fppi[pbin][tbin][phibin] << " " << fpk[pbin][tbin][phibin] << " " << fpp[pbin][tbin][phibin] << " " << fpb[pbin][tbin][phibin] << endl;
0230                         }
0231                   }
0232             }
0233       }
0234 
0235       FILE *filePointer;
0236       filePointer = fopen(output_txt, "w");
0237 
0238       for (int pbin = 0; pbin < 37; pbin++)
0239       {
0240             for (int tbin = 0; tbin < 20; tbin++)
0241             {
0242                   for (int phibin = 0; phibin < 120; phibin++)
0243                   {
0244                         Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1);
0245                         Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1);
0246                         Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1);
0247 
0248                         fprintf(filePointer, "11   -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]);
0249                   }
0250             }
0251       }
0252 
0253       for (int pbin = 0; pbin < 37; pbin++)
0254       {
0255             for (int tbin = 0; tbin < 20; tbin++)
0256             {
0257                   for (int phibin = 0; phibin < 120; phibin++)
0258                   {
0259                         Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1);
0260                         Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1);
0261                         Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1);
0262 
0263                         fprintf(filePointer, "211  1  %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]);
0264                   }
0265             }
0266       }
0267 
0268       for (int pbin = 0; pbin < 37; pbin++)
0269       {
0270             for (int tbin = 0; tbin < 20; tbin++)
0271             {
0272                   for (int phibin = 0; phibin < 120; phibin++)
0273                   {
0274                         Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1);
0275                         Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1);
0276                         Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1);
0277 
0278                         fprintf(filePointer, "321  1  %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]);
0279                   }
0280             }
0281       }
0282 
0283       for (int pbin = 0; pbin < 37; pbin++)
0284       {
0285             for (int tbin = 0; tbin < 20; tbin++)
0286             {
0287                   for (int phibin = 0; phibin < 120; phibin++)
0288                   {
0289                         Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1);
0290                         Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1);
0291                         Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1);
0292 
0293                         fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]);
0294                   }
0295             }
0296       }
0297 
0298       fclose(filePointer);
0299       ofile->Write();
0300       ofile->Close();
0301 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!