benchmarks/lfhcal/LFHCAL_Performance.C

0001 // Code to extract the Tracking Performances
0002 // Shyam Kumar; INFN Bari, Italy
0003 // shyam.kumar@ba.infn.it; shyam.kumar@cern.ch
0004
0005 #include "TGraphErrors.h"
0006 #include "TF1.h"
0007 #include "TRandom.h"
0008 #include "TCanvas.h"
0009 #include "TLegend.h"
0010 #include "TMath.h"
0011 #include "TVector3.h"
0012
0013 #define mpi 0.139  // 1.864 GeV/c^2
0014
0015 void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.0, TString name = "")
0016 {
0017
0018   // style of the plot
0019    gStyle->SetPalette(1);
0020    gStyle->SetOptTitle(1);
0021    gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y");
0022    gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y");
0023    gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
0024    gStyle->SetHistLineWidth(2);
0025    gStyle->SetOptFit(1);
0026    gStyle->SetOptStat(1);
0027
0028    TString dir = "";
0029    TString dist_dir_mom = "mom_resol";
0030
0031    bool debug=true;
0032   // Tree with reconstructed tracks
0033    const int nbins_eta = 5;
0034    int nfiles = 100;
0035    double eta[nbins_eta+1]={1.2,1.5,2,2.5,3,3.5};
0036    double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1};
0037    TH1D *histp[nbins_eta];
0038
0039
0040    for (int i=0; i<nbins_eta; i++){
0041    histp[i] = new TH1D(Form("hist_etabin%d",i),Form("hist_etabin%d",i),600,-1,1);
0042    histp[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom));
0043    histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0044    }
0045
0046    TFile* file = TFile::Open(filename.Data());
0047    if (!file) {printf("file not found !!!"); return;}
0048    TTreeReader myReader("events", file); // name of tree and file
0049    if (debug) cout<<"Filename: "<<file->GetName()<<"\t NEvents: "<<myReader.GetEntries()<<endl;
0050
0051    // MC and Reco information
0052    TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge");
0053    TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x");
0054    TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y");
0055    TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z");
0056    TTreeReaderArray<Float_t> px_mc(myReader, "MCParticles.momentum.x");
0057    TTreeReaderArray<Float_t> py_mc(myReader, "MCParticles.momentum.y");
0058    TTreeReaderArray<Float_t> pz_mc(myReader, "MCParticles.momentum.z");
0059    TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus");
0060    TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG");
0061
0062    TTreeReaderArray<Float_t> pe_lc(myReader, "LFHCALClusters.energy");
0063    TTreeReaderArray<Float_t> px_lc(myReader, "LFHCALClusters.position.x");
0064    TTreeReaderArray<Float_t> py_lc(myReader, "LFHCALClusters.position.y");
0065    TTreeReaderArray<Float_t> pz_lc(myReader, "LFHCALClusters.position.z");
0066    TTreeReaderArray<Float_t> pe_ec(myReader, "EcalEndcapPClusters.energy");
0067    TTreeReaderArray<Float_t> px_ec(myReader, "EcalEndcapPClusters.position.x");
0068    TTreeReaderArray<Float_t> py_ec(myReader, "EcalEndcapPClusters.position.y");
0069    TTreeReaderArray<Float_t> pz_ec(myReader, "EcalEndcapPClusters.position.z");
0070
0071   int count =0;
0072   int matchId = 1; // Always matched track assigned the index 0
0073   while (myReader.Next())
0074     {
0075       std::cout << "events = " << count++ << std::endl;
0076       for (int j = 0; j < pdg.GetSize(); ++j)
0077     {
0078       if (status[j] !=1 && pdg.GetSize()!=1) continue;
0079       Double_t pzmc = pz_mc[j];
0080       Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]);
0081       Double_t pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec
0082       Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2));
0083       Double_t phimc = TMath::ATan2(py_mc[j],px_mc[j]);
0084       std::cout << "neutron p=" << pmc << " pt=" << ptmc << std::endl;
0085
0086       if (fabs(ptmc) < pTcut) continue;
0087
0088       float l_px_tot=0;
0089       float l_py_tot=0;
0090       float l_pz_tot=0;
0091       float l_e_tot=0;
0092
0093       std::cout << "LFHCAL nclus=" << px_lc.GetSize() << " ECAL nclus=" << pe_ec.GetSize() << std::endl;
0094       for (int jl = 0;jl<px_lc.GetSize();jl++)
0095         {
0096           float e = pe_lc[jl];
0097           TVector3 v(px_lc[jl],py_lc[jl],pz_lc[jl]);
0098           v.Print();
0099           float eta = v.PseudoRapidity();
0100           float phi = v.Phi();
0101           float pt = e/cosh(eta);
0102           std::cout << "LFHCAL clus: e=" << e << " eta=" << eta << " pt=" << pt << std::endl;
0103           l_e_tot += e;
0104           l_px_tot += pt*cos(phi);
0105           l_py_tot += pt*sin(phi);
0106           l_pz_tot += pt*sinh(eta);
0107         }
0108
0109       float e_px_tot=0;
0110       float e_py_tot=0;
0111       float e_pz_tot=0;
0112       float e_e_tot=0;
0113
0114       for (int je = 0;je<px_ec.GetSize();je++)
0115         {
0116           float e = pe_ec[je];
0117           TVector3 v(px_ec[je],py_ec[je],pz_ec[je]);
0118           float eta = v.PseudoRapidity();
0119           float phi = v.Phi();
0120           float pt = e/cosh(eta);
0121           std::cout << "ECAL clus: e=" << e << " eta=" << eta << " pt=" << pt << std::endl;
0122           e_e_tot += e;
0123           e_px_tot += pt*cos(phi);
0124           e_py_tot += pt*sin(phi);
0125           e_pz_tot += pt*sinh(eta);
0126         }
0127
0128       std::cout << "LFHCAL e=" <<l_e_tot << " ECAL e=" << e_e_tot << std::endl;
0129       float px_tot = l_px_tot+e_px_tot;
0130       float py_tot = l_py_tot+e_py_tot;
0131       float pz_tot = l_pz_tot+e_pz_tot;
0132       float e_tot = l_e_tot+e_e_tot;
0133
0134       float prec = sqrt(px_tot*px_tot+py_tot*py_tot+pz_tot*pz_tot);
0135       float ptrec = sqrt(px_tot*px_tot+py_tot*py_tot);
0136       float pzrec = pz_tot;
0137
0138       float p_resol = (e_tot-pmc)/pmc;
0139       std::cout << "p_resol = " << p_resol << std::endl;
0140       for (int ibin=0; ibin<nbins_eta; ++ibin){
0141         if(etamc>eta[ibin] && etamc<eta[ibin+1]) histp[ibin]->Fill(p_resol);
0142       }
0143     } // Generated Tracks
0144
0145     }// event loop ends
0146
0147    TFile *fout_mom = new TFile(Form("%s/mom/lfhcal_mom_%1.1f_%s_%s.root",particle.Data(),mom,dist_dir_mom.Data(),particle.Data()),"recreate");
0148    fout_mom->cd();
0149    for (int ibin=0; ibin<nbins_eta; ++ibin) histp[ibin]->Write();
0150    fout_mom->Close();
0151
0152 }
0153
0154
0155
0156