Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:59 

0001 //////////////////////////////////////////////////////////
0002 // This class has been automatically generated on
0003 // Tue May  2 15:43:18 2017 by ROOT version 5.34/19
0004 // from TTree RecJetHF/RecJetHF Tree
0005 // found on file: HCALHF.root
0006 //////////////////////////////////////////////////////////
0007 
0008 #include <TCanvas.h>
0009 #include <TChain.h>
0010 #include <TDirectory.h>
0011 #include <TFile.h>
0012 #include <TGraph.h>
0013 #include <TH1D.h>
0014 #include <TH2D.h>
0015 #include <TPaveStats.h>
0016 #include <TROOT.h>
0017 #include <TStyle.h>
0018 #include <TTree.h>
0019 #include <iostream>
0020 #include <string>
0021 #include <cmath>
0022 #include <map>
0023 #include <fstream>
0024 #include <TF1.h>
0025 #include "TLegend.h"
0026 
0027 #include <iomanip>
0028 #include <sstream>
0029 
0030 class RecJetHF {
0031 public:
0032   struct MyInfo {
0033     double Mom0_F1, Mom1_F1, Mom2_F1,Mom0_F2, Mom1_F2, Mom2_F2;
0034 
0035     MyInfo() {
0036       Mom0_F1 = Mom1_F1 = Mom2_F1 = Mom0_F2 = Mom1_F2 = Mom2_F2 =0.;
0037     }
0038   };
0039 
0040   struct Hists {
0041     TH1D *h1, *h2, *h3, *h4;
0042     Hists() {
0043       h1 = h2 =h3 =h4 =0;
0044     }
0045   };
0046 
0047   RecJetHF(std::string fname, bool ratio=false);
0048   virtual ~RecJetHF();
0049   virtual Int_t    Cut(Long64_t entry);
0050   virtual Int_t    GetEntry(Long64_t entry);
0051   virtual Long64_t LoadTree(Long64_t entry);
0052   virtual void     Init(TTree *tree);
0053   virtual void     Loop();
0054   virtual std::map<unsigned int,RecJetHF::MyInfo> LoopMap();
0055   virtual Bool_t   Notify();
0056   virtual void     Show(Long64_t entry = -1);
0057   std::map<unsigned int, RecJetHF::Hists> MakeRatio(std::map<unsigned int,RecJetHF::MyInfo>& m_);
0058 
0059 private :
0060   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
0061   Int_t           fCurrent; //!current Tree number in a TChain
0062 
0063   // Declaration of leaf types
0064   Int_t           cells;
0065   Int_t           mysubd;
0066   Int_t           depth;
0067   Int_t           ieta;
0068   Int_t           iphi;
0069   Float_t         mom0_F1;
0070   Float_t         mom1_F1;
0071   Float_t         mom2_F1;
0072   Float_t         mom3_F1;
0073   Float_t         mom4_F1;
0074   Float_t         mom0_F2;
0075   Float_t         mom1_F2;
0076   Float_t         mom2_F2;
0077   Float_t         mom3_F2;
0078   Float_t         mom4_F2;
0079   Int_t           trigbit;
0080   Double_t        rnnumber;
0081   
0082   // List of branches
0083   TBranch        *b_cells;   //!
0084   TBranch        *b_mysubd;   //!
0085   TBranch        *b_depth;   //!
0086   TBranch        *b_ieta;   //!
0087   TBranch        *b_iphi;   //!
0088   TBranch        *b_mom0_F1;   //!
0089   TBranch        *b_mom1_F1;   //!
0090   TBranch        *b_mom2_F1;   //!
0091   TBranch        *b_mom3_F1;   //!
0092   TBranch        *b_mom4_F1;   //!
0093   TBranch        *b_mom0_F2;   //!
0094   TBranch        *b_mom1_F2;   //!
0095   TBranch        *b_mom2_F2;   //!
0096   TBranch        *b_mom3_F2;   //!
0097   TBranch        *b_mom4_F2;   //!
0098   TBranch        *b_trigbit;   //!
0099   TBranch        *b_rnnumber;   //!
0100   
0101   TFile          *file;
0102   bool            ratio_;
0103 
0104 };
0105 
0106 RecJetHF::RecJetHF(std::string fname, bool ratio) : fChain(0), ratio_(ratio) {
0107 
0108   file = new TFile(fname.c_str());
0109   TDirectory* directory = (TDirectory*)(file->FindObjectAny("recAnalyzerHF"));
0110   TTree *tree = (TTree*)(directory->FindObjectAny("RecJet"));
0111   Init(tree);
0112 }
0113 
0114 RecJetHF::~RecJetHF() {
0115   if (!fChain) return;
0116   delete fChain->GetCurrentFile();
0117 }
0118 
0119 Int_t RecJetHF::GetEntry(Long64_t entry) {
0120   // Read contents of entry.
0121   if (!fChain) return 0;
0122   return fChain->GetEntry(entry);
0123 }
0124 
0125 Long64_t RecJetHF::LoadTree(Long64_t entry) {
0126   // Set the environment to read one entry
0127   if (!fChain) return -5;
0128   Long64_t centry = fChain->LoadTree(entry);
0129   if (centry < 0) return centry;
0130   if (fChain->GetTreeNumber() != fCurrent) {
0131     fCurrent = fChain->GetTreeNumber();
0132     Notify();
0133   }
0134   return centry;
0135 }
0136 
0137 void RecJetHF::Init(TTree *tree) {
0138   // The Init() function is called when the selector needs to initialize
0139   // a new tree or chain. Typically here the branch addresses and branch
0140   // pointers of the tree will be set.
0141   // It is normally not necessary to make changes to the generated
0142   // code, but the routine can be extended by the user if needed.
0143   // Init() will be called many times when running on PROOF
0144   // (once per file to be processed).
0145 
0146   // Set branch addresses and branch pointers
0147   if (!tree) return;
0148   fChain = tree;
0149   fCurrent = -1;
0150   fChain->SetMakeClass(1);
0151 
0152   fChain->SetBranchAddress("cells", &cells, &b_cells);
0153   fChain->SetBranchAddress("mysubd", &mysubd, &b_mysubd);
0154   fChain->SetBranchAddress("depth", &depth, &b_depth);
0155   fChain->SetBranchAddress("ieta", &ieta, &b_ieta);
0156   fChain->SetBranchAddress("iphi", &iphi, &b_iphi);
0157   fChain->SetBranchAddress("mom0_F1", &mom0_F1, &b_mom0_F1);
0158   fChain->SetBranchAddress("mom1_F1", &mom1_F1, &b_mom1_F1);
0159   fChain->SetBranchAddress("mom2_F1", &mom2_F1, &b_mom2_F1);
0160   fChain->SetBranchAddress("mom3_F1", &mom3_F1, &b_mom3_F1);
0161   fChain->SetBranchAddress("mom4_F1", &mom4_F1, &b_mom4_F1);
0162   fChain->SetBranchAddress("mom0_F2", &mom0_F2, &b_mom0_F2);
0163   fChain->SetBranchAddress("mom1_F2", &mom1_F2, &b_mom1_F2);
0164   fChain->SetBranchAddress("mom2_F2", &mom2_F2, &b_mom2_F2);
0165   fChain->SetBranchAddress("mom3_F2", &mom3_F2, &b_mom3_F2);
0166   fChain->SetBranchAddress("mom4_F2", &mom4_F2, &b_mom4_F2);
0167   fChain->SetBranchAddress("trigbit", &trigbit, &b_trigbit);
0168   fChain->SetBranchAddress("rnnumber", &rnnumber, &b_rnnumber);
0169   Notify();
0170 }
0171 
0172 Bool_t RecJetHF::Notify() {
0173   // The Notify() function is called when a new file is opened. This
0174   // can be either for a new TTree in a TChain or when when a new TTree
0175   // is started when using PROOF. It is normally not necessary to make changes
0176   // to the generated code, but the routine can be extended by the
0177   // user if needed. The return value is currently not used.
0178 
0179   return kTRUE;
0180 }
0181 
0182 void RecJetHF::Show(Long64_t entry) {
0183   // Print contents of entry.
0184   // If entry is not specified, print current entry
0185   if (!fChain) return;
0186   fChain->Show(entry);
0187 }
0188 
0189 Int_t RecJetHF::Cut(Long64_t ) {
0190   // This function may be called from Loop.
0191   // returns  1 if entry is accepted.
0192   // returns -1 otherwise.
0193   return 1;
0194 }
0195 
0196 std::map<unsigned int,RecJetHF::MyInfo> RecJetHF::LoopMap() {
0197   //   In a ROOT session, you can do: 
0198   //      Root > .L RecJetHF.C
0199   //      Root > RecJetHF t
0200   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0201   //      Root > t.Show();       // Show values of entry 12
0202   //      Root > t.Show(16);     // Read and show values of entry 16
0203   //      Root > t.Loop();       // Loop on all entries
0204   //
0205   //     This is the loop skeleton where:
0206   //    jentry is the global entry number in the chain
0207   //    ientry is the entry number in the current Tree
0208   //  Note that the argument to GetEntry must be:
0209   //    jentry for TChain::GetEntry
0210   //    ientry for TTree::GetEntry and TBranch::GetEntry
0211   //
0212   //       To read only selected branches, Insert statements like:
0213   // METHOD1:
0214   //    fChain->SetBranchStatus("*",0);  // disable all branches
0215   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0216   // METHOD2: replace line
0217   //    fChain->GetEntry(jentry);       //read all branches
0218   //by  b_branchname->GetEntry(ientry); //read only this branch
0219   
0220   std::map<unsigned int,RecJetHF::MyInfo> m_;
0221   if (fChain != 0) {
0222     Long64_t nentries = fChain->GetEntriesFast();
0223     Long64_t nbytes = 0, nb = 0;
0224     for (Long64_t jentry=0; jentry<nentries;jentry++) {
0225       Long64_t ientry = LoadTree(jentry);
0226       if (ientry < 0) break;
0227       
0228       nb = fChain->GetEntry(jentry);   nbytes += nb;
0229       unsigned int detId1 = ((mysubd<<20) | ((depth&0x1f)<<14) |
0230                              ((ieta>0)?(0x2000|(ieta<<7)):((-ieta)<<7)) |
0231                              (iphi&0x7f));
0232       std::map<unsigned int, RecJetHF::MyInfo>::iterator mitr = m_.find(detId1);
0233       if (mitr == m_.end()) {
0234     RecJetHF::MyInfo info;
0235         m_[detId1] = info ;
0236         mitr = m_.find(detId1) ;
0237       }
0238       mitr->second.Mom0_F1 += mom0_F1;
0239       mitr->second.Mom1_F1 += mom1_F1;
0240       mitr->second.Mom2_F1 += mom2_F1;
0241       mitr->second.Mom0_F2 += mom0_F2;
0242       mitr->second.Mom1_F2 += mom1_F2;
0243       mitr->second.Mom2_F2 += mom2_F2;
0244     }
0245   }
0246   return m_;
0247 }
0248 
0249 void RecJetHF::Loop() {
0250 
0251   std::map<unsigned int,RecJetHF::MyInfo> m_ = LoopMap();
0252   if (m_.size() == 0) return;
0253 
0254   char fname[80];
0255   
0256   sprintf (fname, "SingleNeutrino2017HF_Prereco.root");
0257   TFile f(fname, "recreate");
0258   std::map<unsigned int,RecJetHF::Hists> hh_ = MakeRatio(m_);
0259   
0260   for (std::map<unsigned int,RecJetHF::Hists>::const_iterator hitr = hh_.begin();
0261        hitr != hh_.end(); ++hitr) {
0262     if (hitr->second.h1  != 0) hitr->second.h1->Write();
0263     if (hitr->second.h2  != 0) hitr->second.h2->Write();
0264     //    if (hitr->second.h3  != 0) hitr->second.h3->Write();
0265     //if (hitr->second.h4  != 0) hitr->second.h4->Write();
0266 
0267   }
0268   f.Close();
0269 }
0270 
0271 std::map<unsigned int, RecJetHF::Hists>
0272 RecJetHF::MakeRatio(std::map<unsigned int,RecJetHF::MyInfo>& m_) {
0273 
0274   char name[80], title[80];
0275   int keta(0), kphi(0);
0276 
0277   std::map<unsigned int, RecJetHF::Hists> hh_ ;
0278   for (std::map<unsigned int,RecJetHF::MyInfo>::iterator mitr = m_.begin();
0279        mitr !=m_.end(); mitr++) {
0280     unsigned int sdet0  = ((mitr->first)>>20)&0x7;
0281 
0282     if (sdet0 == 4) {
0283       unsigned int detId2  = (mitr->first)&0x7fff80 ;
0284       std::map<unsigned int,RecJetHF::Hists>::iterator hitr = hh_.find(detId2);
0285       if (hitr == hh_.end()) {
0286     RecJetHF::Hists hh;
0287         hh_[detId2] = hh;
0288         hitr = hh_.find(detId2);
0289         keta = (((mitr->first)&0x2000) ? (((mitr->first)>>7)&0x3f) :
0290         -(((mitr->first)>>7)&0x3f));
0291         int dept = (((mitr->first)>>14)&0x1f);
0292     sprintf (name,  "histE1%d %d", keta, dept);
0293     if (ratio_) 
0294       sprintf (title, "<E1/(E1+E2)> (i#eta %d depth %d)", keta, dept);
0295     else
0296       sprintf (title, "<E1>/(<E1>+<E2>) (i#eta %d depth %d)", keta, dept);
0297         hitr->second.h1 = new TH1D(name, title, 72,0.,72.);
0298         hitr->second.h1->GetXaxis()->SetTitle("i#phi");
0299         hitr->second.h1->GetXaxis()->CenterTitle();
0300     if (ratio_) 
0301       hitr->second.h1->GetYaxis()->SetTitle("<E1/(E1+E2)>");
0302     else
0303       hitr->second.h1->GetYaxis()->SetTitle("<E1>/(<E1>+<E2>)");
0304         hitr->second.h1->GetYaxis()->CenterTitle();
0305     sprintf (name,  "histE2%d %d", keta, dept);
0306     if (ratio_) 
0307       sprintf (title, "<E2/(E1+<E2)> (i#eta %d depth %d)", keta, dept);
0308     else
0309       sprintf (title, "<E2>/(<E1>+<E2>) (i#eta %d depth %d)", keta, dept);
0310     hitr->second.h2 = new TH1D(name, title, 72,0.,72.);
0311         hitr->second.h2->GetXaxis()->SetTitle("i#phi");
0312         hitr->second.h2->GetXaxis()->CenterTitle();
0313     if (ratio_) 
0314       hitr->second.h2->GetYaxis()->SetTitle("<E2/(E1+E2)>");
0315     else
0316       hitr->second.h2->GetYaxis()->SetTitle("<E2>/(<E1>+<E2>)");
0317         hitr->second.h2->GetYaxis()->CenterTitle();
0318       }
0319       kphi = ((mitr->first)&0x7f); 
0320 
0321       double E1 = mitr->second.Mom1_F1;
0322       double E2 = mitr->second.Mom1_F2;
0323       double mom0 = mitr->second.Mom0_F1;
0324       double mom2_E1 = mitr->second.Mom2_F1;
0325       double mom2_E2 = mitr->second.Mom2_F2;
0326       double err_E1 = std::sqrt((mom2_E1/mom0 - (E1*E1)/(mom0*mom0))/mom0);
0327       double err_E2 = std::sqrt((mom2_E2/mom0 - (E2*E2)/(mom0*mom0))/mom0);
0328       double val1   = ratio_ ? E1/mom0 : E1/(E1+E2);
0329       double val2   = ratio_ ? E2/mom0 : E2/(E1+E2);
0330       double err1   = ratio_ ? err_E1 :
0331     ((E1+E2)* err_E1 - E1*(err_E1+err_E2)) / (pow((E1+E2),2)*mom0);
0332       double err2   = ratio_ ? err_E2 :
0333     ((E1+E2)* err_E2 - E2*(err_E1+err_E2)) / (pow((E1+E2),2)*mom0);
0334       hitr->second.h1->SetBinContent(kphi,val1);
0335       hitr->second.h1->SetBinError(kphi,err1);
0336       hitr->second.h2->SetBinContent(kphi,val2);
0337       hitr->second.h2->SetBinError(kphi,err2);
0338     }
0339   }
0340   return hh_;
0341 }
[ 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!