Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:43:45 

0001 #define TreeAnalysisHcalScale_cxx
0002 #include "TreeAnalysisHcalScale.h"
0003 #include <TStyle.h>
0004 #include <TCanvas.h>
0005 
0006 Bool_t FillChain(TChain *chain, const char *inputFileList);
0007 
0008 int main(Int_t argc, Char_t *argv[]) {
0009   if( argc<4 ){
0010     std::cerr << "Please give 4 arguments \n"
0011               << "runList SeedName" << "\n" 
0012           << "outputFileName" << "\n" 
0013           << "maximum sample size" << "\n"
0014           << "Which sample to do (-1 means all)" << "\n"
0015               << std::endl;
0016     return -1;
0017   }
0018   
0019   const char *inputFileList = argv[1];
0020   const char *outFileName   = argv[2];
0021   double totalTracks        = atof(argv[3]);
0022   int    sample             = atoi(argv[4]);
0023   if (sample < 0 || sample > 5) sample = -1;
0024 
0025   int    cut = 0;
0026 
0027   std::cout << "runList SeeName        " << inputFileList << "\n"
0028         << "outputFileName         " << outFileName << "\n"
0029         << "total number of tracks " << totalTracks 
0030         << " Sample "                << sample << std::endl;
0031 
0032   // Reading Tree                                                        
0033   std::cout << "---------------------" << std::endl;
0034   std::cout << "Reading List of input trees from " << inputFileList << std::endl;
0035   std::string particles[6] = {"pi+","pi-","K+", "K-", "p+", "p-"};
0036   double evFrac[6]         = {0.695, 0.695, 0.163, 0.163, 0.142, 0.142};
0037   std::vector<std::string> particleNames;
0038   std::vector<double>      fraction;
0039   if (sample < 0) {
0040     for (unsigned int i=0; i<6; i++) {
0041       particleNames.push_back(particles[i]);
0042       fraction.push_back(evFrac[i]);
0043     }
0044   } else {
0045     particleNames.push_back(particles[sample]);
0046     fraction.push_back(1.0);
0047   }
0048   
0049   TreeAnalysisHcalScale tree(outFileName, particleNames);
0050 
0051   for (unsigned int i=0; i<particleNames.size(); i++) {
0052     char fileList[200], treeName[200];
0053     sprintf (fileList, "%s%s.txt", inputFileList, particles[i].c_str());
0054     //    sprintf (treeName, "/IsolatedTracksHcalScale/tree");
0055     sprintf (treeName, "/isolatedTracksHcal/tree");
0056     TChain *chain = new TChain(treeName);
0057     std::cout << "try to create a chain for " << fileList << std::endl;
0058     if  (! FillChain(chain, fileList) ) {
0059       std::cerr << "Cannot get the tree " << std::endl;
0060     } else {
0061       unsigned int nmax = (unsigned int)(fraction[i]*totalTracks);
0062       tree.Init(chain);
0063       tree.setParticle(i, nmax);
0064       tree.Loop(cut);
0065       tree.weights[i]= (double)nmax/(double)tree.nIsoTrkTotal;
0066       std::cout << "particle " << i << " nmax " << nmax << " nisotrktotal  " << tree.nIsoTrkTotal << " weight " << tree.weights[i] << std::endl;
0067       tree.clear();
0068       
0069     }
0070   }
0071   tree.AddWeight(particleNames);
0072   
0073   return 0;
0074 }
0075 
0076 Bool_t FillChain(TChain *chain, const char *inputFileList) {
0077 
0078   ifstream infile(inputFileList);
0079   std::string buffer;
0080 
0081   if(!infile.is_open()) {
0082     std::cerr << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
0083     return kFALSE;
0084   }
0085   
0086   //  std::cout << "TreeUtilities : FillChain " << std::endl;
0087   while(1) {
0088     infile >> buffer;
0089     if(!infile.good()) break;
0090 
0091     //    std::cout << "Adding tree from " << buffer.c_str() << std::endl;
0092     chain->Add(buffer.c_str());
0093   }
0094   std::cout << "No. of Entries in this tree : " << chain->GetEntries() << std::endl;
0095   return kTRUE;
0096 }
0097 
0098 TreeAnalysisHcalScale::TreeAnalysisHcalScale(const char *outFileName, std::vector<std::string>& particles) {
0099   
0100   ipBin = nmaxBin = 0;
0101   double tempgen_Eta[NPBins+1] = {20.0, 30.0, 40.0, 60.0, 100.0, 1000.0};
0102   
0103   for(int i=0; i<NPBins+1; i++)  genPartPBins[i]  = tempgen_Eta[i];
0104   BookHistograms(outFileName, particles);
0105 }
0106 
0107 TreeAnalysisHcalScale::~TreeAnalysisHcalScale() {
0108   fout->cd();
0109   fout->Write();
0110   fout->Close();   
0111 }
0112 
0113 Int_t TreeAnalysisHcalScale::Cut(Long64_t entry) {
0114   // This function may be called from Loop.
0115   // returns  1 if entry is accepted.
0116   // returns -1 otherwise.
0117   return 1;
0118 }
0119 
0120 Int_t TreeAnalysisHcalScale::GetEntry(Long64_t entry) {
0121 
0122   // Read contents of entry.
0123   if (!fChain) return 0;
0124   return fChain->GetEntry(entry);
0125 }
0126 
0127 Long64_t TreeAnalysisHcalScale::LoadTree(Long64_t entry) {
0128 
0129   // Set the environment to read one entry
0130   if (!fChain) return -5;
0131   Long64_t centry = fChain->LoadTree(entry);
0132   if (centry < 0) return centry;
0133   if (!fChain->InheritsFrom(TChain::Class()))  return centry;
0134   TChain *chain = (TChain*)fChain;
0135   if (chain->GetTreeNumber() != fCurrent) {
0136     fCurrent = chain->GetTreeNumber();
0137     Notify();
0138   }
0139   return centry;
0140 }
0141 
0142 void TreeAnalysisHcalScale::Init(TChain *tree) {
0143 
0144   // The Init() function is called when the selector needs to initialize
0145   // a new tree or chain. Typically here the branch addresses and branch
0146   // pointers of the tree will be set.
0147   // It is normally not necessary to make changes to the generated
0148   // code, but the routine can be extended by the user if needed.
0149   // Init() will be called many times when running on PROOF
0150   // (once per file to be processed).
0151   // Set object pointer
0152   t_trackP        = 0;
0153   t_trackPt       = 0;
0154   t_trackEta      = 0;
0155   t_trackPhi      = 0;
0156   t_trackHcalEta  = 0;
0157   t_trackHcalPhi  = 0;
0158   t_hCone         = 0;
0159   t_conehmaxNearP = 0;
0160   t_eMipDR        = 0;
0161   t_eMipDR_2        = 0;
0162   t_eECALDR       = 0;
0163   t_eECALDR_2       = 0;
0164   t_eHCALDR       = 0;
0165   t_e11x11_20Sig  = 0;
0166   t_e15x15_20Sig  = 0;
0167 
0168   // Set branch addresses and branch pointers
0169   if (!tree) return;
0170   fChain = tree;
0171   fCurrent = -1;
0172   fChain->SetMakeClass(1);
0173 
0174   //   fChain->SetBranchAddress("t_EvtNo", &t_EvtNo, &b_t_EvtNo);
0175   fChain->SetBranchAddress("t_RunNo",    &t_RunNo,    &b_t_RunNo);
0176   fChain->SetBranchAddress("t_Lumi",     &t_Lumi,     &b_t_Lumi);
0177   fChain->SetBranchAddress("t_Bunch",    &t_Bunch,    &b_t_Bunch);
0178   fChain->SetBranchAddress("t_trackP",   &t_trackP,   &b_t_trackP);
0179   fChain->SetBranchAddress("t_trackPt",  &t_trackPt,  &b_t_trackPt);
0180   fChain->SetBranchAddress("t_trackEta", &t_trackEta, &b_t_trackEta);
0181   fChain->SetBranchAddress("t_trackPhi", &t_trackPhi, &b_t_trackPhi);
0182   fChain->SetBranchAddress("t_trackHcalEta", &t_trackHcalEta, &b_t_trackHcalEta);
0183   fChain->SetBranchAddress("t_trackHcalPhi", &t_trackHcalPhi, &b_t_trackHcalPhi);
0184   fChain->SetBranchAddress("t_hCone",    &t_hCone,     &b_t_hCone);
0185   fChain->SetBranchAddress("t_conehmaxNearP", &t_conehmaxNearP, &b_t_conehmaxNearP);
0186   fChain->SetBranchAddress("t_eMipDR",   &t_eMipDR,    &b_t_eMipDR);
0187   fChain->SetBranchAddress("t_eMipDR_2",   &t_eMipDR_2,    &b_t_eMipDR_2);
0188   fChain->SetBranchAddress("t_eECALDR",  &t_eECALDR,   &b_t_eECALDR);
0189   fChain->SetBranchAddress("t_eECALDR_2",  &t_eECALDR_2,   &b_t_eECALDR_2);
0190   fChain->SetBranchAddress("t_eHCALDR",  &t_eHCALDR,   &b_t_eHCALDR);
0191   fChain->SetBranchAddress("t_e11x11_20Sig",  &t_e11x11_20Sig,   &b_t_e11x11_20Sig);
0192   fChain->SetBranchAddress("t_e15x15_20Sig",  &t_e15x15_20Sig,   &b_t_e15x15_20Sig);
0193   Notify();
0194 }
0195 
0196 void TreeAnalysisHcalScale::Loop(int cut) {
0197 
0198   if (fChain == 0) return;
0199 
0200   Long64_t nentries = fChain->GetEntriesFast();
0201   std::cout << "No. of Entries in tree " << nentries << std::endl;
0202 
0203   Long64_t nEventsGoodRuns=0, nEventsValidPV=0, nEventsPVTracks=0;
0204 
0205   Long64_t nbytes = 0, nb = 0;
0206   std::map<unsigned int, unsigned int> runEvtList, runNTrkList, runNTrkEtaPList,  runNTrkMipList,  runNTrkMipCharIsoList,  runNIsoTrkList;
0207  nIsoTrkTotal=0;
0208 
0209  int nIsoTrkEtaBin[NEtaBins], nGoodTrkEtaBin[NEtaBins];
0210  for(int i=0; i<NEtaBins; i++){
0211    nIsoTrkEtaBin[i]=0;
0212    nGoodTrkEtaBin[i]=0;
0213  }
0214  for (Long64_t jentry=0; jentry<nentries;jentry++) {
0215    
0216    // load tree and get current entry
0217    Long64_t ientry = LoadTree(jentry);
0218    if (ientry < 0) break;
0219    nb = fChain->GetEntry(jentry);   nbytes += nb;
0220    
0221    if( !(jentry%100000) ) {
0222      std::cout << "procesing event " << jentry+1 << std::endl;
0223    }
0224    
0225    bool goodRun=false;
0226    goodRun = true;
0227    bool evtSel = true;
0228    if( runEvtList.find(t_RunNo) != runEvtList.end() ) {
0229      runEvtList[t_RunNo] += 1; 
0230    } else {
0231      runEvtList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,1) );
0232      std::cout << "runNo " << t_RunNo <<" "<<runEvtList[t_RunNo]<<std::endl;
0233    }
0234    
0235    nEventsGoodRuns++;
0236    
0237    if( runNTrkList.find(t_RunNo) != runNTrkList.end() ) {
0238      runNTrkList[t_RunNo] += t_trackP->size(); 
0239    } else {
0240      runNTrkList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,t_trackP->size()) );
0241      std::cout << "runNo " << t_RunNo <<" "<<runNTrkList[t_RunNo]<<std::endl;
0242    }
0243    
0244    unsigned int NIsoTrk=0, NEtaPTrk=0, NMipTrk=0, NMipCharIsoTrk=0;
0245    for(int itrk=0; itrk<t_trackP->size(); itrk++ ){      
0246      int iTrkEtaBin=-1, iTrkMomBin=-1;
0247      if(std::abs((*t_trackHcalEta)[itrk]) < 26)
0248        iTrkEtaBin = ((*t_trackHcalEta)[itrk] + 25);
0249      else continue;
0250      
0251      for(int ip=0;  ip<NPBins;   ip++)  {
0252        if( (*t_trackP)[itrk]>genPartPBins[ip] &&  (*t_trackP)[itrk]<genPartPBins[ip+1] )  iTrkMomBin = ip;
0253      }
0254      if(iTrkMomBin==2) nGoodTrkEtaBin[iTrkEtaBin]++;
0255      
0256      if(iTrkEtaBin>=0  && iTrkMomBin>=0) {  
0257        NEtaPTrk++;
0258        h_trackP[ipBin]          ->Fill((*t_trackP)[itrk]        );  
0259        h_trackPt[ipBin]         ->Fill((*t_trackPt)[itrk]       );  
0260        h_trackEta[ipBin]        ->Fill((*t_trackEta)[itrk]      );  
0261        h_trackPhi[ipBin]        ->Fill((*t_trackPhi)[itrk]      );  
0262        h_trackHcalEta[ipBin]    ->Fill((*t_trackHcalEta)[itrk]  );  
0263        h_trackHcalPhi[ipBin]    ->Fill((*t_trackHcalPhi)[itrk]  );  
0264        
0265        h_hCone[ipBin]           ->Fill( (*t_hCone)[itrk]         );
0266        h_conehmaxNearP[ipBin]   ->Fill( (*t_conehmaxNearP)[itrk] );
0267        h_eMipDR[ipBin]          ->Fill( (*t_eMipDR)[itrk]        );
0268        h_eECALDR[ipBin]         ->Fill( (*t_eECALDR)[itrk]       );
0269        h_eHCALDR[ipBin]         ->Fill( (*t_eHCALDR)[itrk]       );
0270        h_e11x11_20Sig[ipBin]    ->Fill( (*t_e11x11_20Sig)[itrk]  );
0271        h_e15x15_20Sig[ipBin]    ->Fill( (*t_e15x15_20Sig)[itrk]  ); 
0272        
0273        bool ChargedIso = true, IfMip = true;
0274        bool eNeutIso = true, hNeutIso = true;
0275        
0276        if ((*t_conehmaxNearP)[itrk] >= 0.0) ChargedIso = false;
0277        if ((*t_eMipDR)[itrk] >= 1.0)        IfMip      = false;
0278        
0279        if (std::abs((*t_trackEta)[itrk])<1.47 && ((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]) > 0.5) eNeutIso = false;
0280        else if(((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]) > 2.0) eNeutIso = false;
0281        
0282        if (((*t_eHCALDR)[itrk]-(*t_hCone)[itrk]) > 3) hNeutIso = false;
0283        
0284        if(IfMip){
0285      NMipTrk++;
0286      if(ChargedIso) NMipCharIsoTrk++;
0287        }
0288        
0289        if(eNeutIso && hNeutIso && ChargedIso && IfMip){
0290      h_IsotrackPhi[ipBin]        ->Fill((*t_trackPhi)[itrk]);
0291      if(iTrkMomBin==2) nIsoTrkEtaBin[iTrkEtaBin]++;
0292      NIsoTrk++;
0293      h_IsotrackHcalIEta[ipBin]->Fill((*t_trackHcalEta)[itrk]);
0294      
0295      h_Response[ipBin+1][iTrkMomBin][iTrkEtaBin]        ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0296      h_Response_trunc[ipBin+1][iTrkMomBin][iTrkEtaBin]  ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0297      h_Response_E11x11[ipBin+1][iTrkMomBin][iTrkEtaBin]        ->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0298      h_Response_E11x11_trunc[ipBin+1][iTrkMomBin][iTrkEtaBin]  ->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0299      
0300      h_eHcalFrac[ipBin+1][iTrkMomBin][iTrkEtaBin]       ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0301      h_eHcalFrac[ipBin+1][iTrkMomBin][iTrkEtaBin]       ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0302      h_eHcalFrac_trunc[ipBin+1][iTrkMomBin][iTrkEtaBin] ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0303      h_hneutIso[ipBin+1][iTrkMomBin][iTrkEtaBin]        ->Fill((*t_eHCALDR)[itrk] - (*t_hCone)[itrk]);
0304      h_eneutIso[ipBin+1][iTrkMomBin][iTrkEtaBin]        ->Fill((*t_eECALDR)[itrk] - (*t_eMipDR)[itrk]);
0305      h_eneutIsoNxN[ipBin+1][iTrkMomBin][iTrkEtaBin]     ->Fill((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]);
0306      nIsoTrkTotal++;
0307      if (nIsoTrkTotal <= nmaxBin) {
0308        h_Response[0][iTrkMomBin][iTrkEtaBin]            ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0309        h_Response_trunc[0][iTrkMomBin][iTrkEtaBin]      ->Fill(((*t_hCone)[itrk]+(*t_eMipDR_2)[itrk])/(*t_trackP)[itrk]);
0310        h_Response_E11x11[0][iTrkMomBin][iTrkEtaBin] ->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0311        h_Response_E11x11_trunc[0][iTrkMomBin][iTrkEtaBin]->Fill(((*t_hCone)[itrk]+(*t_e11x11_20Sig)[itrk])/(*t_trackP)[itrk]);
0312        
0313        h_eHcalFrac[0][iTrkMomBin][iTrkEtaBin]           ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0314        h_eHcalFrac_trunc[0][iTrkMomBin][iTrkEtaBin]     ->Fill((*t_hCone)[itrk]/(*t_trackP)[itrk]);
0315        h_hneutIso[0][iTrkMomBin][iTrkEtaBin]            ->Fill((*t_eHCALDR)[itrk] - (*t_hCone)[itrk]);
0316        h_eneutIso[0][iTrkMomBin][iTrkEtaBin]            ->Fill((*t_eECALDR)[itrk] - (*t_eMipDR)[itrk]);
0317        h_eneutIsoNxN[0][iTrkMomBin][iTrkEtaBin]         ->Fill((*t_e15x15_20Sig)[itrk]-(*t_e11x11_20Sig)[itrk]);
0318      }
0319        }
0320      }
0321    } // loop over tracks in the event 
0322    if( runNIsoTrkList.find(t_RunNo) != runNIsoTrkList.end() ) {
0323      runNIsoTrkList[t_RunNo] += NIsoTrk;
0324    } else {
0325      runNIsoTrkList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NIsoTrk) );
0326      std::cout << "runNo " << t_RunNo <<" "<<runNIsoTrkList[t_RunNo]<<std::endl;
0327    }
0328    if( runNTrkEtaPList.find(t_RunNo) != runNTrkEtaPList.end() ) {
0329      runNTrkEtaPList[t_RunNo] += NEtaPTrk;
0330    } else {
0331      runNTrkEtaPList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NEtaPTrk) );
0332      std::cout << "runNo " << t_RunNo <<" "<<runNTrkEtaPList[t_RunNo]<<std::endl;
0333    }
0334    
0335     if( runNTrkMipList.find(t_RunNo) != runNTrkMipList.end() ) {
0336       runNTrkMipList[t_RunNo] += NMipTrk;
0337     } else {
0338       runNTrkMipList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NMipTrk) );
0339       std::cout << "runNo " << t_RunNo <<" "<<runNTrkMipList[t_RunNo]<<std::endl;
0340     }
0341     
0342     if( runNTrkMipCharIsoList.find(t_RunNo) != runNTrkMipCharIsoList.end() ) {
0343       runNTrkMipCharIsoList[t_RunNo] += NMipCharIsoTrk;
0344     } else {
0345       runNTrkMipCharIsoList.insert( std::pair<unsigned int, unsigned int>(t_RunNo,NMipCharIsoTrk) );
0346       std::cout << "runNo " << t_RunNo <<" "<<runNTrkMipCharIsoList[t_RunNo]<<std::endl;
0347     }
0348           
0349   
0350   } //loop over tree entries
0351 
0352   std::cout <<"Number of entries in tree "<< nentries <<" # of isolated tracks "<< nIsoTrkTotal 
0353         <<" should be at least " << nmaxBin << std::endl;
0354   
0355  
0356   std::cout << "saved runEvtList " << std::endl;
0357   std::map<unsigned int, unsigned int>::iterator runEvtListItr = runEvtList.begin();
0358   for(runEvtListItr=runEvtList.begin(); runEvtListItr != runEvtList.end(); runEvtListItr++) {
0359     std::cout<<runEvtListItr->first << " "<< runEvtListItr->second << std::endl;
0360   }
0361 
0362   std::cout << "Number of tracks in runs " << std::endl;
0363   std::map<unsigned int, unsigned int>::iterator runNTrkListItr = runNTrkList.begin();
0364   for(runNTrkListItr=runNTrkList.begin(); runNTrkListItr != runNTrkList.end(); runNTrkListItr++) {
0365     std::cout<<runNTrkListItr->first << " "<< runNTrkListItr->second << std::endl;
0366   }
0367   
0368   std::cout << "Number of tracks in eta range in runs " << std::endl;
0369   std::map<unsigned int, unsigned int>::iterator runNTrkEtaPListItr = runNTrkEtaPList.begin();
0370   for(runNTrkEtaPListItr=runNTrkEtaPList.begin(); runNTrkEtaPListItr != runNTrkEtaPList.end(); runNTrkEtaPListItr++) {
0371     std::cout<<runNTrkEtaPListItr->first << " "<< runNTrkEtaPListItr->second << std::endl;
0372   }
0373   
0374   std::cout << "Number of Mip tracks in runs " << std::endl;
0375   std::map<unsigned int, unsigned int>::iterator runNTrkMipListItr = runNTrkMipList.begin();
0376   for(runNTrkMipListItr=runNTrkMipList.begin(); runNTrkMipListItr != runNTrkMipList.end(); runNTrkMipListItr++) {
0377     std::cout<<runNTrkMipListItr->first << " "<< runNTrkMipListItr->second << std::endl;
0378    }
0379 
0380 
0381   std::cout << "Number of Charged isolated Mip tracks in runs " << std::endl;
0382   std::map<unsigned int, unsigned int>::iterator runNTrkMipCharIsoListItr = runNTrkMipCharIsoList.begin();
0383   for(runNTrkMipCharIsoListItr=runNTrkMipCharIsoList.begin(); runNTrkMipCharIsoListItr != runNTrkMipCharIsoList.end(); runNTrkMipCharIsoListItr++) {
0384     std::cout<<runNTrkMipCharIsoListItr->first << " "<< runNTrkMipCharIsoListItr->second << std::endl;
0385    }
0386 
0387   std::cout << "Number of isolated tracks in runs " << std::endl;
0388   std::map<unsigned int, unsigned int>::iterator runNIsoTrkListItr = runNIsoTrkList.begin();
0389   for(runNIsoTrkListItr=runNIsoTrkList.begin(); runNIsoTrkListItr != runNIsoTrkList.end(); runNIsoTrkListItr++) {
0390     std::cout<<runNIsoTrkListItr->first << " "<< runNIsoTrkListItr->second << std::endl;
0391   } 
0392   for(int i=0; i<NEtaBins;i++){
0393     std::cout<< "Number of tracks in ieta " << i-25 << ": " <<  nIsoTrkEtaBin[i] << std::endl;
0394     if(nGoodTrkEtaBin[i]!=0){
0395       double frac = (double)nIsoTrkEtaBin[i]/(double)nGoodTrkEtaBin[i];
0396       h_FracIsotrackHcalIEta[ipBin]->Fill(i-25, frac);
0397     }
0398   }
0399 }
0400 
0401 Bool_t TreeAnalysisHcalScale::Notify() {
0402   // The Notify() function is called when a new file is opened. This
0403   // can be either for a new TTree in a TChain or when when a new TTree
0404   // is started when using PROOF. It is normally not necessary to make changes
0405   // to the generated code, but the routine can be extended by the
0406   // user if needed. The return value is currently not used.
0407 
0408   return kTRUE;
0409 }
0410 
0411 void TreeAnalysisHcalScale::Show(Long64_t entry) {
0412   // Print contents of entry.
0413   // If entry is not specified, print current entry
0414   if (!fChain) return;
0415   fChain->Show(entry);
0416 }
0417 
0418 void TreeAnalysisHcalScale::BookHistograms(const char *outFileName, std::vector<std::string>& particles) {
0419 
0420   fout = new TFile(outFileName, "RECREATE");
0421 
0422   fout->cd();
0423 
0424   char hname[1000], htit[1000];
0425 
0426   for (unsigned int i=0; i<particles.size(); i++) {
0427     sprintf (hname, "trackEtaAll%s", particles[i].c_str());
0428     sprintf (htit,  "Eta:    All Tracks for %s", particles[i].c_str());
0429     h_trackEtaAll[i]     = new TH1F(hname, htit,    100,  -3.0,   3.0    );
0430     h_trackEtaAll[i]     ->Sumw2();
0431     sprintf (hname, "trackP%s", particles[i].c_str());
0432     sprintf (htit,  "P:    Tracks in eta region for %s", particles[i].c_str());
0433     h_trackP[i]          = new TH1F(hname, htit,   5000,  0.0,   1000   );
0434     h_trackP[i]          ->Sumw2();
0435     sprintf (hname, "trackPt%s", particles[i].c_str());
0436     sprintf (htit,  "Pt:   Tracks in eta region for %s", particles[i].c_str());
0437     h_trackPt[i]         = new TH1F(hname, htit,   5000,  0.0,   1000   );
0438     h_trackPt[i]         ->Sumw2();
0439     sprintf (hname, "trackEta%s", particles[i].c_str());
0440     sprintf (htit,  "Eta:  Tracks in eta region for %s", particles[i].c_str());
0441     h_trackEta[i]        = new TH1F(hname, htit,    100,  -3.0,   3.0    );
0442     h_trackEta[i]        ->Sumw2();
0443     sprintf (hname, "trackPhi%s", particles[i].c_str());
0444     sprintf (htit,  "Phi:  Tracks in eta region for %s", particles[i].c_str());
0445     h_trackPhi[i]        = new TH1F(hname, htit,    100,  -4.0,   4.0    );
0446     h_trackPhi[i]       ->Sumw2();
0447     sprintf (hname, "IsotrackPhi%s", particles[i].c_str());
0448     sprintf (htit,  "Phi:  Isolated tracks for %s", particles[i].c_str());
0449     h_IsotrackPhi[i]        = new TH1F(hname, htit,    100,  -4.0,   4.0    );
0450     h_IsotrackPhi[i]       ->Sumw2();
0451     sprintf (hname, "trackHcalEta%s", particles[i].c_str());
0452     sprintf (htit,  "iEta:  Track Hit at Hcal for %s", particles[i].c_str());
0453     h_trackHcalEta[i]    = new TH1F(hname, htit,    100,  -50.0,  50.0   );
0454     h_trackHcalEta[i]    ->Sumw2();
0455     sprintf (hname, "IsotrackHcalIEta%s", particles[i].c_str());
0456     sprintf (htit,  "iEta:  Isolated Track Hit at Hcal for %s", particles[i].c_str());
0457     h_IsotrackHcalIEta[i]    = new TH1F(hname, htit,    60,  -30.0,  30.0   );
0458     h_IsotrackHcalIEta[i]    ->Sumw2();
0459     sprintf (hname, "FracIsotrackHcalIEta%s", particles[i].c_str());
0460     sprintf (htit,  "iEta:  Fraction of Isolated Track Hit at Hcal for %s", particles[i].c_str());
0461     h_FracIsotrackHcalIEta[i]    = new TH1F(hname, htit,    60,  -30.0,  30.0   );
0462     h_FracIsotrackHcalIEta[i]    ->Sumw2();
0463 
0464     sprintf (hname, "trackHcalPhi%s", particles[i].c_str());
0465     sprintf (htit,  "iPhi:  Track hit at Hcal for %s", particles[i].c_str());
0466     h_trackHcalPhi[i]    = new TH1F(hname, htit,     100,   0.0,   100.   );
0467     h_trackHcalPhi[i]    ->Sumw2();
0468     sprintf (hname, "h_hcone%s", particles[i].c_str());
0469     sprintf (htit,  "Energy in Hcal in the cone for %s", particles[i].c_str());
0470     h_hCone[i]           = new TH1F(hname, htit,    5000, -2.0,   1000.0 );
0471     h_hCone[i]          ->Sumw2();
0472     sprintf (hname, "h_conehmaxNearP%s", particles[i].c_str());
0473     sprintf (htit,  "Max energy in charge isolation for %s", particles[i].c_str());
0474     h_conehmaxNearP[i]   = new TH1F(hname, htit,    5000, -2.0,   1000.0 );
0475     h_conehmaxNearP[i]   ->Sumw2();
0476     sprintf (hname, "h_eMipDR%s", particles[i].c_str());
0477     sprintf (htit,  "Energy in the MIP region for %s", particles[i].c_str());
0478     h_eMipDR[i]          = new TH1F(hname, htit,    5000, -2.0,   1000.0 );
0479     h_eMipDR[i]         ->Sumw2();
0480     sprintf (hname, "h_eECALDR%s", particles[i].c_str());
0481     sprintf (htit,  "Energy in ECAL iso region for %s", particles[i].c_str());
0482     h_eECALDR[i]         = new TH1F(hname, htit,    5000, -2.0,   1000.0 );
0483     h_eECALDR[i]         ->Sumw2();
0484     sprintf (hname, "h_eHCALDR%s", particles[i].c_str());
0485     sprintf (htit,  "Energy in HCAL iso region for %s", particles[i].c_str());
0486     h_eHCALDR[i]         = new TH1F(hname, htit,    5000, -2.0,   1000.0 );
0487     h_eHCALDR[i]         ->Sumw2();
0488     sprintf (hname, "h_e11x11_20Sig%s", particles[i].c_str());
0489     sprintf (htit,  "Energy in 11x11 for %s", particles[i].c_str());
0490     h_e11x11_20Sig[i]    = new TH1F(hname ,htit,    5000, -2.0,   1000.0 );
0491     h_e11x11_20Sig[i]    ->Sumw2();
0492     sprintf (hname, "h_e15x15_20Sig%s", particles[i].c_str());
0493     sprintf (htit,  "Energy in 15x15 for %s", particles[i].c_str());
0494     h_e15x15_20Sig[i]    = new TH1F(hname, htit,    5000, -2.0,   1000.0 );
0495     h_e15x15_20Sig[i]    ->Sumw2();
0496   }
0497 
0498   TDirectory *d_HcalFrac     = fout->mkdir( "HcalFrac"    );
0499   TDirectory *d_eNeutIso     = fout->mkdir( "eNeutIso"    );
0500   TDirectory *d_hNeutIso     = fout->mkdir( "hNeutIso"    );
0501   TDirectory *d_eNeutIsoNxN  = fout->mkdir( "eNeutIsoNxN" );
0502   TDirectory *d_Response     = fout->mkdir( "Response" );
0503 
0504   for (unsigned int i=0; i<particles.size()+1; i++) {
0505     char part[10];
0506     if (i == 0) sprintf (part, "all");
0507     else        sprintf (part, "%s", particles[i-1].c_str());
0508     for (int ieta=0; ieta<NEtaBins; ieta++) {
0509       int iEta = ieta;
0510       for (int ip=0; ip<NPBins; ip++) {
0511     double lowMom=-5.0, highMom= 5.0;
0512     lowMom  = genPartPBins[ip];
0513     highMom = genPartPBins[ip+1];
0514           
0515     d_Response->cd();
0516     sprintf(hname, "h_Response_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0517     sprintf(htit,  "Response (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0518     h_Response[i][ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0519     h_Response[i][ip][ieta] ->Sumw2();
0520 
0521     sprintf(hname, "h_Response_trunc_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0522     sprintf(htit,  "Response_trunc (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0523     h_Response_trunc[i][ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0524     h_Response_trunc[i][ip][ieta] ->Sumw2();
0525 
0526     sprintf(hname, "h_Response_E11x11_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0527     sprintf(htit,  "Response_E11x11 (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0528     h_Response_E11x11[i][ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0529     h_Response_E11x11[i][ip][ieta] ->Sumw2();
0530 
0531     sprintf(hname, "h_Response_E11x11_trunc_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0532     sprintf(htit,  "Response_E11x11_trunc (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0533     h_Response_E11x11_trunc[i][ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0534     h_Response_E11x11_trunc[i][ip][ieta] ->Sumw2();
0535 
0536     d_HcalFrac->cd();
0537     sprintf(hname, "h_eHcalFrac_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0538     sprintf(htit,  "eHcalFrac (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0539     h_eHcalFrac[i][ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0540     h_eHcalFrac[i][ip][ieta] ->Sumw2();
0541 
0542     sprintf(hname, "h_eHcalFrac_trunc_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0543     sprintf(htit,  "eHcalFrac_trunc (|i#eta|=%i),(%3.2f<|#p|<%3.2f) for %s", iEta, lowMom, highMom, part);
0544     h_eHcalFrac_trunc[i][ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0545     h_eHcalFrac_trunc[i][ip][ieta] ->Sumw2();
0546       
0547     d_hNeutIso->cd();
0548     sprintf(hname, "h_hneutIso_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0549     sprintf(htit,  "hneutIso (|i#eta|=%i),(%3.2f<|#p|<%3.2f for %s)", iEta, lowMom, highMom, part);
0550     h_hneutIso[i][ip][ieta] = new TH1F(hname, htit, 1000, -5.0, 20.0);
0551     h_hneutIso[i][ip][ieta] ->Sumw2();
0552       
0553     d_eNeutIso->cd();
0554     sprintf(hname, "h_eneutIso_etaBin%i_pBin%i_particle%s", ieta,ip,part);
0555     sprintf(htit,  "eneutIso (|i#eta|=%i),(%3.2f<|#p|<%3.2f for %s)", iEta, lowMom, highMom, part);
0556     h_eneutIso[i][ip][ieta] = new TH1F(hname, htit, 1000, -1.0, 25.0);
0557     h_eneutIso[i][ip][ieta] ->Sumw2();
0558       
0559     d_eNeutIsoNxN->cd();
0560     sprintf(hname, "h_eneutIsoNxN_etaBin%i_pBin%i_particle%s",ieta,ip,part);
0561     sprintf(htit,  "eneutIsoNxN (|i#eta|=%i),(%3.2f<|#p|<%3.2f for %s)", iEta, lowMom, highMom, part);
0562     h_eneutIsoNxN[i][ip][ieta] = new TH1F(hname, htit, 1000, -1.0, 10.0);
0563     h_eneutIsoNxN[i][ip][ieta] ->Sumw2();
0564       }
0565     }  
0566   }
0567   for (int ieta=0; ieta<NEtaBins; ieta++) {
0568     int iEta = ieta;
0569     for (int ip=0; ip<NPBins; ip++) {
0570       double lowMom=-5.0, highMom= 5.0;
0571       lowMom  = genPartPBins[ip];
0572       highMom = genPartPBins[ip+1];
0573       
0574       d_HcalFrac->cd();
0575       sprintf(hname, "h_eHcalFrac_all_etaBin%i_pBin%i", ieta,ip);
0576       sprintf(htit,  "eHcalFrac allWeighted (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0577       h_eHcalFrac_all[ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0578       h_eHcalFrac_all[ip][ieta] ->Sumw2();
0579       
0580       sprintf(hname, "h_eHcalFrac_trunc_all_etaBin%i_pBin%i", ieta,ip);
0581       sprintf(htit,  "eHcalFrac_trunc allWeigthed (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0582       h_eHcalFrac_trunc_all[ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0583       h_eHcalFrac_trunc_all[ip][ieta] ->Sumw2();
0584 
0585       d_Response->cd();
0586       sprintf(hname, "h_Response_all_etaBin%i_pBin%i", ieta,ip);
0587       sprintf(htit,  "Response allWeighted (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0588       h_Response_all[ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0589       h_Response_all[ip][ieta] ->Sumw2();
0590       
0591       sprintf(hname, "h_Response_trunc_all_etaBin%i_pBin%i", ieta,ip);
0592       sprintf(htit,  "Response_trunc allWeigthed (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0593       h_Response_trunc_all[ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0594       h_Response_trunc_all[ip][ieta] ->Sumw2();
0595 
0596       sprintf(hname, "h_Response_E11x11_all_etaBin%i_pBin%i", ieta,ip);
0597       sprintf(htit,  "Response_E11x11 allWeighted (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0598       h_Response_E11x11_all[ip][ieta] = new TH1F(hname, htit, 500, -1.0, 4.0);
0599       h_Response_E11x11_all[ip][ieta] ->Sumw2();
0600       
0601       sprintf(hname, "h_Response_E11x11_trunc_all_etaBin%i_pBin%i", ieta,ip);
0602       sprintf(htit,  "Response_E11x11_trunc allWeigthed (|i#eta|=%i),(%3.2f<|#p|<%3.2f)", iEta, lowMom, highMom);
0603       h_Response_E11x11_trunc_all[ip][ieta] = new TH1F(hname, htit, 500, 0.2, 4.0);
0604       h_Response_E11x11_trunc_all[ip][ieta] ->Sumw2();
0605 
0606     }
0607   }
0608   fout->cd();
0609   
0610 }
0611               
0612 void TreeAnalysisHcalScale::clear() {
0613    if (!fChain) return;
0614    delete fChain->GetCurrentFile();
0615 }
0616               
0617 void TreeAnalysisHcalScale::setParticle(unsigned int ip, unsigned int nmax) {
0618   ipBin = ip;
0619   nmaxBin = nmax;
0620 }
0621 
0622 void TreeAnalysisHcalScale::AddWeight( std::vector<std::string> particleNames){
0623 
0624  for (unsigned int i=0; i<particleNames.size(); i++) {
0625    char fileList[200], treeName[200];
0626    std::cout << "particle " << i << " weight " << weights[i] << std::endl; 
0627    for (int ieta=0; ieta<NEtaBins; ieta++) {
0628      for (int ip=0; ip<NPBins; ip++) {
0629        h_eHcalFrac_all[ip][ieta]->Add(h_eHcalFrac[i+1][ip][ieta], weights[i]);
0630        h_eHcalFrac_trunc_all[ip][ieta]->Add(h_eHcalFrac_trunc[i+1][ip][ieta], weights[i]);
0631 
0632        h_Response_all[ip][ieta]->Add(h_Response[i+1][ip][ieta], weights[i]);
0633        h_Response_trunc_all[ip][ieta]->Add(h_Response_trunc[i+1][ip][ieta], weights[i]);
0634 
0635        h_Response_E11x11_all[ip][ieta]->Add(h_Response_E11x11[i+1][ip][ieta], weights[i]);
0636        h_Response_E11x11_trunc_all[ip][ieta]->Add(h_Response_E11x11_trunc[i+1][ip][ieta], weights[i]);
0637 
0638      }
0639    }
0640  }
0641 }
0642
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!