JetsAndHF/Job_submission_RCF/analysis.C

0001 #ifdef __CINT__
0002
0003 #pragma link off all globals;
0004 #pragma link off all classes;
0005 #pragma link off all functions;
0006
0007 #pragma link C++ class PlotFile;
0008 #endif
0009
0010 #ifndef __CINT__
0011 #include <stdio.h>
0012 #include <stdlib.h>
0013 #include <fstream>
0014 #include <iostream>
0015 #include <iomanip>
0016 #include <string>
0017 #include <sys/types.h>
0018 #include <sys/stat.h>
0019 #include <dirent.h>
0020 #include "math.h"
0021 #include "string.h"
0022
0023 #include "TROOT.h"
0024 #include "TFile.h"
0025 #include "TChain.h"
0026 #include "TH1F.h"
0027 #include "TH1D.h"
0028 #include "TH2D.h"
0029 #include "TH3D.h"
0030 #include "THnSparse.h"
0031 #include "TStyle.h"
0032 #include "TCanvas.h"
0033 #include "TProfile.h"
0034 #include "TTree.h"
0035 #include "TNtuple.h"
0036 #include "TRandom3.h"
0037 #include "TMath.h"
0038 #include "TSystem.h"
0039 #include "TUnixSystem.h"
0040 #include "TVector2.h"
0041 #include "TVector3.h"
0042 #include "TLorentzVector.h"
0043 #include "TTreeReader.h"
0044 #include "TTreeReaderValue.h"
0045 #include "TTreeReaderArray.h"
0046 #include "TLatex.h"
0047 #endif
0048
0049 using namespace std;
0050
0051 //==================
0052 void analysis(TString listname="", TString outname="")
0053 {
0054   if(listname=="" || outname=="")
0055     {
0056       printf("[e] Missing input list name or output file name.\n");
0057       return;
0058     }
0059
0060   printf("[i] Input list: %s\n", listname.Data());
0061   printf("[i] Output file: %s\n", outname.Data());
0062
0063   TChain *chain = new TChain("events");
0064
0065   int nfiles = 0;
0066   char filename[512];
0067   ifstream *inputstream = new ifstream;
0068   inputstream->open(listname.Data());
0069   if(!inputstream)
0070     {
0071       printf("[e] Cannot open file list: %s\n", listname.Data());
0072     }
0073   while(inputstream->good())
0074     {
0075       inputstream->getline(filename, 512);
0076       if(inputstream->good())
0077     {
0078       TFile *ftmp = TFile::Open(filename, "read");
0079       if(!ftmp||!(ftmp->IsOpen())||!(ftmp->GetNkeys()))
0080         {
0081           printf("[e] Could you open file: %s\n", filename);
0082         }
0083       else
0084         {
0085           cout<<"[i] Add "<<nfiles<<"th file: "<<filename<<endl;
0086           chain->Add(filename);
0087           nfiles++;
0088         }
0089     }
0090     }
0091   inputstream->close();
0092   printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0093
0094   TH1F *hEventStat = new TH1F("hEventStat", "Event statistics", 8, 0, 8);
0095   hEventStat->GetXaxis()->SetBinLabel(1, "All events");
0096
0097   TH1F *hMcMult = new TH1F("hMcMult", "MC multiplicity (|#eta| < 3.5);N_{MC}", 50, 0, 50);
0098   TH1F *hMcVtxX = new TH1F("hMcVtxX", "x position of MC vertex;x (mm)", 100, -5.05, 4.95);
0099   TH1F *hMcVtxY = new TH1F("hMcVtxY", "y position of MC vertex;y (mm)", 500, -5.01, 4.99);
0100   TH1F *hMcVtxZ = new TH1F("hMcVtxZ", "z position of MC vertex;z (mm)", 400, -200, 200);
0101
0102   TH1F *hNRecoJet = new TH1F("hNRecoJet", "Number of reconstructed jets with |#eta| < 2.5 and p_{T} > 2;N", 10, 0, 10);
0103   TH1F *hRecoJetPt = new TH1F("hRecoJetPt", "Reconstructed jets p_{T};p_{T} (GeV/c)", 50, 0, 50);
0104
0105   TTreeReader treereader(chain);
0106
0107   // MC
0108   TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0109   TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0110   TTreeReaderArray<float> mcPartCharge = {treereader, "MCParticles.charge"};
0111   TTreeReaderArray<unsigned int> mcPartParent_begin = {treereader, "MCParticles.parents_begin"};
0112   TTreeReaderArray<unsigned int> mcPartParent_end = {treereader, "MCParticles.parents_end"};
0113   TTreeReaderArray<int> mcPartParent_index = {treereader, "_MCParticles_parents.index"};
0114   TTreeReaderArray<unsigned int> mcPartDaughter_begin = {treereader, "MCParticles.daughters_begin"};
0115   TTreeReaderArray<unsigned int> mcPartDaughter_end = {treereader, "MCParticles.daughters_end"};
0116   TTreeReaderArray<int> mcPartDaughter_index = {treereader, "_MCParticles_daughters.index"};
0117   TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0118   TTreeReaderArray<double> mcPartVx = {treereader, "MCParticles.vertex.x"};
0119   TTreeReaderArray<double> mcPartVy = {treereader, "MCParticles.vertex.y"};
0120   TTreeReaderArray<double> mcPartVz = {treereader, "MCParticles.vertex.z"};
0121   TTreeReaderArray<double> mcMomPx = {treereader, "MCParticles.momentum.x"};
0122   TTreeReaderArray<double> mcMomPy = {treereader, "MCParticles.momentum.y"};
0123   TTreeReaderArray<double> mcMomPz = {treereader, "MCParticles.momentum.z"};
0124   TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0125   TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0126   TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0127
0128   // reconstructed jets
0129   TTreeReaderArray<int> recoJetType = {treereader, "ReconstructedChargedJets.type"};
0130   TTreeReaderArray<float> recoJetE = {treereader, "ReconstructedChargedJets.energy"};
0131   TTreeReaderArray<float> recoJetMomPx = {treereader, "ReconstructedChargedJets.momentum.x"};
0132   TTreeReaderArray<float> recoJetMomPy = {treereader, "ReconstructedChargedJets.momentum.y"};
0133   TTreeReaderArray<float> recoJetMomPz = {treereader, "ReconstructedChargedJets.momentum.z"};
0134   TTreeReaderArray<float> recoJetMass = {treereader, "ReconstructedChargedJets.mass"};
0135
0136   int nevents = 0;
0137   while(treereader.Next())
0138     {
0139       nevents++;
0140       if(nevents%1000==0) printf("\n[i] New event %d\n",nevents);
0141       hEventStat->Fill(0.5);
0142
0143
0144       //===== find MC primary vertex
0145       //===== find 4-momentum of incoming electron and p/A
0146       int nMCPart = mcPartMass.GetSize();
0147       TVector3 vertex_mc(-999., -999., -999.);
0148       for(int imc=0; imc<nMCPart; imc++)
0149     {
0150       int status = mcPartGenStatus[imc];
0151       int pdg = mcPartPdg[imc];
0152
0153       if(status == 4)
0154         {
0155           double mcE = sqrt(pow(mcMomPx[imc],2)+pow(mcMomPy[imc],2)+pow(mcMomPz[imc],2)+pow(mcPartMass[imc],2));
0156           if(pdg == 11)
0157         {
0158           vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0159         }
0160         }
0161     }
0162       hMcVtxX->Fill(vertex_mc.x());
0163       hMcVtxY->Fill(vertex_mc.y());
0164       hMcVtxZ->Fill(vertex_mc.z());
0165
0166       //===== Loop over MC primary particles
0167       int nMcPart = 0;
0168       for(int imc=0; imc<nMCPart; imc++)
0169     {
0170       if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0171         {
0172           double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));
0173           if(dist < 1e-4)
0174         {
0175           // count charged particles within |eta| < 3.5
0176           TVector3 mc_mom(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);
0177           double mcEta = mc_mom.PseudoRapidity();
0178           if(fabs(mcEta) < 3.5) nMcPart++;
0179         }
0180         }
0181     }
0182       hMcMult->Fill(nMcPart);
0183
0184       //===== select di-jet events
0185       int nRecoJets = recoJetType.GetSize();
0186       int nGoodJet = 0;
0187       for (int i=0; i<nRecoJets; i++)
0188     {
0189       // select jets within acceptance
0190       TVector3 jetMom(recoJetMomPx[i],recoJetMomPy[i],recoJetMomPz[i]);
0191       if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0192       hRecoJetPt->Fill(jetMom.Perp());
0193       if(jetMom.Perp()>2) nGoodJet++;
0194     }
0195       hNRecoJet->Fill(nGoodJet);
0196     }
0197
0198   TFile *outfile = new TFile(outname.Data(), "recreate");
0199
0200   hEventStat->Write();
0201   hMcMult->Write();
0202   hMcVtxX->Write();
0203   hMcVtxY->Write();
0204   hMcVtxZ->Write();
0205
0206   hNRecoJet->Write();
0207   hRecoJetPt->Write();
0208
0209   outfile->Close();
0210
0211   return;
0212 }