MIPP(E907): Acceptance.cc Source File

00001 //Standast includes
00002 #include <iostream>
00003 #include <string>
00004 #include <sstream>
00005 #include <cassert>
00006 #include <vector>
00007 #include <fstream>
00008 #include <stdio.h>
00009 #include <iomanip>
00010 #include <cstdlib>
00011 
00012 //ROOT Includes
00013 #include "TRint.h"
00014 #include "TROOT.h"
00015 #include "TFile.h"
00016 #include "TTree.h"
00017 #include "TBranch.h"
00018 #include "TRandom.h"
00019 #include <TH1F.h>
00020 #include <TH2F.h>
00021 #include <TCanvas.h>
00022 #include <TChain.h>
00023 
00024 //MIPP Includes
00025 #include "MIPPEventSummary/MIPPEventSummary.h"
00026 #include "MIPPEventSummary/MIPPMCSummary.h"
00027 #include "DSTAnalysis/DSTUtil.h"
00028 //#include "ScintUtils/ScintPed.h"
00029 #include "NumericalMethods/NMRandom.h"
00030 
00031 using namespace std;
00032 int getPid(MIPPTPCSummary* itpc) {
00033   
00034   int maxLikelyhoodIndex = 0;
00035   
00036   // checking for regular pids
00037   for (int iPart = 1; iPart < PIDDef::kNPID; iPart++){
00038     
00039     if (itpc->LL[iPart] >
00040     itpc->LL[maxLikelyhoodIndex]){
00041       
00042       maxLikelyhoodIndex = iPart;
00043     }
00044     else if (itpc->LL[iPart] == 
00045          itpc->LL[maxLikelyhoodIndex]){
00046       //printf("Two likelyhoods are identical...\n");  
00047     }
00048   }
00049   
00050   return maxLikelyhoodIndex;
00051 }
00052 
00053 int getPidCkov(MIPPDCkovSummary* ickov){
00054   
00055   int maxLikelyhoodIndex = 0;
00056   
00057   // checking for regular pids
00058   for (int iPart = 1; iPart < PIDDef::kNPID; iPart++){
00059     
00060     if (ickov->LL[iPart] >
00061     ickov->LL[maxLikelyhoodIndex]){
00062       
00063       maxLikelyhoodIndex = iPart;
00064     }
00065     else if (ickov->LL[iPart] == 
00066          ickov->LL[maxLikelyhoodIndex]){
00067       
00068       //printf("Two likelyhoods are identical...\n");
00069     }
00070   }
00071   
00072   return maxLikelyhoodIndex;
00073 }
00074 
00075 
00076 int getPidToF(MIPPToFSummary* itof){
00077   
00078   int maxLikelyhoodIndex = 0;
00079   
00080   // checking for regular pids
00081   for (int iPart = 1; iPart < PIDDef::kNPID; iPart++){
00082     
00083     if (itof->LL[iPart] >
00084     itof->LL[maxLikelyhoodIndex]){
00085       
00086       maxLikelyhoodIndex = iPart;
00087     }
00088     else if (itof->LL[iPart] == 
00089          itof->LL[maxLikelyhoodIndex]){
00090       
00091       //printf("Two likelyhoods are identical...\n");
00092     }
00093   }
00094   
00095   return maxLikelyhoodIndex;
00096 }
00097 
00098 int getPidRICH(MIPPRICHSummary* iRICH){
00099   
00100   int maxLikelyhoodIndex = 0;
00101   
00102   // checking for regular pids
00103   for (int iPart = 1; iPart < PIDDef::kNPID; iPart++){
00104     
00105     if (iRICH->LL[iPart] > iRICH->LL[maxLikelyhoodIndex]){
00106       
00107       maxLikelyhoodIndex = iPart;
00108     }
00109     else if (iRICH->LL[iPart] == iRICH->LL[maxLikelyhoodIndex]){
00110       
00111       //printf("Two likelyhoods are identical...\n");  
00112     }
00113   }
00114   
00115   return maxLikelyhoodIndex;
00116 }
00117 
00118 //-----------------------------------------------------------------------------
00119 int main (int argc, char** argv){
00120   if (argc <2){
00121     cout << "Not Enough Arguments!" << endl;
00122     exit(-1);
00123   }
00124   // Extract file names from the command-line arguments
00125   vector<string> list;
00126   DSTUtil::GetFileList(argc,argv,list);
00127   if (list.empty()){
00128     cout << argv[0] << "Found NO files to work with!" << endl;
00129     return -1;
00130   }
00131   //...................................................................
00132   
00133   MIPPEventSummary* fEvt = new MIPPEventSummary;
00134   MIPPMCSummary* fMC = new MIPPMCSummary;
00135   
00136   TFile* fout = TFile::Open("dstAcceptance.root", "recreate");
00137   
00138   //...................................................................
00139   //Constants and arrays
00140   
00141   int const knybin = 6;
00142   int const knxbin = 21;
00143   double binsx[knxbin+1] = {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.4, 1.9, 2.5, 3.2, 4.0, 4.9, 5.9, 7.0, 14.0, 22.0, 31.0, 41.0, 56.0};
00144   double binsy[knybin+1] = {0.7660444431, 0.8660254038, 0.9396926208, 0.9781476007, 0.9902680687, 0.9975640503, 1.0};//Angle Bins
00145   
00146   //---------------------------------------------------------------------------
00147   //Begin Histogram Definitions
00148   TH2F* hReco = new TH2F("ThetaMomReco","Theta vs Momentum for Reco",knxbin,0,knxbin,knybin,0,6.0); 
00149   hReco->Sumw2();
00150   TH2F* hReco2 = new TH2F("ThetaMomReco2","Theta vs Momentum for Reco2",knxbin,0,knxbin,knybin,0,6.0);
00151   hReco2->Sumw2();
00152   TH2F* hRecoPosPi = new TH2F("RecoPosPi","Theta vs Momentum for RecoPosPi",knxbin,0,knxbin,knybin,0,6.0);
00153   hRecoPosPi->Sumw2();
00154   TH2F* hRecoPosK = new TH2F("RecoPosK","Theta vs Momentum for RecoPosK",knxbin,0,knxbin,knybin,0,6.0);
00155   hRecoPosK->Sumw2();
00156   TH2F* hRecoPosP = new TH2F("RecoPosP","Theta vs Momentum for RecoPosP",knxbin,0,knxbin,knybin,0,6.0);
00157   hRecoPosP->Sumw2();
00158   TH2F* hRecoNegPi = new TH2F("RecoNegPi","Theta vs Momentum for RecoNegPi",knxbin,0,knxbin,knybin,0,6.0);
00159   hRecoNegPi->Sumw2();
00160   TH2F* hRecoNegK = new TH2F("RecoNegK","Theta vs Momentum for RecoNegK",knxbin,0,knxbin,knybin,0,6.0);
00161   hRecoNegK->Sumw2();
00162   TH2F* hRecoNegP = new TH2F("RecoNegP","Theta vs Momentum for RecoNegP",knxbin,0,knxbin,knybin,0,6.0);
00163   hRecoNegP->Sumw2();
00164   TH2F* hTruth = new TH2F("Truth","Theta vs Momentum for Truth",knxbin,0,knxbin,knybin,0,6.0); 
00165   hTruth->Sumw2();
00166   TH2F* hTruthPosPi = new TH2F("TruthPosPi","Theta vs Momentum for TruthPosPi",knxbin,0,knxbin,knybin,0,6.0); 
00167   hTruthPosPi->Sumw2();
00168   TH2F* hTruthPosK = new TH2F("TruthPosK","Theta vs Momentum for TruthPosK",knxbin,0,knxbin,knybin,0,6.0); 
00169   hTruthPosK->Sumw2();
00170   TH2F* hTruthPosP = new TH2F("TruthPosP","Theta vs Momentum for TruthPosP",knxbin,0,knxbin,knybin,0,6.0); 
00171   hTruthPosP->Sumw2();
00172   TH2F* hTruthNegPi = new TH2F("TruthNegPi","Theta vs Momentum for TruthNegPi",knxbin,0,knxbin,knybin,0,6.0); 
00173   hTruthNegPi->Sumw2();
00174   TH2F* hTruthNegK = new TH2F("TruthNegK","Theta vs Momentum for TruthNegK",knxbin,0,knxbin,knybin,0,6.0); 
00175   hTruthNegK->Sumw2();
00176   TH2F* hTruthNegP = new TH2F("TruthNegP","Theta vs Momentum for TruthNegP",knxbin,0,knxbin,knybin,0,6.0); 
00177   hTruthNegP->Sumw2();
00178   TH2F* hAcceptance = new TH2F("Accpt","Detector Acceptance",knxbin,0,knxbin,knybin,0,6.0);
00179   hAcceptance->Sumw2();
00180   TH2F* hAcceptPPi = new TH2F("AccptPPi","Detector AcceptancePPi",knxbin,0,knxbin,knybin,0,6.0);
00181   hAcceptPPi->Sumw2();
00182   TH2F* hAcceptPK = new TH2F("AccptPK","Detector AcceptancePK",knxbin,0,knxbin,knybin,0,6.0);
00183   hAcceptPK->Sumw2();
00184   TH2F* hAcceptPP = new TH2F("AccptPP","Detector AcceptancePP",knxbin,0,knxbin,knybin,0,6.0);
00185   hAcceptPP->Sumw2();
00186   TH2F* hAcceptNPi = new TH2F("AccptNPi","Detector AcceptanceNPi",knxbin,0,knxbin,knybin,0,6.0);
00187   hAcceptNPi->Sumw2();
00188   TH2F* hAcceptNK = new TH2F("AccptNK","Detector AcceptanceNK",knxbin,0,knxbin,knybin,0,6.0);
00189   hAcceptNK->Sumw2();
00190   TH2F* hAcceptNP = new TH2F("AccptNP","Detector AcceptanceNP",knxbin,0,knxbin,knybin,0,6.0);
00191   hAcceptNP->Sumw2();
00192   TH2F* hTrackPi = new TH2F("TrackPi","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00193   hTrackPi->Sumw2();
00194   TH2F* hTrackK = new TH2F("TrackK","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00195   hTrackK->Sumw2();
00196   TH2F* hTrackP = new TH2F("TrackP","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00197   hTrackP->Sumw2();
00198   TH2F* hTrackAPi = new TH2F("TrackAPi","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00199   hTrackAPi->Sumw2();
00200   TH2F* hTrackAK = new TH2F("TrackAK","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00201   hTrackAK->Sumw2();
00202   TH2F* hTrackAP = new TH2F("TrackAP","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00203   hTrackAP->Sumw2();
00204   TH2F* hTrackEffPi = new TH2F("TrackEffPi","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00205   hTrackEffPi->Sumw2();
00206   TH2F* hTrackEffK = new TH2F("TrackEffK","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00207   hTrackEffK->Sumw2();
00208   TH2F* hTrackEffP = new TH2F("TrackEffP","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00209   hTrackEffP->Sumw2();
00210   TH2F* hTrackEffAPi = new TH2F("TrackEffAPi","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00211   hTrackEffAPi->Sumw2();
00212   TH2F* hTrackEffAK = new TH2F("TrackEffAK","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00213   hTrackEffAK->Sumw2();
00214   TH2F* hTrackEffAP = new TH2F("TrackEffAP","Track Matched",knxbin,0,knxbin,knybin,0,6.0);
00215   hTrackEffAP->Sumw2();
00216 
00217   
00218   //  TH1F* hHistoCut = new TH1F("cuts", "Number of Events after Cuts", 5, -0.5, 4.5);
00219   TH1F* hHistoPZ = new TH1F("pZ", "Momentum Z",100, -100, 60);
00220   TH1F* hHistoPT = new TH1F("pT", "Transverse Momentum", 100, -0.5,1.5);
00221   TH1F* hHistoPtot = new TH1F("Ptot", "Momentum",500, -100, 100);
00222   TH1F* hHistotht = new TH1F("tht", "cos tht", 500, -1.5,1.5);
00223   TH1F* hHistoTgt = new TH1F("tgt", "Target Position", 10000,-835,-825);
00224   TH2F* hDedxMomPi = new TH2F("DedxMomPi","Energy Loss vs Momentum",100, 0, 1.2, 500, 0, 50);
00225   TH2F* hDedxMomK = new TH2F("DedxMomK","Energy Loss vs Momentum",100, 0, 1.2, 500, 0, 50);
00226   TH2F* hDedxMomP = new TH2F("DedxMomP","Energy Loss vs Momentum",100, 0, 1.2, 500, 0, 50);
00227   TH1F* hRecoPID = new TH1F("recopid","Reconstructed Particles", 8, -0.5, 7.5);
00228   TH1F* hRecoPIDAll = new TH1F("recopidall","Reconstructed Particles", 8, -0.5, 7.5);
00229   TH1F* hRecoPID1 = new TH1F("recopid1","Reconstructed Particles", 8, -0.5, 7.5);
00230   TH1F* hTruthPID = new TH1F("truthpid","Truth Particles", 8, -0.5, 7.5);
00231   TH1F* hTruthPID1 = new TH1F("truthpid1","Truth Particles", 8, -0.5, 7.5);
00232   //End  histogram definitions
00233   //--------------------------------------------------------------------------
00234   int binPtot, binTeta;// binPtotR, binTetaR;
00235   //------------------------------------------------------------------------- 
00236   // Loop over all files that were inserted into the file name list
00237   for (unsigned int ifi = 0; ifi <list.size() && !gStop; ++ifi) {
00238     TFile* f = TFile::Open(list[ifi].c_str());
00239     if (!f) {
00240       cout << "Couldn't open " << list[ifi] << endl;
00241       continue;
00242     }
00243     // Get event tree from the file
00244     TTree* evtTree = DSTUtil::GetEvtTree(f);
00245     if (!evtTree) {
00246       cout << "Couldn't find event tree in " << list[ifi] << endl;
00247       continue;
00248     }
00249     int nEnt = evtTree->GetEntries();
00250     
00251     cout << "DST " << ifi << "/" << list.size()
00252          << " with " << nEnt << " entries" << endl;
00253     
00254     evtTree->SetBranchAddress("fvtxconfit.", &fEvt);
00255     evtTree->SetBranchAddress("fMCTruth.", &fMC);
00256     fout->cd();
00257     
00258     //------------------------------------------------------------------------
00259     //Beginning of the for loop to process track counting both truth and reco.
00260     for (int ent = 0; ent < nEnt; ++ent){
00261       evtTree->GetEntry(ent);
00262       if(fEvt->NBTrk() != 1) continue;
00263       for(int i = 0; i < fMC->NTrk(); ++i){ 
00264     MIPPMCTrackSummary* mctrk = fMC->GetTrk(i);
00265     MIPPMCVertexSummary* v=fMC->GetVtx(mctrk->vtxindex);
00266     hHistoTgt->Fill(v->x[2]);
00267     if(!(mctrk->pid == PIDDef::kPion || mctrk->pid == PIDDef::kKaon || mctrk->pid == PIDDef::kProton))continue; // Select only pion kaon and proton
00268     if(v->x[2]<-829 || v->x[2]>-822) continue;// Target position cut
00269     
00270     float pX = mctrk->p[0];
00271     float pY = mctrk->p[1];
00272     float pZ = mctrk->p[2];
00273     hHistoPZ->Fill(pZ);
00274     float pT = sqrt(pX*pX+pY*pY);
00275     hHistoPT->Fill(pT);
00276     double Ptot = sqrt(pX*pX + pY*pY + pZ*pZ);
00277     
00278     if(pZ == 0) continue;// Cut for avoiding indefinite value for pT/fabs(pZ)
00279     double costeta = cos(atan(pT/fabs(pZ)));
00280     
00281     hHistoPtot->Fill(Ptot);
00282     hHistotht->Fill(costeta);
00283     
00284     // Get track ptot and find the appropriate bin
00285     binPtot = -1;
00286     for (int j = -1; j < knxbin; ++j){ 
00287       if (Ptot < binsx[j+1]) {
00288         binPtot = j; 
00289         break;
00290       }
00291     }
00292     // Get track cos(teta) and find the appropriate bin
00293     binTeta = -1;
00294     for (int j = -1; j < knybin; ++j){ 
00295       if (costeta < binsy[j+1]){
00296         binTeta = j; 
00297         break;
00298       }
00299     }
00300     hTruth->Fill(binPtot,binTeta);//Fills truth plot
00301     if(mctrk->q > 0) {
00302       if(mctrk->pid == 2) hTruthPosPi->Fill(binPtot,binTeta);
00303       if(mctrk->pid==3) hTruthPosK->Fill(binPtot,binTeta);
00304       if(mctrk->pid==4) hTruthPosP->Fill(binPtot,binTeta);
00305     }
00306     if(mctrk->q < 0){
00307       if(mctrk->pid == 2) hTruthNegPi->Fill(binPtot,binTeta);
00308       if(mctrk->pid==3) hTruthNegK->Fill(binPtot,binTeta);
00309       if(mctrk->pid==4) hTruthNegP->Fill(binPtot,binTeta);
00310     }
00311     
00312     //Now look at reco tracks
00313     int rcindex = mctrk->recoindex;
00314     int pidTPC = -1;
00315     
00316     if(rcindex>-1){
00317       MIPPTrackSummary* trk = fEvt->GetTrk(rcindex);
00318       MIPPRICHSummary* rich = fEvt->GetRICH(trk->irich);
00319       MIPPTPCSummary* tpc = fEvt->GetTPC(trk->itpc);
00320       MIPPToFSummary* tof = fEvt->GetToF(trk->itof[0]);
00321       MIPPDCkovSummary* ckov = fEvt->GetDCkov(trk->idckov);
00322       if(mctrk->pid ==2)hDedxMomPi->Fill(Ptot, tpc->dedx);
00323       if(mctrk->pid ==3)hDedxMomK->Fill(Ptot, tpc->dedx);
00324       if(mctrk->pid ==4)hDedxMomP->Fill(Ptot, tpc->dedx);
00325       double PTot = trk->ptot;
00326       
00327       if(PTot<=1) {
00328         if (tpc) pidTPC = getPid(tpc);
00329       }
00330       else if( PTot<=2.6 && PTot>1.0) {
00331         if(tof)  pidTPC = getPidToF(tof);
00332       }
00333       else if(Ptot<= 17.0 && Ptot>2.6) {
00334         if(ckov) pidTPC = getPidCkov(ckov);
00335       }
00336       else {
00337         if(rich) pidTPC = getPidRICH(rich);
00338       }
00339       if(mctrk->q > 0){
00340         if(mctrk->pid == 2 && (pidTPC == 2 || pidTPC == 1)) hRecoPosPi->Fill(binPtot,binTeta);
00341         if(mctrk->pid == 3 && pidTPC == 3) hRecoPosK->Fill(binPtot,binTeta);
00342         if(mctrk->pid == 4 && pidTPC == 4 ) hRecoPosP->Fill(binPtot,binTeta);
00343         if(mctrk->pid == 2) hTrackPi->Fill(binPtot,binTeta);
00344         if(mctrk->pid == 3) hTrackK->Fill(binPtot,binTeta);
00345         if(mctrk->pid == 4) hTrackP->Fill(binPtot,binTeta);
00346       }
00347       if(mctrk->q < 0){
00348         if(mctrk->pid == 2 && (pidTPC == 2 || pidTPC == 1)) hRecoNegPi->Fill(binPtot,binTeta);
00349         if(mctrk->pid == 3 && pidTPC == 3) hRecoNegK->Fill(binPtot,binTeta);
00350         if(mctrk->pid == 4 && pidTPC == 4) hRecoNegP->Fill(binPtot,binTeta);
00351         if(mctrk->pid == 2) hTrackAPi->Fill(binPtot,binTeta);
00352         if(mctrk->pid == 3) hTrackAK->Fill(binPtot,binTeta);
00353         if(mctrk->pid == 4) hTrackAP->Fill(binPtot,binTeta);
00354       }
00355       hReco->Fill(binPtot,binTeta);//Fills reco plot
00356       if(mctrk->q > 0){
00357         if(binPtot ==0 && mctrk->pid ==4) hRecoPID->Fill(pidTPC);
00358         if(binPtot ==1 && mctrk->pid ==4) hRecoPID1->Fill(pidTPC);
00359         if(binPtot ==0 && mctrk->pid ==4) hTruthPID->Fill(mctrk->pid);
00360         if(binPtot ==1 && mctrk->pid ==4) hTruthPID1->Fill(mctrk->pid);
00361         if(mctrk->pid ==4) hRecoPIDAll->Fill(pidTPC);
00362         
00363       }
00364     }
00365 
00366       }
00367       //End of the for loop to process track counting for both truth and reco.   
00368       hTrackEffPi->Divide(hTrackPi,hTruthPosPi,1,1);
00369       hTrackEffK->Divide(hTrackK,hTruthPosK,1,1);
00370       hTrackEffP->Divide(hTrackP,hTruthPosP,1,1);
00371       hTrackEffAPi->Divide(hTrackAPi,hTruthNegPi,1,1);
00372       hTrackEffAK->Divide(hTrackAK,hTruthNegK,1,1);
00373       hTrackEffAP->Divide(hTrackAP,hTruthNegP,1,1);
00374       hAcceptance->Divide(hReco,hTruth,1,1);//Fills acceptance plot
00375       hAcceptPPi->Divide(hRecoPosPi,hTruthPosPi,1,1);
00376       hAcceptPK->Divide(hRecoPosK,hTruthPosK,1,1);
00377       hAcceptPP->Divide(hRecoPosP,hTruthPosP,1,1);
00378       hAcceptNPi->Divide(hRecoNegPi,hTruthNegPi,1,1);
00379       hAcceptNK->Divide(hRecoNegK,hTruthNegK,1,1);
00380       hAcceptNP->Divide(hRecoNegP,hTruthNegP,1,1);
00381       
00382     }
00383     //End of the entries loop
00384     f->Close(); 
00385   }
00386   //End of the all files loop
00387   fout->Write();
00388   delete fout;
00389   return 0;
00390 }
00391 //End of the Main Function