Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-29 10:23:23

0001 void readTreeQA(const char *input, const char *output, const char *output_txt)
0002 {
0003 
0004     TFile *fa = new TFile(input);
0005 
0006     TTree *ta = (TTree *)fa->Get("T");
0007 
0008     TFile *ofile = TFile::Open(output, "recreate");
0009 
0010     Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount;
0011     Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle;
0012     Double_t zero = pow(10., -50.);
0013 
0014     ta->SetBranchAddress("partPDG", &partPDG);
0015     ta->SetBranchAddress("partMom", &partMom);
0016     ta->SetBranchAddress("partEta", &partEta);
0017     ta->SetBranchAddress("partTheta", &partTheta);
0018     ta->SetBranchAddress("partPhi", &partPhi);
0019     ta->SetBranchAddress("partVertX", &partVertX);
0020     ta->SetBranchAddress("partVertY", &partVertY);
0021     ta->SetBranchAddress("partVertZ", &partVertZ);
0022     ta->SetBranchAddress("recoPDG", &recoPDG);
0023     ta->SetBranchAddress("recoPhotons", &recoPhotons);
0024     ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle);
0025 
0026     TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1);
0027     TH1D *h_t = new TH1D("theta", "", 10, 2.7, 3.1);
0028     TH1D *h_phi = new TH1D("phi", "", 12, -1.0 * TMath::Pi(), TMath::Pi());
0029 
0030     TProfile2D *recoPhiVsEtaElecPhoton[4];
0031     TH2D *partPhiVsEtaElec[4];
0032     TH2D *partPhiVsEtaElecMatch[4];
0033     TH2D *partPhiVsEtaElecUnMatch[4];
0034     TH2D *partPhiVsEtaElecMatchPhoton[4];
0035     TH1D *recoElecPhoton[4];
0036 
0037     TProfile2D *recoPhiVsEtaPionPhoton[4];
0038     TH2D *partPhiVsEtaPion[4];
0039     TH2D *partPhiVsEtaPionMatch[4];
0040     TH2D *partPhiVsEtaPionUnMatch[4];
0041     TH1D *recoPionPhoton[4];
0042 
0043     TProfile2D *recoPhiVsEtaKaonPhoton[4];
0044     TH2D *partPhiVsEtaKaon[4];
0045     TH2D *partPhiVsEtaKaonMatch[4];
0046     TH2D *partPhiVsEtaKaonUnMatch[4];
0047     TH1D *recoKaonPhoton[4];
0048 
0049     TProfile2D *recoPhiVsEtaProtonPhoton[4];
0050     TH2D *partPhiVsEtaProton[4];
0051     TH2D *partPhiVsEtaProtonMatch[4];
0052     TH2D *partPhiVsEtaProtonUnMatch[4];
0053     TH1D *recoProtonPhoton[4];
0054 
0055     for (int i = 0; i < 4; i++)
0056     {
0057         recoPhiVsEtaElecPhoton[i] = new TProfile2D(Form("recoPhiVsEtaElecPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0058         partPhiVsEtaElec[i] = new TH2D(Form("partPhiVsEtaElecMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0059         partPhiVsEtaElecMatch[i] = new TH2D(Form("partPhiVsEtaElecMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0060         partPhiVsEtaElecUnMatch[i] = new TH2D(Form("partPhiVsEtaElecUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0061         partPhiVsEtaElecMatchPhoton[i] = new TH2D(Form("partPhiVsEtaElecMatchPhotonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0062         recoElecPhoton[i] = new TH1D(Form("recoElecPhoton_%d", i), "", 25, -0.1, 24.9);
0063 
0064         recoPhiVsEtaPionPhoton[i] = new TProfile2D(Form("recoPhiVsEtaPionPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0065         partPhiVsEtaPion[i] = new TH2D(Form("partPhiVsEtaPionMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0066         partPhiVsEtaPionMatch[i] = new TH2D(Form("partPhiVsEtaPionMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0067         partPhiVsEtaPionUnMatch[i] = new TH2D(Form("partPhiVsEtaPionUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0068         recoPionPhoton[i] = new TH1D(Form("recoPionPhoton_%d", i), "", 25, -0.1, 24.9);
0069 
0070         recoPhiVsEtaKaonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaKaonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0071         partPhiVsEtaKaon[i] = new TH2D(Form("partPhiVsEtaKaonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0072         partPhiVsEtaKaonMatch[i] = new TH2D(Form("partPhiVsEtaKaonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0073         partPhiVsEtaKaonUnMatch[i] = new TH2D(Form("partPhiVsEtaKaonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0074         recoKaonPhoton[i] = new TH1D(Form("recoKaonPhoton_%d", i), "", 25, -0.1, 24.9);
0075 
0076         recoPhiVsEtaProtonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaProtonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0077         partPhiVsEtaProton[i] = new TH2D(Form("partPhiVsEtaProtonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0078         partPhiVsEtaProtonMatch[i] = new TH2D(Form("partPhiVsEtaProtonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0079         partPhiVsEtaProtonUnMatch[i] = new TH2D(Form("partPhiVsEtaProtonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi());
0080         recoProtonPhoton[i] = new TH1D(Form("recoProtonPhoton_%d", i), "", 25, -0.1, 24.9);
0081     }
0082 
0083     TH2D *numPhotonsVsPhiEta25Mom10Pion = new TH2D("numPhotonsVsPhiEta25Mom10Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.);
0084     TH2D *numPhotonsVsPhiEta25Mom4Pion = new TH2D("numPhotonsVsPhiEta25Mom4Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.);
0085 
0086     static Float_t fee[15][10][12], fepi[15][10][12], fek[15][10][12], fep[15][10][12], feb[15][10][12];
0087     static Float_t fpie[15][10][12], fpipi[15][10][12], fpik[15][10][12], fpip[15][10][12], fpib[15][10][12];
0088     static Float_t fke[15][10][12], fkpi[15][10][12], fkk[15][10][12], fkp[15][10][12], fkb[15][10][12];
0089     static Float_t fpe[15][10][12], fppi[15][10][12], fpk[15][10][12], fpp[15][10][12], fpb[15][10][12];
0090 
0091     for (int pbin = 0; pbin < 15; pbin++)
0092     {
0093         for (int tbin = 0; tbin < 10; tbin++)
0094         {
0095             for (int phibin = 0; phibin < 12; phibin++)
0096             {
0097                 fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0;
0098                 fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0;
0099                 fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0;
0100                 fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0;
0101             }
0102         }
0103     }
0104 
0105     Int_t nEntries = (Int_t)ta->GetEntries();
0106     for (int i = 0; i < nEntries; i++)
0107     {
0108         if (i % 100000 == 0)
0109             cout << "Processing Event " << i << endl;
0110         ta->GetEntry(i);
0111 
0112         Int_t pbin = h_p->FindBin(partMom) - 1;
0113         Int_t tbin = h_t->FindBin(partTheta) - 1;
0114         Int_t phibin = h_phi->FindBin(partPhi) - 1;
0115 
0116         Int_t momSelect = -1;
0117         if (partMom > 1.0 && partMom < 3.0)
0118             momSelect = 0;
0119         if (partMom > 3.0 && partMom < 5.0)
0120             momSelect = 1;
0121         if (partMom > 5.0 && partMom < 7.0)
0122             momSelect = 2;
0123         if (partMom > 7.0 && partMom < 9.0)
0124             momSelect = 3;
0125 
0126         // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1
0127 
0128         if (partPDG == 11) // Look at Thrown Electrons
0129         {
0130             if (recoPhotons < 4)
0131             {
0132                 feb[pbin][tbin][phibin]++;
0133             }
0134             else
0135             {
0136                 if (recoPDG == 11)
0137                     fee[pbin][tbin][phibin]++;
0138                 if (recoPDG == 211)
0139                     fepi[pbin][tbin][phibin]++;
0140                 if (recoPDG == 321)
0141                     fek[pbin][tbin][phibin]++;
0142                 if (recoPDG == 2212)
0143                     fep[pbin][tbin][phibin]++;
0144             }
0145             if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2 && partPhi < 2.5)
0146             {
0147                 recoElecPhoton[momSelect]->Fill(recoPhotons);
0148                 recoPhiVsEtaElecPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons);
0149                 partPhiVsEtaElec[momSelect]->Fill(partEta, partPhi);
0150                 if (recoPDG == 11)
0151                     partPhiVsEtaElecMatch[momSelect]->Fill(partEta, partPhi);
0152                 if (recoPDG != 11)
0153                     partPhiVsEtaElecUnMatch[momSelect]->Fill(partEta, partPhi);
0154                 if (recoPDG == 11 && recoPhotons > 0)
0155                     partPhiVsEtaElecMatchPhoton[momSelect]->Fill(partEta, partPhi);
0156             }
0157         }
0158 
0159         if (partPDG == 211) // Look at Thrown Pions
0160         {
0161             if (recoPhotons < 4)
0162             {
0163                 fpib[pbin][tbin][phibin]++;
0164             }
0165             else
0166             {
0167                 if (recoPDG == 11)
0168                     fpie[pbin][tbin][phibin]++;
0169                 if (recoPDG == 211)
0170                     fpipi[pbin][tbin][phibin]++;
0171                 if (recoPDG == 321)
0172                     fpik[pbin][tbin][phibin]++;
0173                 if (recoPDG == 2212)
0174                     fpip[pbin][tbin][phibin]++;
0175             }
0176             if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2 && partPhi < 2.5)
0177             {
0178                 recoPionPhoton[momSelect]->Fill(recoPhotons);
0179                 recoPhiVsEtaPionPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons);
0180                 partPhiVsEtaPion[momSelect]->Fill(partEta, partPhi);
0181                 if (recoPDG == 211)
0182                     partPhiVsEtaPionMatch[momSelect]->Fill(partEta, partPhi);
0183                 if (recoPDG != 211)
0184                     partPhiVsEtaPionUnMatch[momSelect]->Fill(partEta, partPhi);
0185             }
0186             if (partEta < -2.49 && partEta > -2.51)
0187             {
0188                 if (partMom > 9.5 && partMom < 10.5)
0189                 {
0190                     numPhotonsVsPhiEta25Mom10Pion->Fill(partPhi, recoPhotons);
0191                 }
0192                 if (partMom > 3.5 && partMom < 4.5)
0193                 {
0194                     numPhotonsVsPhiEta25Mom4Pion->Fill(partPhi, recoPhotons);
0195                 }
0196             }
0197         }
0198 
0199         if (partPDG == 321) // Look at Thrown Kaons
0200         {
0201             if (recoPhotons < 4)
0202             {
0203                 fkb[pbin][tbin][phibin]++;
0204             }
0205             else
0206             {
0207                 if (recoPDG == 11)
0208                     fke[pbin][tbin][phibin]++;
0209                 if (recoPDG == 211)
0210                     fkpi[pbin][tbin][phibin]++;
0211                 if (recoPDG == 321)
0212                     fkk[pbin][tbin][phibin]++;
0213                 if (recoPDG == 2212)
0214                     fkp[pbin][tbin][phibin]++;
0215             }
0216             if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2 && partPhi < 2.5)
0217             {
0218                 recoKaonPhoton[momSelect]->Fill(recoPhotons);
0219                 recoPhiVsEtaKaonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons);
0220                 partPhiVsEtaKaon[momSelect]->Fill(partEta, partPhi);
0221                 if (recoPDG == 321)
0222                     partPhiVsEtaKaonMatch[momSelect]->Fill(partEta, partPhi);
0223                 if (recoPDG != 321)
0224                     partPhiVsEtaKaonUnMatch[momSelect]->Fill(partEta, partPhi);
0225             }
0226         }
0227 
0228         if (partPDG == 2212) // Look at Thrown Protons
0229         {
0230             if (recoPhotons < 4)
0231             {
0232                 fpb[pbin][tbin][phibin]++;
0233             }
0234             else
0235             {
0236                 if (recoPDG == 11)
0237                     fpe[pbin][tbin][phibin]++;
0238                 if (recoPDG == 211)
0239                     fppi[pbin][tbin][phibin]++;
0240                 if (recoPDG == 321)
0241                     fpk[pbin][tbin][phibin]++;
0242                 if (recoPDG == 2212)
0243                     fpp[pbin][tbin][phibin]++;
0244             }
0245             if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2. && partPhi < 2.5)
0246             {
0247                 recoProtonPhoton[momSelect]->Fill(recoPhotons);
0248                 recoPhiVsEtaProtonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons);
0249                 partPhiVsEtaProton[momSelect]->Fill(partEta, partPhi);
0250                 if (recoPDG == 2212)
0251                     partPhiVsEtaProtonMatch[momSelect]->Fill(partEta, partPhi);
0252                 if (recoPDG != 2212)
0253                     partPhiVsEtaProtonUnMatch[momSelect]->Fill(partEta, partPhi);
0254             }
0255         }
0256     }
0257 
0258     TH2D *h[15][10][12];
0259     for (int pbin = 0; pbin < 15; pbin++)
0260     {
0261         for (int tbin = 0; tbin < 10; tbin++)
0262         {
0263             for (int phibin = 0; phibin < 12; phibin++)
0264             {
0265                 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.);
0266             }
0267         }
0268     }
0269 
0270     for (int pbin = 0; pbin < 15; pbin++)
0271     {
0272         for (int tbin = 0; tbin < 10; tbin++)
0273         {
0274             for (int phibin = 0; phibin < 12; phibin++)
0275             {
0276 
0277                 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;
0278                 fee[pbin][tbin][phibin] /= fe_tot;
0279                 fepi[pbin][tbin][phibin] /= fe_tot;
0280                 fek[pbin][tbin][phibin] /= fe_tot;
0281                 fep[pbin][tbin][phibin] /= fe_tot;
0282                 feb[pbin][tbin][phibin] /= fe_tot;
0283 
0284                 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;
0285                 fpie[pbin][tbin][phibin] /= fpi_tot;
0286                 fpipi[pbin][tbin][phibin] /= fpi_tot;
0287                 fpik[pbin][tbin][phibin] /= fpi_tot;
0288                 fpip[pbin][tbin][phibin] /= fpi_tot;
0289                 fpib[pbin][tbin][phibin] /= fpi_tot;
0290 
0291                 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;
0292                 fke[pbin][tbin][phibin] /= fk_tot;
0293                 fkpi[pbin][tbin][phibin] /= fk_tot;
0294                 fkk[pbin][tbin][phibin] /= fk_tot;
0295                 fkp[pbin][tbin][phibin] /= fk_tot;
0296                 fkb[pbin][tbin][phibin] /= fk_tot;
0297 
0298                 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;
0299                 fpe[pbin][tbin][phibin] /= fp_tot;
0300                 fppi[pbin][tbin][phibin] /= fp_tot;
0301                 fpk[pbin][tbin][phibin] /= fp_tot;
0302                 fpp[pbin][tbin][phibin] /= fp_tot;
0303                 fpb[pbin][tbin][phibin] /= fp_tot;
0304 
0305                 h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]);
0306                 h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]);
0307                 h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]);
0308                 h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]);
0309                 h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]);
0310 
0311                 h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]);
0312                 h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]);
0313                 h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]);
0314                 h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]);
0315                 h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]);
0316 
0317                 h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]);
0318                 h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]);
0319                 h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]);
0320                 h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]);
0321                 h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]);
0322 
0323                 h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]);
0324                 h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]);
0325                 h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]);
0326                 h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]);
0327                 h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]);
0328 
0329                 Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1);
0330                 Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1);
0331                 Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1);
0332 
0333                 if (pbin == 0 && tbin == 0 && phibin == 0)
0334                 {
0335                     cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl;
0336                     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;
0337                     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;
0338                     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;
0339                     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;
0340                 }
0341             }
0342         }
0343     }
0344 
0345     FILE *filePointer;
0346     filePointer = fopen(output_txt, "w");
0347 
0348     for (int pbin = 0; pbin < 15; pbin++)
0349     {
0350         for (int tbin = 0; tbin < 10; tbin++)
0351         {
0352             for (int phibin = 0; phibin < 12; phibin++)
0353             {
0354                 Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1);
0355                 Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1);
0356                 Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1);
0357 
0358                 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]);
0359             }
0360         }
0361     }
0362 
0363     for (int pbin = 0; pbin < 15; pbin++)
0364     {
0365         for (int tbin = 0; tbin < 10; tbin++)
0366         {
0367             for (int phibin = 0; phibin < 12; phibin++)
0368             {
0369                 Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1);
0370                 Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1);
0371                 Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1);
0372 
0373                 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]);
0374             }
0375         }
0376     }
0377 
0378     for (int pbin = 0; pbin < 15; pbin++)
0379     {
0380         for (int tbin = 0; tbin < 10; tbin++)
0381         {
0382             for (int phibin = 0; phibin < 12; phibin++)
0383             {
0384                 Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1);
0385                 Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1);
0386                 Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1);
0387 
0388                 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]);
0389             }
0390         }
0391     }
0392 
0393     for (int pbin = 0; pbin < 15; pbin++)
0394     {
0395         for (int tbin = 0; tbin < 10; tbin++)
0396         {
0397             for (int phibin = 0; phibin < 12; phibin++)
0398             {
0399                 Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1);
0400                 Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1);
0401                 Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1);
0402 
0403                 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]);
0404             }
0405         }
0406     }
0407 
0408     fclose(filePointer);
0409     ofile->Write();
0410     ofile->Close();
0411 }