From ec4b6a6f9de9751cc882c2e1f52fd08009bdbf11 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Fri, 17 Nov 2017 20:35:52 +0100 Subject: [PATCH 1/3] Update calibration code using MIPs and IsoTracks --- .../HcalCalibAlgos/macros/AnalyzeLepTree.C | 786 ++++++++++++++++++ 1 file changed, 786 insertions(+) create mode 100644 Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C diff --git a/Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C b/Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C new file mode 100644 index 0000000000000..e932a45d816de --- /dev/null +++ b/Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C @@ -0,0 +1,786 @@ +////////////////////////////////////////////////////////////////////////////// +// This class analyzes the "Lep Tree" created by HBHEOfflineAnalyzer +// It has two constructors using either a pointer to the tree chain or +// the file name where the tree resides +// There are 2 additional arguments: +// mode (a packed integer) with bits specifying some action +// Bit 0 : 0 will produce set of histograms of energy deposit, +// active length, active length corrected energy; +// 1 will produce plots of p, nvx and scatter plots for +// each ieta +// 1 : 0 ignores the information on number of vertex +// 1 groups ranges of # of vertex 0:15:20:25:30:100 +// 2 : 0 ignores the information on track momentum +// 1 separate plots for certain momentum range +// (the range depends on ieta) +// 3 : 0 separate plot for each depth +// 1 sums up all depths +// 4-5: 0 no check on iphi +// 1 separate plot for each iphi +// 2 separate plot for each RBX +// modeLHC (integer) specifies the detector condition +// 0 Run1 detector (till 2016) +// 1 Plan36 detector +// 2 Phase1 detector +// 3 Plan1 detector +// 4 Phase2 detector +// +// AnalyzeLepTree a1(tree, mode, modeLHC); +// AnalyzeLepTree a1(fname, mode, modeLHC); +// a1.Loop(); +// a1.writeHisto(outfile); +// a1.plotHisto(drawStatBox,type,save); +// +// outfile (string) Name of the file where histograms to be written +// drawStatBox (bool) If Statbox to be drawn or not +// type (int) Each bit says what plots to be drawn +// If bit 0 of "mode" is 1 +// (0: momentum for each ieta; +// 1: number of vertex for each ieta; +// 2: 2D plot for nvx vs p for each ieta; +// 3: number of vertex for each ieta & depth; +// 4: momentum for each ieta & depth) +// If bit 0 of "mode" is 0 +// plots for all or specific depth & phi if +// "depth" as of (type/16)&15 is 0 or the specified +// value and "phi" as of (type/256)&127 is 0 or +// the specified value. Bit 0 set plots the energy +// distribution; 1 set plots active length corrected +// energy; 2 set plots charge; 3 set plots active +// length corrected charge distributions +// save (bool) If plots to be saved as pdf file or not +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +class AnalyzeLepTree { +public : + AnalyzeLepTree(TTree *tree, int mode=0, int modeLHC=3); + AnalyzeLepTree(std::string fname, int mode=0, int modeLHC=3); + virtual ~AnalyzeLepTree(); + virtual Int_t Cut(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry); + virtual Long64_t LoadTree(Long64_t entry); + virtual void Init(TTree *tree); + virtual void Loop(); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); + void writeHisto(std::string outfile); + std::vector plotHisto(bool drawStatBox, int type, bool save=false); +private: + void bookHisto(); + void plotHisto(std::map hists, + std::vector& cvs, bool save); + void plotHisto(std::map hists, int phi0, + int depth0, std::vector& cvs, + bool save); + TCanvas* plotHisto(TH1D* hist); + void plot2DHisto(std::map hists, + std::vector& cvs, bool save); + int getRBX(int eta); + int getPBin(int eta); + int getVxBin(); + int getDepthBin(int depth); + int getPhiBin(int eta); + int nDepthBins(int eta, int phi); + int nPhiBins(int eta); + int nPBins(int eta); + int nVxBins(); + unsigned int packID(int zside, int eta, int phi, int depth, + int nvx, int ipbin); + void unpackID(unsigned int id, int& zside, int& eta, + int& phi, int& depth, int& nvx, int& ipbin); +private: + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain + + // Declaration of leaf types + Int_t t_ieta; + Int_t t_iphi; + Int_t t_nvtx; + Double_t t_p; + std::vector *t_ene; + std::vector *t_charge; + std::vector *t_actln; + std::vector *t_depth; + + // List of branches + TBranch *b_t_ieta; //! + TBranch *b_t_iphi; //! + TBranch *b_t_nvtx; //! + TBranch *b_t_p; //! + TBranch *b_t_ene; //! + TBranch *b_t_charge; //! + TBranch *b_t_actln; //! + TBranch *b_t_depth; //! + + int mode_, modeLHC_; + std::vector npvbin_, iprange_; + std::vector prange_[5]; + std::map h_p_, h_nv_; + std::map h_pnv_; + std::map h_p2_, h_nv2_; + std::map h_Energy_, h_Ecorr_, h_Charge_, h_Chcorr_; +}; + +AnalyzeLepTree::AnalyzeLepTree(TTree *tree, int mode1, + int mode2) : mode_(mode1), modeLHC_(mode2) { + Init(tree); +} + +AnalyzeLepTree::AnalyzeLepTree(std::string fname, int mode1, + int mode2) : mode_(mode1), modeLHC_(mode2) { + new TFile(fname.c_str()); + TTree *tree = (TTree*)gDirectory->Get("Lep_Tree"); + Init(tree); +} + +AnalyzeLepTree::~AnalyzeLepTree() { + if (!fChain) return; + delete fChain->GetCurrentFile(); +} + +Int_t AnalyzeLepTree::GetEntry(Long64_t entry) { + // Read contents of entry. + if (!fChain) return 0; + return fChain->GetEntry(entry); +} + +Long64_t AnalyzeLepTree::LoadTree(Long64_t entry) { + // Set the environment to read one entry + if (!fChain) return -5; + Long64_t centry = fChain->LoadTree(entry); + if (centry < 0) return centry; + if (!fChain->InheritsFrom(TChain::Class())) return centry; + TChain *chain = (TChain*)fChain; + if (chain->GetTreeNumber() != fCurrent) { + fCurrent = chain->GetTreeNumber(); + Notify(); + } + return centry; +} + +void AnalyzeLepTree::Init(TTree *tree) { + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set object pointer + t_ene = 0; + t_charge = 0; + t_actln = 0; + t_depth = 0; + fChain = tree; + // Set branch addresses and branch pointers + if (!tree) return; + fCurrent = -1; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta); + fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi); + fChain->SetBranchAddress("t_nvtx", &t_nvtx, &b_t_nvtx); + fChain->SetBranchAddress("t_p", &t_p, &b_t_p); + fChain->SetBranchAddress("t_ene", &t_ene, &b_t_ene); + fChain->SetBranchAddress("t_charge", &t_charge, &b_t_charge); + fChain->SetBranchAddress("t_actln", &t_actln, &b_t_actln); + fChain->SetBranchAddress("t_depth", &t_depth, &b_t_depth); + Notify(); + + bookHisto(); +} + +Bool_t AnalyzeLepTree::Notify() { + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + return kTRUE; +} + +void AnalyzeLepTree::Show(Long64_t entry) { + // Print contents of entry. + // If entry is not specified, print current entry + if (!fChain) return; + fChain->Show(entry); +} + +Int_t AnalyzeLepTree::Cut(Long64_t ) { + // This function may be called from Loop. + // returns 1 if entry is accepted. + // returns -1 otherwise. + return 1; +} + +void AnalyzeLepTree::Loop() { + // In a ROOT session, you can do: + // Root > .L AnalyzeLepTree.C + // Root > AnalyzeLepTree t + // Root > t.GetEntry(12); // Fill t data members with entry number 12 + // Root > t.Show(); // Show values of entry 12 + // Root > t.Show(16); // Read and show values of entry 16 + // Root > t.Loop(); // Loop on all entries + // + + // This is the loop skeleton where: + // jentry is the global entry number in the chain + // ientry is the entry number in the current Tree + // Note that the argument to GetEntry must be: + // jentry for TChain::GetEntry + // ientry for TTree::GetEntry and TBranch::GetEntry + // + // To read only selected branches, Insert statements like: + // METHOD1: + // fChain->SetBranchStatus("*",0); // disable all branches + // fChain->SetBranchStatus("branchname",1); // activate branchname + // METHOD2: replace line + // fChain->GetEntry(jentry); //read all branches + //by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + Long64_t nentries = fChain->GetEntriesFast(); + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + // if (Cut(ientry) < 0) continue; + int zside = (t_ieta > 0) ? 1 : -1; + int eta = (t_ieta > 0) ? t_ieta : -t_ieta; + int phi = getPhiBin(eta); + int pbin = getPBin(eta); + int vbin = getVxBin(); + if ((mode_/1)%2 == 1) { + unsigned int id0 = packID(zside,eta,1,1,1,1); + std::map::iterator it1 = h_p_.find(id0); + if (it1 != h_p_.end()) (it1->second)->Fill(t_p); + std::map::iterator it2 = h_nv_.find(id0); + if (it2 != h_nv_.end()) (it2->second)->Fill(t_nvtx); + std::map::iterator it3 = h_pnv_.find(id0); + if (it3 != h_pnv_.end()) (it3->second)->Fill(t_nvtx,t_p); + unsigned int id1 = packID(zside,eta,1,1,1,pbin); + std::map::iterator it4 = h_p2_.find(id1); + if (it4 != h_p2_.end()) (it4->second)->Fill(t_p); + unsigned int id2 = packID(zside,eta,1,1,vbin,1); + std::map::iterator it5 = h_nv2_.find(id2); + if (it5 != h_nv2_.end()) (it5->second)->Fill(t_nvtx); + } else { + if ((mode_/8)%2 == 0) { + for (unsigned int k=0; ksize(); ++k) { + int depth = (*t_depth)[k]; + unsigned int id = packID(zside,eta,phi,depth+1,vbin,pbin); + double ene = (*t_ene)[k]; + double charge = (*t_charge)[k]; + double actl = (*t_actln)[k]; + if (ene > 0 && actl > 0 && charge > 0) { + std::map::iterator it1 = h_Energy_.find(id); + if (it1 != h_Energy_.end()) (it1->second)->Fill(ene); + std::map::iterator it2 = h_Ecorr_.find(id); + if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl); + std::map::iterator it3 = h_Charge_.find(id); + if (it3 != h_Charge_.end()) (it3->second)->Fill(charge); + std::map::iterator it4 = h_Chcorr_.find(id); + if (it4 != h_Chcorr_.end()) (it4->second)->Fill(charge/actl); + /* + if ((eta>20 && (t_iphi > 35)) || (t_iphi > 71)) std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth+1 << ":" << vbin << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_Charge_.end()) << ":" << (it4 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl; + if ((it1 == h_Energy_.end()) || (it2 == h_Ecorr_.end()) || (it3 == h_Charge_.end()) || (it4 == h_Chcorr_.end())) std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth+1 << ":" << vbin << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_Charge_.end()) << ":" << (it4 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl; + */ + } + } + } else { + double ene(0), actl(0), charge(0); + unsigned int id = packID(zside,eta,phi,1,vbin,pbin); + for (unsigned int k=0; ksize(); ++k) { + if ((*t_ene)[k] > 0 && (*t_actln)[k] > 0 && (*t_charge)[k] > 0) { + ene += (*t_ene)[k]; + charge += (*t_charge)[k]; + actl += (*t_actln)[k]; + } + } + std::map::iterator it1 = h_Energy_.find(id); + if (it1 != h_Energy_.end()) (it1->second)->Fill(ene); + std::map::iterator it2 = h_Ecorr_.find(id); + if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl); + std::map::iterator it3 = h_Charge_.find(id); + if (it3 != h_Charge_.end()) (it3->second)->Fill(charge); + std::map::iterator it4 = h_Chcorr_.find(id); + if (it4 != h_Chcorr_.end()) (it4->second)->Fill(charge/actl); + } + } + } +} + +void AnalyzeLepTree::bookHisto() { + + int npvbin[6] = {0, 15, 20, 25, 30, 100}; + npvbin_.clear(); + for (int i=0; i<5; ++i) prange_[i].clear(); + for (int i=0; i<6; ++i) npvbin_.push_back(npvbin[i]); + int ipbrng[30]= {0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,2,2,2,2,2,3,3,3,4,4,4,4}; + for (int i=0; i<30; ++i) iprange_.push_back(ipbrng[i]); + double prange0[7] = {0,30,45,55,75,100,500}; + double prange1[7] = {0,50,75,100,125,150,500}; + double prange2[7] = {0,60,75,100,125,150,500}; + double prange3[7] = {0,100,125,150,175,200,500}; + double prange4[7] = {0,125,150,175,200,250,500}; + for (int i=0; i<7; ++i) { + prange_[0].push_back(prange0[i]); + prange_[1].push_back(prange1[i]); + prange_[2].push_back(prange2[i]); + prange_[3].push_back(prange3[i]); + prange_[4].push_back(prange4[i]); + } + + char name[100], title[200]; + if ((mode_/1)%2 == 1) { + h_p_.clear(); h_nv_.clear(); h_pnv_.clear(); h_nv2_.clear(); h_p2_.clear(); + for (int ieta=-26; ieta<=26; ++ieta) { + if (ieta != 0) { + int zside = (ieta>0) ? 1 : -1; + int eta = (ieta>0) ? ieta : -ieta; + unsigned int id0 = packID(zside,eta,1,1,1,1); + sprintf(name, "peta%d", ieta); + sprintf(title,"Momentum for i#eta = %d (GeV)", ieta); + h_p_[id0] = new TH1D(name, title, 100, 0.0, 500.0); + sprintf(name, "nveta%d", ieta); + sprintf(title,"Number of Vertex for i#eta = %d", ieta); + h_nv_[id0] = new TH1D(name, title, 100, 0.0, 100.0); + sprintf(name, "pnveta%d", ieta); + sprintf(title,"i#eta = %d", ieta); + TH2D* h2 = new TH2D(name, title, 100, 0.0, 100.0, 100, 0.0, 500.0); + h2->GetXaxis()->SetTitle("Number of Vertex"); + h2->GetYaxis()->SetTitle("Momentum (GeV)"); + h_pnv_[id0] = h2; + char etas[10]; + sprintf (etas, "i#eta=%d", ieta); + char name[100], title[200]; + for (int pbin=0; pbin= 0 && eta < (int)(iprange_.size())) ? + iprange_[eta]-1 : iprange_[0]; + sprintf (ps, "p=%d:%d", (int)prange_[np][pbin], (int)prange_[np][pbin+1]); + }; + unsigned int id = packID(zside,eta,1,1,1,pbin); + sprintf (name,"pvx%d111%d",ieta,pbin); + sprintf (title,"Momentum for %s %s", etas, ps); + h_p2_[id] = new TH1D(name,title,100,0.0,500.0); + h_p2_[id]->Sumw2(); + } + for (int vbin=0; vbinSumw2(); + } + } + } + } else { + h_Energy_.clear(); h_Ecorr_.clear(); h_Charge_.clear(); h_Chcorr_.clear(); + for (int ieta=-26; ieta<=26; ++ieta) { + if (ieta != 0) { + int zside = (ieta>0) ? 1 : -1; + int eta = (ieta>0) ? ieta : -ieta; + char etas[20]; + sprintf (etas, "i#eta=%d", ieta); + for (int iphi=0; iphi= 0 && eta < (int)(iprange_.size())) ? + iprange_[eta]-1 : iprange_[0]; + sprintf (ps, "p=%d:%d", (int)prange_[np][pbin], (int)prange_[np][pbin+1]); + } else { + sprintf (ps, "all p"); + }; + for (int vbin=0; vbinSumw2(); + sprintf (name,"EcorE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin); + sprintf (title,"Active length corrected energy for %s %s %s %s %s (GeV/cm)", etas, phis, deps, ps, vtx); + h_Ecorr_[id] = new TH1D(name,title,4000,0.0,10.0); + h_Ecorr_[id]->Sumw2(); + sprintf (name,"ChrgE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin); + sprintf (title,"Measured charge for %s %s %s %s %s (cm)", etas, phis, deps, ps, vtx); + h_Charge_[id] = new TH1D(name,title,20,0.0,20.0); + h_Charge_[id]->Sumw2(); + sprintf (name,"ChcorE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin); + sprintf (title,"Active length corrected charge for %s %s %s %s %s (cm)", etas, phis, deps, ps, vtx); + h_Chcorr_[id] = new TH1D(name,title,20,0.0,20.0); + h_Chcorr_[id]->Sumw2(); + } + } + } + } + } + } + } +} + +void AnalyzeLepTree::writeHisto(std::string fname) { + + TFile* output_file = TFile::Open(fname.c_str(),"RECREATE"); + output_file->cd(); + if ((mode_/1)%2 == 1) { + for (std::map::const_iterator itr = h_p_.begin(); + itr != h_p_.end(); ++itr) (itr->second)->Write(); + for (std::map::const_iterator itr = h_nv_.begin(); + itr != h_nv_.end(); ++itr) (itr->second)->Write(); + for (std::map::const_iterator itr = h_pnv_.begin(); + itr != h_pnv_.end(); ++itr) (itr->second)->Write(); + for (std::map::const_iterator itr = h_p2_.begin(); + itr != h_p2_.end(); ++itr) (itr->second)->Write(); + for (std::map::const_iterator itr = h_nv2_.begin(); + itr != h_nv2_.end(); ++itr) (itr->second)->Write(); + } else { + for (int ieta=-26; ieta<=26; ++ieta) { + if (ieta != 0) { + char dirname[50]; + sprintf (dirname, "DieMuonEta%d", ieta); + TDirectory* d_output = output_file->mkdir(dirname); + d_output->cd(); + int zside = (ieta>0) ? 1 : -1; + int eta = (ieta>0) ? ieta : -ieta; + for (int iphi=0; iphi::const_iterator itr; + itr = h_Energy_.find(id); + if (itr != h_Energy_.end()) (itr->second)->Write(); + itr = h_Ecorr_.find(id); + if (itr != h_Ecorr_.end()) (itr->second)->Write(); + itr = h_Charge_.find(id); + if (itr != h_Charge_.end()) (itr->second)->Write(); + itr = h_Chcorr_.find(id); + if (itr != h_Chcorr_.end()) (itr->second)->Write(); + } + } + } + } + } + output_file->cd(); + } + } +} + +std::vector AnalyzeLepTree::plotHisto(bool drawStatBox, int type, + bool save) { + + std::vector cvs; + gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(0); + if (drawStatBox) { + gStyle->SetOptStat(111110); gStyle->SetOptFit(1); + } else { + gStyle->SetOptStat(0); gStyle->SetOptFit(0); + } + + if ((mode_/1)%2 == 1) { + if (type%2 > 0) plotHisto(h_p_, cvs, save); + if ((type/2)%2 > 0) plotHisto(h_nv_, cvs, save); + if ((type/4)%2 > 0) plot2DHisto(h_pnv_, cvs, save); + if ((type/8)%2 > 0) plotHisto(h_nv2_, cvs, save); + if ((type/16)%2 > 0) plotHisto(h_p2_, cvs, save); + } else { + int depth0 = (type/16)&15; + int phi0 = (type/256)&127; + bool doEn = ((type/1)%2 > 0); + bool doEnL = ((type/2)%2 > 0); + bool doChg = ((type/4)%2 > 0); + bool doChL = ((type/8)%2 > 0); + if (doEn) plotHisto(h_Energy_, phi0, depth0, cvs, save); + if (doEnL) plotHisto(h_Ecorr_, phi0, depth0, cvs, save); + if (doChg) plotHisto(h_Charge_, phi0, depth0, cvs, save); + if (doChL) plotHisto(h_Chcorr_, phi0, depth0, cvs, save); + } + return cvs; +} + +void AnalyzeLepTree::plotHisto(std::map hists, + std::vector& cvs, bool save) { + + for (std::map::const_iterator itr = hists.begin(); + itr != hists.end(); ++itr) { + TH1D* hist = (itr->second); + if (hist != 0) { + TCanvas* pad = plotHisto(hist); + cvs.push_back(pad); + if (save) { + char name[100]; + sprintf (name, "c_%s.pdf", pad->GetName()); + pad->Print(name); + } + } + } +} + +void AnalyzeLepTree::plotHisto(std::map hists, int phi0, + int depth0, std::vector& cvs, + bool save) { + + for (std::map::const_iterator itr=hists.begin(); + itr != hists.end(); ++itr) { + int zside, eta, phi, depth, pbin, vbin; + unpackID(itr->first,zside,eta,phi,depth,vbin,pbin); + TH1D* hist = itr->second; + if (((depth == depth0) || (depth0 == 0)) && + ((phi == phi0) || (phi0 == 0)) && (hist != 0)) { + TCanvas* pad = plotHisto(hist); + cvs.push_back(pad); + if (save) { + char name[100]; + sprintf (name, "c_%s.pdf", pad->GetName()); + pad->Print(name); + } + } + } +} + +TCanvas* AnalyzeLepTree::plotHisto(TH1D* hist) { + + TCanvas *pad = new TCanvas(hist->GetName(), hist->GetName(), 700, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + hist->GetXaxis()->SetTitleSize(0.04); + hist->GetXaxis()->SetTitle(hist->GetTitle()); + hist->GetYaxis()->SetTitle("Tracks"); + hist->GetYaxis()->SetLabelOffset(0.005); + hist->GetYaxis()->SetTitleSize(0.04); + hist->GetYaxis()->SetLabelSize(0.035); + hist->GetYaxis()->SetTitleOffset(1.30); + hist->SetMarkerStyle(20); + hist->SetLineWidth(2); + hist->Draw(); + pad->Update(); + TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats"); + if (st1 != NULL) { + st1->SetY1NDC(0.70); st1->SetY2NDC(0.90); + st1->SetX1NDC(0.65); st1->SetX2NDC(0.90); + } + pad->Modified(); + pad->Update(); + return pad; +} + +void AnalyzeLepTree::plot2DHisto(std::map hists, + std::vector& cvs, bool save) { + char name[100]; + for (std::map::const_iterator itr = hists.begin(); + itr != hists.end(); ++itr) { + TH2D* hist = (itr->second); + if (hist != 0) { + TCanvas *pad = new TCanvas(hist->GetName(), hist->GetName(), 700, 700); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + hist->GetXaxis()->SetTitleSize(0.04); + hist->GetYaxis()->SetLabelOffset(0.005); + hist->GetYaxis()->SetTitleSize(0.04); + hist->GetYaxis()->SetLabelSize(0.035); + hist->GetYaxis()->SetTitleOffset(1.30); + hist->Draw("COLZ"); + pad->Update(); + TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats"); + if (st1 != NULL) { + st1->SetY1NDC(0.65); st1->SetY2NDC(0.90); + st1->SetX1NDC(0.65); st1->SetX2NDC(0.90); + } + pad->Modified(); + pad->Update(); + cvs.push_back(pad); + if (save) { + sprintf (name, "c_%s.pdf", pad->GetName()); + pad->Print(name); + } + } + } +} + +int AnalyzeLepTree::getRBX(int eta) { + int rbx(1); + int phi = (eta > 20) ? (2*t_iphi + 1) : (t_iphi + 1); + if (phi > 2 && phi < 71) rbx = (phi+5)/4; + return rbx; +} + +int AnalyzeLepTree::getPBin(int eta) { + int bin(0); + if ((mode_/4)%2 == 1) { + int np = (eta >= 0 && eta < (int)(iprange_.size())) ? iprange_[eta] : iprange_[0]; + for (unsigned int k=1; k 20) ? (2*t_iphi + 1) : (t_iphi + 1); + } else if ((mode_/16)%4 == 2) { + bin = getRBX(eta); + } + return bin; +} + +int AnalyzeLepTree::nDepthBins(int eta, int phi) { + // Run 1 scenario + int nDepthR1[29]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,3,1,2,2,2,2,2,2,2,2,2,3,3,2}; + // Run 2 scenario from 2018 + int nDepthR2[29]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,4,3,5,6,6,6,6,6,6,6,7,7,7,3}; + // Run 3 scenario + int nDepthR3[29]={4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,5,6,6,6,6,6,6,6,7,7,7,3}; + // Run 4 scenario + int nDepthR4[29]={4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,7,7,7,7,7,7,7,7,7,7,7,7,7}; + // modeLHC = 0 --> corresponds to Run 1 (valid till 2016) + // = 1 --> corresponds to Run 2 (2018 geometry) + // = 2 --> corresponds to Run 3 (post LS2) + // = 3 --> corresponds to 2017 (Plan 1) + // = 4 --> corresponds to Run4 (post LS3) + int nbin(0); + if (modeLHC_ == 0) { + nbin = nDepthR1[eta-1]; + } else if (modeLHC_ == 1) { + nbin = nDepthR2[eta-1]; + } else if (modeLHC_ == 2) { + nbin = nDepthR3[eta-1]; + } else if (modeLHC_ == 3) { + if (phi > 0) { + if (eta >= 16 && phi >= 63 && phi <= 66) { + nbin = nDepthR2[eta-1]; + } else { + nbin = nDepthR1[eta-1]; + } + } else { + if (eta >= 16) { + nbin = (nDepthR2[eta-1] > nDepthR1[eta-1]) ? nDepthR2[eta-1] : nDepthR1[eta-1]; + } else { + nbin = nDepthR1[eta-1]; + } + } + } else { + if (eta > 0 && eta < 30) { + nbin = nDepthR4[eta-1]; + } else { + nbin = nDepthR4[28]; + } + } + return nbin; +} + +int AnalyzeLepTree::nPhiBins(int eta) { + int nphi = (eta <= 20) ? 72 : 36; + if (modeLHC_ == 4 && eta > 16) nphi = 360; + if ((mode_/16)%4 == 0) nphi = 1; + else if ((mode_/16)%4 == 2) nphi = 18; + return nphi; +} + +int AnalyzeLepTree::nPBins(int eta) { + int bin(1); + if ((mode_/4)%2 == 1) { + int np = (eta >= 0 && eta < (int)(iprange_.size())) ? iprange_[eta]-1 : iprange_[0]; + bin = (int)(prange_[np].size())-1; + } + return bin; +} + +int AnalyzeLepTree::nVxBins() { + int nvx = ((mode_/2)%2 == 1) ? ((int)(npvbin_.size())-1) : 1; + return nvx; +} + +unsigned int AnalyzeLepTree::packID(int zside, int eta, int phi, int depth, + int nvx, int ipbin) { + unsigned int id = (zside > 0) ? 1 : 0; + id |= (((nvx&7)<<19) | ((ipbin&7)<<16) | ((depth&7)<<13) | ((eta&31)<<8) | + ((phi&127)<<1)); + return id; +} + +void AnalyzeLepTree::unpackID(unsigned int id, int& zside, int& eta, int& phi, + int& depth, int& nvx, int& ipbin) { + zside = (id%2 == 0) ? -1 : 1; + phi = (id >> 1) & 127; + eta = (id >> 8) & 31; + depth = (id >> 13) & 7; + ipbin = (id >> 16) & 7; + nvx = (id >> 19) & 7; +} From 3ed91f1aa242209e7b36fe5f9f8a803398dd60dc Mon Sep 17 00:00:00 2001 From: Sunanda Date: Fri, 17 Nov 2017 20:35:55 +0100 Subject: [PATCH 2/3] Update calibration code using MIPs and IsoTracks --- .../HcalCalibAlgos/macros/CalibMonitor.C | 121 ++++++++-- .../macros/HBHEMuonOfflineAnalyzer.C | 211 +++++++++++------- .../HcalCalibAlgos/macros/OfflineMain.C | 32 ++- .../HcalCalibAlgos/test/HcalIsoTrkAnalyzer.cc | 17 +- 4 files changed, 253 insertions(+), 128 deletions(-) diff --git a/Calibration/HcalCalibAlgos/macros/CalibMonitor.C b/Calibration/HcalCalibAlgos/macros/CalibMonitor.C index c03c11449b9fc..6fc338e87550c 100644 --- a/Calibration/HcalCalibAlgos/macros/CalibMonitor.C +++ b/Calibration/HcalCalibAlgos/macros/CalibMonitor.C @@ -1196,6 +1196,7 @@ public : std::vector *t_hltbits; std::vector *t_ietaAll; std::vector *t_ietaGood; + std::vector *t_trackType; // List of branches TBranch *b_t_Tracks; //! @@ -1209,8 +1210,10 @@ public : TBranch *b_t_hltbits; //! TBranch *b_t_ietaAll; //! TBranch *b_t_ietaGood; //! + TBranch *b_t_trackType; //! - GetEntries(const std::string& fname, const std::string& dirname); + GetEntries(const std::string& fname, const std::string& dirname, + const unsigned int bit1, const unsigned int bit2); virtual ~GetEntries(); virtual Int_t Cut(Long64_t entry); virtual Int_t GetEntry(Long64_t entry); @@ -1221,17 +1224,20 @@ public : virtual void Show(Long64_t entry = -1); private: - TH1I *h_tk[3], *h_eta[2]; - TH1D *h_eff; + unsigned int bit_[2]; + TH1I *h_tk[3], *h_eta[4]; + TH1D *h_eff[3]; }; -GetEntries::GetEntries(const std::string& fname, const std::string& dirnm) { +GetEntries::GetEntries(const std::string& fname, const std::string& dirnm, + const unsigned int bit1, const unsigned int bit2) { TFile *file = new TFile(fname.c_str()); TDirectory *dir = (TDirectory*)file->FindObjectAny(dirnm.c_str()); std::cout << fname << " file " << file << " " << dirnm << " " << dir << std::endl; TTree *tree = (TTree*)dir->Get("EventInfo"); std::cout << "CalibTree " << tree << std::endl; + bit_[0] = bit1; bit_[1] = bit2; Init(tree); } @@ -1274,6 +1280,7 @@ void GetEntries::Init(TTree *tree) { t_hltbits = 0; t_ietaAll = 0; t_ietaGood = 0; + t_trackType = 0; t_L1Bit = false; fChain = tree; fCurrent = -1; @@ -1290,6 +1297,7 @@ void GetEntries::Init(TTree *tree) { fChain->SetBranchAddress("t_hltbits", &t_hltbits, &b_t_hltbits); fChain->SetBranchAddress("t_ietaAll", &t_ietaAll, &b_t_ietaAll); fChain->SetBranchAddress("t_ietaGood", &t_ietaGood, &b_t_ietaGood); + fChain->SetBranchAddress("t_trackType", &t_trackType, &b_t_trackType); Notify(); h_tk[0] = new TH1I("Track0", "# of tracks produced", 2000, 0, 2000); @@ -1297,7 +1305,11 @@ void GetEntries::Init(TTree *tree) { h_tk[2] = new TH1I("Track2", "# of tracks saved in tree", 2000, 0, 2000); h_eta[0] = new TH1I("Eta0", "i#eta (All Tracks)", 60, -30, 30); h_eta[1] = new TH1I("Eta1", "i#eta (Good Tracks)", 60, -30, 30); - h_eff = new TH1D("Eta2", "i#eta (Selection Efficiency)", 60, -30, 30); + h_eta[2] = new TH1I("Eta2", "i#eta (Loose Isolated Tracks)",60, -30, 30); + h_eta[3] = new TH1I("Eta3", "i#eta (Tight Isolated Tracks)",60, -30, 30); + h_eff[0] = new TH1D("Eff0", "i#eta (Selection Efficiency)", 60, -30, 30); + h_eff[1] = new TH1D("Eff1", "i#eta (Loose Isolation Efficiency)", 60,-30,30); + h_eff[2] = new TH1D("Eta2", "i#eta (Tight Isolation Efficiency)", 60,-30,30); } Bool_t GetEntries::Notify() { @@ -1354,6 +1366,9 @@ void GetEntries::Loop() { Long64_t nbytes = 0, nb = 0; int kount(0), selected(0); int l1(0), hlt(0), loose(0), tight(0); + int allHLT[3] = {0,0,0}; + int looseHLT[3] = {0,0,0}; + int tightHLT[3] = {0, 0, 0}; for (Long64_t jentry=0; jentrysize(); ++k) { + if ((*t_hltbits)[k] > 0) { + pass = true; + for (int i=0; i<2; ++i) if (k == bit_[i]) passt[i] = true; + } + } + if (pass) { + ++allHLT[0]; + for (int i=0; i<2; ++i) if (passt[i]) ++allHLT[i+1]; + if (t_TracksLoose > 0) { + ++looseHLT[0]; + for (int i=0; i<2; ++i) if (passt[i]) ++looseHLT[i+1]; + } + if (t_TracksTight > 0) { + ++tightHLT[0]; + for (int i=0; i<2; ++i) if (passt[i]) ++tightHLT[i+1]; + } + } for (unsigned int k=0; ksize(); ++k) h_eta[0]->Fill((*t_ietaAll)[k]); - for (unsigned int k=0; ksize(); ++k) + for (unsigned int k=0; ksize(); ++k) { h_eta[1]->Fill((*t_ietaGood)[k]); + if (t_trackType->size() > 0) { + if ((*t_trackType)[k] > 0) h_eta[2]->Fill((*t_ietaGood)[k]); + if ((*t_trackType)[k] > 1) h_eta[3]->Fill((*t_ietaGood)[k]); + } + } } double ymaxk(0); - for (int i=1; i<=h_eff->GetNbinsX(); ++i) { + for (int i=1; i<=h_eff[0]->GetNbinsX(); ++i) { double rat(0), drat(0); if (h_eta[0]->GetBinContent(i) > ymaxk) ymaxk = h_eta[0]->GetBinContent(i); if ((h_eta[1]->GetBinContent(i) > 0)&&(h_eta[0]->GetBinContent(i) > 0)) { @@ -1386,8 +1425,26 @@ void GetEntries::Loop() { drat= rat*std::sqrt(pow((h_eta[1]->GetBinError(i)/h_eta[1]->GetBinContent(i)),2) + pow((h_eta[0]->GetBinError(i)/h_eta[0]->GetBinContent(i)),2)); } - h_eff->SetBinContent(i,rat); - h_eff->SetBinError(i,drat); + h_eff[0]->SetBinContent(i,rat); + h_eff[0]->SetBinError(i,drat); + if ((h_eta[1]->GetBinContent(i) > 0)&&(h_eta[2]->GetBinContent(i) > 0)) { + rat = h_eta[2]->GetBinContent(i)/h_eta[1]->GetBinContent(i); + drat= rat*std::sqrt(pow((h_eta[2]->GetBinError(i)/h_eta[2]->GetBinContent(i)),2) + + pow((h_eta[1]->GetBinError(i)/h_eta[1]->GetBinContent(i)),2)); + } else { + rat = drat = 0; + } + h_eff[1]->SetBinContent(i,rat); + h_eff[1]->SetBinError(i,drat); + if ((h_eta[1]->GetBinContent(i) > 0)&&(h_eta[3]->GetBinContent(i) > 0)) { + rat = h_eta[3]->GetBinContent(i)/h_eta[1]->GetBinContent(i); + drat= rat*std::sqrt(pow((h_eta[3]->GetBinError(i)/h_eta[3]->GetBinContent(i)),2) + + pow((h_eta[1]->GetBinError(i)/h_eta[1]->GetBinContent(i)),2)); + } else { + rat = drat = 0; + } + h_eff[1]->SetBinContent(i,rat); + h_eff[1]->SetBinError(i,drat); } std::cout << "===== " << kount << " events passed trigger of which " << selected << " events get selected =====\n" << std::endl; @@ -1395,11 +1452,20 @@ void GetEntries::Loop() { << " events passed HLT and " << loose << ":" << tight << " events have at least 1 track candidate with loose:tight" << " isolation cut =====\n" << std::endl; + for (int i=0; i<3; ++i) { + char tbit[20]; + if (i == 0) sprintf (tbit, "Any"); + else sprintf (tbit, "%3d", bit_[i-1]); + std::cout << "Satisfying HLT trigger bit " << tbit << " Kount " + << allHLT[i] << " Loose " << looseHLT[i] << " Tight " + << tightHLT[i] << std::endl; + } + gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); gStyle->SetOptStat(1110); gStyle->SetOptTitle(0); - int color[3] = {kBlack, kRed, kBlue}; - int lines[3] = {1, 2, 3}; + int color[5] = {kBlack, kRed, kBlue, kMagenta, kCyan}; + int lines[5] = {1, 2, 3, 4, 5}; TCanvas *pad1 = new TCanvas("c_track", "c_track", 500, 500); pad1->SetRightMargin(0.10); pad1->SetTopMargin(0.10); @@ -1464,12 +1530,13 @@ void GetEntries::Loop() { pad2->SetTopMargin(0.10); pad2->SetFillColor(kWhite); pad2->SetLogy(); - std::string titl2[2] = {"All Tracks", "Selected Tracks"}; - TLegend *legend2 = new TLegend(0.11, 0.82, 0.50, 0.89); + std::string titl2[4] = {"All Tracks", "Selected Tracks", "Loose Isolation", + "Tight Isolation"}; + TLegend *legend2 = new TLegend(0.11, 0.75, 0.50, 0.89); legend2->SetFillColor(kWhite); i2 = (int)(0.001*ymaxk) + 1; ymax = 1000.0*i2; - for (int k=0; k<2; ++k) { + for (int k=0; k<4; ++k) { h_eta[k]->GetYaxis()->SetRangeUser(1,ymax); h_eta[k]->SetLineColor(color[k]); h_eta[k]->SetMarkerColor(color[k]); @@ -1485,7 +1552,7 @@ void GetEntries::Loop() { pad2->Update(); ymax = 0.90; //double ymin = 0.10; - for (int k=0; k<2; ++k) { + for (int k=0; k<4; ++k) { TPaveStats* st1 = (TPaveStats*)h_eta[k]->GetListOfFunctions()->FindObject("stats"); if (st1 != NULL) { st1->SetLineColor(color[k]); @@ -1500,18 +1567,28 @@ void GetEntries::Loop() { legend2->Draw("same"); pad2->Update(); + std::string titl3[3] = {"Selection", "Loose Isolation", "Tight Isolation"}; TCanvas *pad3 = new TCanvas("c_effi", "c_effi", 500, 500); + TLegend *legend3 = new TLegend(0.11, 0.785, 0.50, 0.89); pad3->SetRightMargin(0.10); pad3->SetTopMargin(0.10); pad3->SetFillColor(kWhite); pad3->SetLogy(); - h_eff->SetStats(0); - h_eff->SetMarkerStyle(20); - h_eff->GetXaxis()->SetTitle("i#eta"); - h_eff->GetYaxis()->SetTitle("Efficiency"); - h_eff->GetYaxis()->SetLabelOffset(0.005); - h_eff->GetYaxis()->SetTitleOffset(1.20); - h_eff->Draw(); + for (int k=0; k<3; ++k) { + h_eff[k]->SetStats(0); + h_eff[k]->SetMarkerStyle(20); + h_eff[k]->SetLineColor(color[k]); + h_eff[k]->SetMarkerColor(color[k]); + h_eff[k]->GetXaxis()->SetTitle("i#eta"); + h_eff[k]->GetYaxis()->SetTitle("Efficiency"); + h_eff[k]->GetYaxis()->SetLabelOffset(0.005); + h_eff[k]->GetYaxis()->SetTitleOffset(1.20); + if (k == 0) h_eff[k]->Draw(""); + else h_eff[k]->Draw("same"); + legend3->AddEntry(h_eff[k],titl3[k].c_str(),"l"); + } pad3->Modified(); pad3->Update(); + legend3->Draw("same"); + pad3->Update(); } diff --git a/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C b/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C index 58192569d20b2..82c1d33dbd5b7 100644 --- a/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C +++ b/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C @@ -1,3 +1,32 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// HBHEMuonOfflineAnalyzer h1(tree,outfile,flag,mode,maxDHB,maxDHE,runLo, +// runHi,etaMin,etaMax); +// HBHEMuonOfflineAnalyzer h1(infile,outfile,flag,mode,maxDHB,maxDHE,runLo, +// runHi,etaMin,etaMax); +// h1.Loop() +// +// tree TTree* Pointer to the tree chain +// infile const char* Name of the input file +// outfile const char* Name of the output file +// (dyll_PU20_25_output_10.root) +// flag int Flag of 2 digits: to where o decides if corrected +// (1) or default (0) energy to be used; t decides +// if all depths to be merged (1) or not (0) +// (default is 0) +// mode int Geometry file used 0:(defined by maxDHB/HE); +// 1 (Run 1; valid till 2016); 2 (Run 2; 2018); +// 3 (Run 3; post LS2); 4 (2017 Plan 1); +// 5 (Run 4; post LS3); default (4) +// maxDHB int Maximum number of depths for HB (4) +// maxDHE int Maximum number of depths for HE (7) +// runLO int Minimum run number (1) +// runHI int Maximum run number (99999999) +// etaMin int Minimum (absolute) eta value (1) +// etaMax int Maximum (absolute) eta value (29) +// +/////////////////////////////////////////////////////////////////////////////// + #include #include #include @@ -287,8 +316,8 @@ public : TBranch *b_all_triggers; //! - HBHEMuonOfflineAnalyzer(TTree *tree=0, const char *outfile="dyll_PU20_25_output_10.root", int flag=0, int mode=4, int maxDHB=5, int maxDHE=7, int runLo=297494, int runHi=99999999); - HBHEMuonOfflineAnalyzer(const char *infile, const char *outfile="dyll_PU20_25_output_10.root", int flag=0, int mode=4, int maxDHB=5, int maxDHE=7, int runLo=297494, int runHi=99999999); + HBHEMuonOfflineAnalyzer(TTree *tree=0, const char *outfile="dyll_PU20_25_output_10.root", int flag=0, int mode=4, int maxDHB=4, int maxDHE=7, int runLo=1, int runHi=99999999, int etaMin=1, int etaMax=29); + HBHEMuonOfflineAnalyzer(const char *infile, const char *outfile="dyll_PU20_25_output_10.root", int flag=0, int mode=4, int maxDHB=4, int maxDHE=7, int runLo=1, int runHi=99999999, int etaMin=1, int etaMax=29); // mode of LHC is kept 1 for 2017 scenario as no change in depth segmentation // mode of LHC is 0 for 2019 virtual ~HBHEMuonOfflineAnalyzer(); @@ -296,7 +325,7 @@ public : virtual Int_t GetEntry(Long64_t entry); virtual Long64_t LoadTree(Long64_t entry); virtual void Init(TTree *tree, int flag, int mode, int maxDHB, int maxDHE, - int runLo, int runHi); + int runLo, int runHi, int etaMin, int etaMax); virtual void Loop(); virtual Bool_t Notify(); virtual void Show(Long64_t entry = -1); @@ -324,7 +353,8 @@ private: //3x16x72x2 + 5x4x72x2 + 5x9x36x2 static const int maxHist=20000;//13032; static const int nCut_ = 1; - int modeLHC, maxDepthHB_, maxDepthHE_, maxDepth_, runLo_, runHi_; + int modeLHC_, maxDepthHB_, maxDepthHE_, maxDepth_, runLo_, runHi_; + int etaMin_, etaMax_; bool useCorrect_, mergeDepth_; int nHist, nDepths[maxEta], nDepthsPhi[maxEta],indxEta[maxEta][maxDep][maxPhi]; TFile *output_file; @@ -332,7 +362,7 @@ private: TTree *outtree_; int t_ieta, t_iphi, t_nvtx; double t_p; - std::vector t_ene, t_actl, t_encor; + std::vector t_ene, t_actln, t_charge; std::vector t_depth; TH1D *h_Pt_Muon[3], *h_Eta_Muon[3], *h_Phi_Muon[3], *h_P_Muon[3]; @@ -362,9 +392,10 @@ HBHEMuonOfflineAnalyzer::HBHEMuonOfflineAnalyzer(TTree *tree, const char* outFileName, int flag, int mode, int maxDHB, int maxDHE, - int runLo, int runHi) { + int runLo, int runHi, + int etaMin, int etaMax) { - Init(tree, flag, mode, maxDHB, maxDHE, runLo, runHi); + Init(tree, flag, mode, maxDHB, maxDHE, runLo, runHi, etaMin, etaMax); //Now book histograms BookHistograms(outFileName); @@ -374,14 +405,15 @@ HBHEMuonOfflineAnalyzer::HBHEMuonOfflineAnalyzer(const char* infile, const char* outFileName, int flag, int mode, int maxDHB, int maxDHE, - int runLo, int runHi) { + int runLo, int runHi, + int etaMin, int etaMax) { TFile *f = new TFile(infile); TDirectory *dir = (TDirectory*)f->Get("hcalHBHEMuon"); TTree *tree(0); dir->GetObject("TREE",tree); - Init(tree, flag, mode, maxDHB, maxDHE, runLo, runHi); + Init(tree, flag, mode, maxDHB, maxDHE, runLo, runHi, etaMin, etaMax); //Now book histograms BookHistograms(outFileName); @@ -418,14 +450,21 @@ Long64_t HBHEMuonOfflineAnalyzer::LoadTree(Long64_t entry) { } void HBHEMuonOfflineAnalyzer::Init(TTree *tree, int flag, int mode, int maxDHB, - int maxDHE, int runLo, int runHi) { + int maxDHE, int runLo, int runHi, int etaMin, + int etaMax) { - modeLHC = mode; + modeLHC_ = mode; maxDepthHB_ = maxDHB; maxDepthHE_ = maxDHE; maxDepth_ = (maxDepthHB_ > maxDepthHE_) ? maxDepthHB_ : maxDepthHE_; runLo_ = runLo; runHi_ = runHi; + etaMin_ = (etaMin > 0) ? etaMin : 1; + etaMax_ = (etaMax <= 29) ? etaMax : 29; + if (etaMax_ <= etaMin_) { + if (etaMax_ == 29) etaMin_ = etaMax_ - 1; + else etaMax_ = etaMin_ + 1; + } useCorrect_ = ((flag%10) > 0); mergeDepth_ = (((flag/10)%10) > 0); @@ -721,7 +760,7 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (debug_) std::cout << "Run " << Run_No << " Event " << Event_No << " Muons " << pt_of_muon->size() << std::endl; for (unsigned int ml = 0; ml< pt_of_muon->size(); ml++) { - t_ene.clear(); t_actl.clear(); t_encor.clear(); t_depth.clear(); + t_ene.clear(); t_charge.clear(); t_actln.clear(); t_depth.clear(); if(debug_) std::cout << "ecal_det_id " << ecal_detID->at(ml) << std::endl; @@ -760,7 +799,7 @@ void HBHEMuonOfflineAnalyzer::Loop() { else select = LooseMuon(ml); - if (select) { + if (select && ((eta+1) >= etaMin_) && ((eta+1) <= etaMax_)) { // h_P_Muon[cut]->Fill(p_of_muon->at(ml)); h_Pt_Muon[cut]->Fill(pt_of_muon->at(ml)); h_Eta_Muon[cut]->Fill(eta_of_muon->at(ml)); @@ -886,8 +925,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { h_charge_bg[cut][ind]->Fill(chargeBG); t_ene.push_back(en2); - t_actl.push_back(energyFill); - t_encor.push_back(en2/energyFill); + t_charge.push_back(chargeS); + t_actln.push_back(energyFill); t_depth.push_back(0); outtree_->Fill(); @@ -972,7 +1011,7 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (debug_) std::cout<<"enter ok1"<at(ml)==1) { - if(energyFill > 0) { + if (energyFill > 0) { h_Hot_MuonEnergy_hcal_HotCell[cut][ind]->Fill(en2); h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[cut][ind]->Fill(en2/energyFill); h_active_length_Fill[cut][ind]->Fill(energyFill); @@ -980,15 +1019,15 @@ void HBHEMuonOfflineAnalyzer::Loop() { h_charge_signal[cut][ind]->Fill(chargeS); h_charge_bg[cut][ind]->Fill(chargeBG); t_ene.push_back(en2); - t_actl.push_back(energyFill); - t_encor.push_back(en2/energyFill); + t_charge.push_back(chargeS); + t_actln.push_back(energyFill); // added depth vector AmanKalsi t_depth.push_back(dep); fillTree = true; } else { t_ene.push_back(-999.0); - t_actl.push_back(-999.0); - t_encor.push_back(-999.0); + t_charge.push_back(-999.0); + t_actln.push_back(-999.0); t_depth.push_back(-999.0); } if(debug_) std::cout<<"enter hot cell"<cd(); outtree_ = new TTree("Lep_Tree","Lep_Tree"); - outtree_->Branch("t_ieta", &t_ieta); - outtree_->Branch("t_iphi", &t_iphi); - outtree_->Branch("t_nvtx", &t_nvtx); - outtree_->Branch("t_p", &t_p); - outtree_->Branch("t_ene", &t_ene); - outtree_->Branch("t_actl", &t_actl); - outtree_->Branch("t_encor", &t_encor); - outtree_->Branch("t_depth", &t_depth); + outtree_->Branch("t_ieta", &t_ieta); + outtree_->Branch("t_iphi", &t_iphi); + outtree_->Branch("t_nvtx", &t_nvtx); + outtree_->Branch("t_p", &t_p); + outtree_->Branch("t_ene", &t_ene); + outtree_->Branch("t_charge", &t_charge); + outtree_->Branch("t_actln", &t_actln); + outtree_->Branch("t_depth", &t_depth); std::string type[]={"tight"};//,"soft","loose"}; char name[128], title[500]; nHist = 0; - for (int eta=0; eta<29; ++eta) { + for (int eta=etaMin_; eta<=etaMax_; ++eta) { - int nDepth = NDepthBins(eta+1,-1); - int nPhi = NPhiBins(eta+1); + int nDepth = NDepthBins(eta,-1); + int nPhi = NPhiBins(eta); //std::cout<<"problem 2"<cd(); - for (int eta=0; eta<29; ++eta) { - int nDepth = NDepthBins(eta+1,-1); - int nPhi = NPhiBins(eta+1); + for (int eta=etaMin_; eta<=etaMax_; ++eta) { + int nDepth = NDepthBins(eta,-1); + int nPhi = NPhiBins(eta); //sprintf(name, "Dir_muon_type_%s_ieta%d",type[i].c_str(), eta); //d_output_file[i][eta]= output_file->mkdir(name); //output_file->cd(name); @@ -1217,68 +1256,68 @@ void HBHEMuonOfflineAnalyzer::BookHistograms(const char* fname) { for (int depth=0; depthSumw2(); - sprintf (name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength", (eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (title, "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell) divided by Active Length", (eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength", eta, (depth+1), PHI0, type[i].c_str()); + sprintf (title, "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell) divided by Active Length", eta, (depth+1), PHI0, type[i].c_str()); h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih] = new TH1D(name, title,4000,0.0,10.0); h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Sumw2(); - sprintf (name, "h_active_length_Fill_%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (title, "active_length%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_active_length_Fill_%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); + sprintf (title, "active_length%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); h_active_length_Fill[i][ih] = new TH1D(name, title,20,0,20); h_active_length_Fill[i][ih]->Sumw2(); - sprintf (name, "h_p_muon_in_%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (title, "p_muon_in%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_p_muon_in_%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); + sprintf (title, "p_muon_in%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); h_p_muon_ineta[i][ih] = new TH1D(name, title, 500,0,500); h_p_muon_ineta[i][ih]->Sumw2(); - sprintf (name, "h_charge_signal_in_%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (name, "charge_signal_in_%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_charge_signal_in_%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); + sprintf (name, "charge_signal_in_%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); h_charge_signal[i][ih] = new TH1D(name, title, 500,0,500); h_charge_signal[i][ih]->Sumw2(); - sprintf (name, "h_charge_bg_in_%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (name, "charge_bg_in_%d_%d_%d_%s",(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_charge_bg_in_%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); + sprintf (name, "charge_bg_in_%d_%d_%d_%s", eta, (depth+1), PHI0, type[i].c_str()); h_charge_bg[i][ih] = new TH1D(name, title, 500,0,500); h_charge_bg[i][ih]->Sumw2(); ih++; - sprintf (name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell", -(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (title, "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell)", -(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell", -eta, (depth+1), PHI0, type[i].c_str()); + sprintf (title, "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell)", -eta, (depth+1), PHI0, type[i].c_str()); h_Hot_MuonEnergy_hcal_HotCell[i][ih] = new TH1D(name, title,4000,0.0,10.0); h_Hot_MuonEnergy_hcal_HotCell[i][ih] ->Sumw2(); - sprintf (name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength", -(eta+1), (depth+1),PHI0, type[i].c_str()); - sprintf (title, "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi=%d) for extrapolated %s muons (Hot Cell) divided by Active Length", -(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength", -eta, (depth+1),PHI0, type[i].c_str()); + sprintf (title, "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi=%d) for extrapolated %s muons (Hot Cell) divided by Active Length", -eta, (depth+1), PHI0, type[i].c_str()); h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih] = new TH1D(name, title,4000,0.0,10.0); h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Sumw2(); - sprintf (name, "h_active_length_Fill_%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (title, "active_length%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_active_length_Fill_%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); + sprintf (title, "active_length%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); h_active_length_Fill[i][ih] = new TH1D(name, title,20,0,20); h_active_length_Fill[i][ih]->Sumw2(); - sprintf (name, "h_p_muon_in_%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (title, "p_muon_in%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_p_muon_in_%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); + sprintf (title, "p_muon_in%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); h_p_muon_ineta[i][ih] = new TH1D(name, title, 500,0,500); h_p_muon_ineta[i][ih]->Sumw2(); - sprintf (name, "h_charge_signal_in_%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (name, "charge_signal_in_%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_charge_signal_in_%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); + sprintf (name, "charge_signal_in_%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); h_charge_signal[i][ih] = new TH1D(name, title, 500,0,500); h_charge_signal[i][ih]->Sumw2(); - sprintf (name, "h_charge_bg_in_%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); - sprintf (name, "charge_bg_in_%d_%d_%d_%s",-(eta+1), (depth+1), PHI0, type[i].c_str()); + sprintf (name, "h_charge_bg_in_%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); + sprintf (name, "charge_bg_in_%d_%d_%d_%s",-eta, (depth+1), PHI0, type[i].c_str()); h_charge_bg[i][ih] = new TH1D(name, title, 500,0,500); h_charge_bg[i][ih]->Sumw2(); @@ -1406,16 +1445,16 @@ void HBHEMuonOfflineAnalyzer::WriteHistograms() { std::cout<<"WriteHistograms"<Write();; output_file->cd(); - for (int eta=0; eta<29; ++eta) { - int nDepth = NDepthBins(eta+1,-1); - int nPhi = NPhiBins(eta+1); - sprintf(name, "Dir_muon_type_%s_ieta%d",type[i].c_str(), eta+1); - d_output_file[i][eta]= output_file->mkdir(name); + for (int eta=etaMin_; eta<=etaMax_; ++eta) { + int nDepth = NDepthBins(eta,-1); + int nPhi = NPhiBins(eta); + sprintf(name, "Dir_muon_type_%s_ieta%d",type[i].c_str(), eta); + d_output_file[i][eta-1]= output_file->mkdir(name); //output_file->cd(name); - d_output_file[i][eta]->cd(); + d_output_file[i][eta-1]->cd(); for (int depth=0; depthWrite(); @@ -1509,15 +1548,18 @@ int HBHEMuonOfflineAnalyzer::NDepthBins(int eta, int phi) { int nDepthR2[29]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,4,3,5,6,6,6,6,6,6,6,7,7,7,3}; // Run 3 scenario int nDepthR3[29]={4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,5,6,6,6,6,6,6,6,7,7,7,3}; + // Run 4 scenario + int nDepthR4[29]={4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,7,7,7,7,7,7,7,7,7,7,7,7,7}; // for 2019 scenario multi depth segmentation // int nDepth[29]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,5,5,5,5,5,5,5,5,5,5,5,5,5}; - // modeLHC = 0 --> nbin defined maxDepthHB/HE - // = 1 --> corresponds to Run 1 (valid till 2016) - // = 2 --> corresponds to Run 2 (2018 geometry) - // = 3 --> corresponds to Run 3 (post LS2) - // = 4 --> corresponds to 2017 (Plan 1) + // modeLHC_ = 0 --> nbin defined maxDepthHB/HE + // = 1 --> corresponds to Run 1 (valid till 2016) + // = 2 --> corresponds to Run 2 (2018 geometry) + // = 3 --> corresponds to Run 3 (post LS2) + // = 4 --> corresponds to 2017 (Plan 1) + // = 5 --> corresponds to Run 4 (post LS3) int nbin(0); - if (modeLHC == 0) { + if (modeLHC_ == 0) { if (eta<=15) { nbin = maxDepthHB_; } else if (eta == 16) { @@ -1525,13 +1567,13 @@ int HBHEMuonOfflineAnalyzer::NDepthBins(int eta, int phi) { } else { nbin = maxDepthHE_; } - } else if (modeLHC == 1) { + } else if (modeLHC_ == 1) { nbin = nDepthR1[eta-1]; - } else if (modeLHC == 2) { + } else if (modeLHC_ == 2) { nbin = nDepthR2[eta-1]; - } else if (modeLHC == 3) { + } else if (modeLHC_ == 3) { nbin = nDepthR3[eta-1]; - } else { + } else if (modeLHC_ == 4) { if (phi > 0) { if (eta >= 16 && phi >= 63 && phi <= 66) { nbin = nDepthR2[eta-1]; @@ -1545,11 +1587,18 @@ int HBHEMuonOfflineAnalyzer::NDepthBins(int eta, int phi) { nbin = nDepthR1[eta-1]; } } + } else { + if (eta > 0 && eta < 30) { + nbin = nDepthR4[eta-1]; + } else { + nbin = nDepthR4[28]; + } } return nbin; } int HBHEMuonOfflineAnalyzer::NPhiBins(int eta) { int nphi = (eta <= 20) ? 72 : 36; + if (modeLHC_ == 5 && eta > 16) nphi = 360; return nphi; } diff --git a/Calibration/HcalCalibAlgos/macros/OfflineMain.C b/Calibration/HcalCalibAlgos/macros/OfflineMain.C index 8c18f2057ccb7..b15ddf160912a 100644 --- a/Calibration/HcalCalibAlgos/macros/OfflineMain.C +++ b/Calibration/HcalCalibAlgos/macros/OfflineMain.C @@ -35,13 +35,18 @@ int main(Int_t argc, Char_t *argv[]) { << "Output FIle Name " << outFileName << std::endl << "Process Name " << processName << std::endl; - int flag(0), mode(4), maxDHB(5), maxDHE(7), runLo(0), runHi(9999999); - if (argc>4) flag = atoi(argv[4]); - if (argc>5) mode = atoi(argv[5]); - if (argc>6) maxDHB = atoi(argv[6]); - if (argc>7) maxDHE = atoi(argv[7]); - if (argc>8) runLo = atoi(argv[8]); - if (argc>9) runHi = atoi(argv[9]); + int flag(0), mode(4), maxDHB(5), maxDHE(7), runLo(1), int runHi(9999999); + int etaMin(1), etaMax(29); + if (argc>4) flag = atoi(argv[4]); + if (argc>5) mode = atoi(argv[5]); + if (argc>6) maxDHB = atoi(argv[6]); + if (argc>7) maxDHE = atoi(argv[7]); + if (argc>8) runLo = atoi(argv[8]); + if (argc>9) runHi = atoi(argv[9]); + if (argc>8) runLo = atoi(argv[8]); + if (argc>9) runHi = atoi(argv[9]); + if (argc>10) etaMin = atoi(argv[10]); + if (argc>11) etaMax = atoi(argv[11]); char treeName[400]; sprintf (treeName, "%s/TREE", processName); std::cout << "try to create a chain for " << treeName << std::endl; @@ -49,7 +54,7 @@ int main(Int_t argc, Char_t *argv[]) { if (FillChain(chain, inputFileList)) { - HBHEMuonOfflineAnalyzer* tree = new HBHEMuonOfflineAnalyzer(chain, outFileName, flag, mode, maxDHB, maxDHE, runLo, runHi); + HBHEMuonOfflineAnalyzer* tree = new HBHEMuonOfflineAnalyzer(chain, outFileName, flag, mode, maxDHB, maxDHE, runLo, runHi, etaMin, etaMax); tree->Loop(); } @@ -59,9 +64,6 @@ int main(Int_t argc, Char_t *argv[]) { Bool_t FillChain(TChain *chain, const char* inputFileList) { ifstream infile(inputFileList); -#if 0 - std::string buffer; -#endif if (!infile.is_open()) { std::cerr << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl; return kFALSE; @@ -75,15 +77,7 @@ Bool_t FillChain(TChain *chain, const char* inputFileList) { std::cout << "Adding " << buffer << " to chain" << std::endl; chain->Add(buffer); } -#if 0 - while (1) { - infile >> buffer; - if (!infile.good()) break; - chain->Add(buffer.c_str()); - } -#endif std::cout << "No. of Entries in this tree : " << chain->GetEntries() << std::endl; infile.close(); return kTRUE; } - diff --git a/Calibration/HcalCalibAlgos/test/HcalIsoTrkAnalyzer.cc b/Calibration/HcalCalibAlgos/test/HcalIsoTrkAnalyzer.cc index 173baf480316e..28d7f450725e3 100644 --- a/Calibration/HcalCalibAlgos/test/HcalIsoTrkAnalyzer.cc +++ b/Calibration/HcalCalibAlgos/test/HcalIsoTrkAnalyzer.cc @@ -171,7 +171,7 @@ class HcalIsoTrkAnalyzer : public edm::one::EDAnalyzer *t_ietaAll, *t_ietaGood; + std::vector *t_ietaAll, *t_ietaGood, *t_trackType; }; static const bool useL1EventSetup(false); @@ -449,7 +449,7 @@ void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const std::vector vecL1, vecL3; t_Tracks = trkCollection->size(); t_TracksProp = trkCaloDirections.size(); - t_ietaAll->clear(); t_ietaGood->clear(); + t_ietaAll->clear(); t_ietaGood->clear(); t_trackType->clear(); #ifdef EDM_ML_DEBUG edm::LogVerbatim("HcalIsoTrack") << "# of propagated tracks " << t_TracksProp << " out of " << t_Tracks << " with Trigger " @@ -705,8 +705,10 @@ void HcalIsoTrkAnalyzer::beginJob() { t_hltbits = new std::vector(); t_ietaAll = new std::vector(); t_ietaGood = new std::vector(); + t_trackType = new std::vector(); tree2->Branch("t_ietaAll", "std::vector", &t_ietaAll); tree2->Branch("t_ietaGood", "std::vector", &t_ietaGood); + tree2->Branch("t_trackType", "std::vector", &t_trackType); tree2->Branch("t_hltbits", "std::vector", &t_hltbits); } @@ -967,12 +969,15 @@ std::array HcalIsoTrkAnalyzer::fillTree(std::vector< math::XYZTLorentzVec if (t_p>pTrackMin_) { tree->Fill(); nSave++; + int type(0); if (t_eMipDR < 1.0) { - if (t_hmaxNearP < 2.0) ++nTight; - if (t_hmaxNearP < 10.0) ++nLoose; + if (t_hmaxNearP < 10.0) { ++nLoose; type = 1;} + if (t_hmaxNearP < 2.0) { ++nTight; type = 2;} } - if (t_p > 40.0 && t_p <= 60.0 && t_selectTk) + if (t_p > 40.0 && t_p <= 60.0 && t_selectTk) { t_ietaGood->emplace_back(t_ieta); + t_trackType->emplace_back(type); + } #ifdef EDM_ML_DEBUG for (unsigned int k=0; ksize(); k++) { edm::LogVerbatim("HcalIsoTrack") << "trigger bit is = " @@ -983,7 +988,7 @@ std::array HcalIsoTrkAnalyzer::fillTree(std::vector< math::XYZTLorentzVec } } } - std::array i3{ {nSave,nTight,nLoose} }; + std::array i3{ {nSave,nLoose,nTight} }; return i3; } From b7fda3063675cf5b822bfb4a7aa2a2c5ff2fe4e6 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sun, 19 Nov 2017 05:17:54 +0100 Subject: [PATCH 3/3] Update the macros --- .../HcalCalibAlgos/macros/AnalyzeLepTree.C | 78 +++++++++++++------ .../macros/HBHEMuonOfflineAnalyzer.C | 56 ++++++++++--- 2 files changed, 100 insertions(+), 34 deletions(-) diff --git a/Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C b/Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C index e932a45d816de..8f9a7da1cbf24 100644 --- a/Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C +++ b/Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C @@ -114,6 +114,7 @@ private: Int_t t_nvtx; Double_t t_p; std::vector *t_ene; + std::vector *t_enec; std::vector *t_charge; std::vector *t_actln; std::vector *t_depth; @@ -124,6 +125,7 @@ private: TBranch *b_t_nvtx; //! TBranch *b_t_p; //! TBranch *b_t_ene; //! + TBranch *b_t_enec; //! TBranch *b_t_charge; //! TBranch *b_t_actln; //! TBranch *b_t_depth; //! @@ -135,6 +137,7 @@ private: std::map h_pnv_; std::map h_p2_, h_nv2_; std::map h_Energy_, h_Ecorr_, h_Charge_, h_Chcorr_; + std::map h_EnergyC_, h_EcorrC_; }; AnalyzeLepTree::AnalyzeLepTree(TTree *tree, int mode1, @@ -185,6 +188,7 @@ void AnalyzeLepTree::Init(TTree *tree) { // Set object pointer t_ene = 0; + t_enec = 0; t_charge = 0; t_actln = 0; t_depth = 0; @@ -199,6 +203,7 @@ void AnalyzeLepTree::Init(TTree *tree) { fChain->SetBranchAddress("t_nvtx", &t_nvtx, &b_t_nvtx); fChain->SetBranchAddress("t_p", &t_p, &b_t_p); fChain->SetBranchAddress("t_ene", &t_ene, &b_t_ene); + fChain->SetBranchAddress("t_enec", &t_enec, &b_t_enec); fChain->SetBranchAddress("t_charge", &t_charge, &b_t_charge); fChain->SetBranchAddress("t_actln", &t_actln, &b_t_actln); fChain->SetBranchAddress("t_depth", &t_depth, &b_t_depth); @@ -289,41 +294,51 @@ void AnalyzeLepTree::Loop() { int depth = (*t_depth)[k]; unsigned int id = packID(zside,eta,phi,depth+1,vbin,pbin); double ene = (*t_ene)[k]; + double enec = (*t_enec)[k]; double charge = (*t_charge)[k]; double actl = (*t_actln)[k]; if (ene > 0 && actl > 0 && charge > 0) { std::map::iterator it1 = h_Energy_.find(id); - if (it1 != h_Energy_.end()) (it1->second)->Fill(ene); + if (it1 != h_Energy_.end()) (it1->second)->Fill(ene); std::map::iterator it2 = h_Ecorr_.find(id); - if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl); - std::map::iterator it3 = h_Charge_.find(id); - if (it3 != h_Charge_.end()) (it3->second)->Fill(charge); - std::map::iterator it4 = h_Chcorr_.find(id); - if (it4 != h_Chcorr_.end()) (it4->second)->Fill(charge/actl); + if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl); + std::map::iterator it3 = h_EnergyC_.find(id); + if (it3 != h_EnergyC_.end()) (it3->second)->Fill(enec); + std::map::iterator it4 = h_EcorrC_.find(id); + if (it4 != h_EcorrC_.end()) (it4->second)->Fill(enec/actl); + std::map::iterator it5 = h_Charge_.find(id); + if (it5 != h_Charge_.end()) (it5->second)->Fill(charge); + std::map::iterator it6 = h_Chcorr_.find(id); + if (it6 != h_Chcorr_.end()) (it6->second)->Fill(charge/actl); /* - if ((eta>20 && (t_iphi > 35)) || (t_iphi > 71)) std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth+1 << ":" << vbin << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_Charge_.end()) << ":" << (it4 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl; - if ((it1 == h_Energy_.end()) || (it2 == h_Ecorr_.end()) || (it3 == h_Charge_.end()) || (it4 == h_Chcorr_.end())) std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth+1 << ":" << vbin << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_Charge_.end()) << ":" << (it4 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl; + if ((eta>20 && (t_iphi > 35)) || (t_iphi > 71)) std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth+1 << ":" << vbin << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_EnergyC_.end()) << ":" << (it4 != h_EcorrC_.end()) << ":" << (it5 != h_Charge_.end()) << ":" << (it6 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl; + if ((it1 == h_Energy_.end()) || (it2 == h_Ecorr_.end()) || (it3 == h_EnergyC_.end()) || (it4 == h_EcorrC_.end()) || (it5 == h_Charge_.end()) || (it6 == h_Chcorr_.end())) std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth+1 << ":" << vbin << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_Charge_.end()) << ":" << (it4 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl; */ } } } else { - double ene(0), actl(0), charge(0); + double ene(0), enec(0), actl(0), charge(0); unsigned int id = packID(zside,eta,phi,1,vbin,pbin); for (unsigned int k=0; ksize(); ++k) { - if ((*t_ene)[k] > 0 && (*t_actln)[k] > 0 && (*t_charge)[k] > 0) { + if ((*t_ene)[k] > 0 && (*t_actln)[k] > 0) { ene += (*t_ene)[k]; + enec += (*t_enec)[k]; charge += (*t_charge)[k]; actl += (*t_actln)[k]; } } std::map::iterator it1 = h_Energy_.find(id); - if (it1 != h_Energy_.end()) (it1->second)->Fill(ene); + if (it1 != h_Energy_.end()) (it1->second)->Fill(ene); std::map::iterator it2 = h_Ecorr_.find(id); - if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl); - std::map::iterator it3 = h_Charge_.find(id); - if (it3 != h_Charge_.end()) (it3->second)->Fill(charge); - std::map::iterator it4 = h_Chcorr_.find(id); - if (it4 != h_Chcorr_.end()) (it4->second)->Fill(charge/actl); + if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl); + std::map::iterator it3 = h_EnergyC_.find(id); + if (it3 != h_EnergyC_.end()) (it3->second)->Fill(ene); + std::map::iterator it4 = h_Ecorr_.find(id); + if (it4 != h_EcorrC_.end()) (it4->second)->Fill(ene/actl); + std::map::iterator it5 = h_Charge_.find(id); + if (it5 != h_Charge_.end()) (it5->second)->Fill(charge); + std::map::iterator it6 = h_Chcorr_.find(id); + if (it6 != h_Chcorr_.end()) (it6->second)->Fill(charge/actl); } } } @@ -400,6 +415,7 @@ void AnalyzeLepTree::bookHisto() { } } else { h_Energy_.clear(); h_Ecorr_.clear(); h_Charge_.clear(); h_Chcorr_.clear(); + h_EnergyC_.clear(); h_EcorrC_.clear(); for (int ieta=-26; ieta<=26; ++ieta) { if (ieta != 0) { int zside = (ieta>0) ? 1 : -1; @@ -449,6 +465,14 @@ void AnalyzeLepTree::bookHisto() { sprintf (title,"Active length corrected energy for %s %s %s %s %s (GeV/cm)", etas, phis, deps, ps, vtx); h_Ecorr_[id] = new TH1D(name,title,4000,0.0,10.0); h_Ecorr_[id]->Sumw2(); + sprintf (name,"EdepCE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin); + sprintf (title,"Response Corrected deposited energy for %s %s %s %s %s (GeV)", etas, phis, deps, ps, vtx); + h_EnergyC_[id] = new TH1D(name,title,4000,0.0,10.0); + h_EnergyC_[id]->Sumw2(); + sprintf (name,"EcorCE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin); + sprintf (title,"Response and active length corrected energy for %s %s %s %s %s (GeV/cm)", etas, phis, deps, ps, vtx); + h_EcorrC_[id] = new TH1D(name,title,4000,0.0,10.0); + h_EcorrC_[id]->Sumw2(); sprintf (name,"ChrgE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin); sprintf (title,"Measured charge for %s %s %s %s %s (cm)", etas, phis, deps, ps, vtx); h_Charge_[id] = new TH1D(name,title,20,0.0,20.0); @@ -505,13 +529,17 @@ void AnalyzeLepTree::writeHisto(std::string fname) { unsigned int id = packID(zside,eta,phi,depth+1,vbin,pbin); std::map::const_iterator itr; itr = h_Energy_.find(id); - if (itr != h_Energy_.end()) (itr->second)->Write(); + if (itr != h_Energy_.end()) (itr->second)->Write(); itr = h_Ecorr_.find(id); - if (itr != h_Ecorr_.end()) (itr->second)->Write(); + if (itr != h_Ecorr_.end()) (itr->second)->Write(); + itr = h_EnergyC_.find(id); + if (itr != h_EnergyC_.end()) (itr->second)->Write(); + itr = h_EcorrC_.find(id); + if (itr != h_EcorrC_.end()) (itr->second)->Write(); itr = h_Charge_.find(id); - if (itr != h_Charge_.end()) (itr->second)->Write(); + if (itr != h_Charge_.end()) (itr->second)->Write(); itr = h_Chcorr_.find(id); - if (itr != h_Chcorr_.end()) (itr->second)->Write(); + if (itr != h_Chcorr_.end()) (itr->second)->Write(); } } } @@ -548,10 +576,12 @@ std::vector AnalyzeLepTree::plotHisto(bool drawStatBox, int type, bool doEnL = ((type/2)%2 > 0); bool doChg = ((type/4)%2 > 0); bool doChL = ((type/8)%2 > 0); - if (doEn) plotHisto(h_Energy_, phi0, depth0, cvs, save); - if (doEnL) plotHisto(h_Ecorr_, phi0, depth0, cvs, save); - if (doChg) plotHisto(h_Charge_, phi0, depth0, cvs, save); - if (doChL) plotHisto(h_Chcorr_, phi0, depth0, cvs, save); + if (doEn) plotHisto(h_Energy_, phi0, depth0, cvs, save); + if (doEn) plotHisto(h_EnergyC_, phi0, depth0, cvs, save); + if (doEnL) plotHisto(h_Ecorr_, phi0, depth0, cvs, save); + if (doEnL) plotHisto(h_EcorrC_, phi0, depth0, cvs, save); + if (doChg) plotHisto(h_Charge_, phi0, depth0, cvs, save); + if (doChL) plotHisto(h_Chcorr_, phi0, depth0, cvs, save); } return cvs; } diff --git a/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C b/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C index 82c1d33dbd5b7..2a4bd5905abdc 100644 --- a/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C +++ b/Calibration/HcalCalibAlgos/macros/HBHEMuonOfflineAnalyzer.C @@ -362,7 +362,7 @@ private: TTree *outtree_; int t_ieta, t_iphi, t_nvtx; double t_p; - std::vector t_ene, t_actln, t_charge; + std::vector t_ene, t_enec, t_actln, t_charge; std::vector t_depth; TH1D *h_Pt_Muon[3], *h_Eta_Muon[3], *h_Phi_Muon[3], *h_P_Muon[3]; @@ -760,7 +760,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (debug_) std::cout << "Run " << Run_No << " Event " << Event_No << " Muons " << pt_of_muon->size() << std::endl; for (unsigned int ml = 0; ml< pt_of_muon->size(); ml++) { - t_ene.clear(); t_charge.clear(); t_actln.clear(); t_depth.clear(); + t_ene.clear(); t_enec.clear(); t_charge.clear(); t_actln.clear(); + t_depth.clear(); if(debug_) std::cout << "ecal_det_id " << ecal_detID->at(ml) << std::endl; @@ -848,10 +849,13 @@ void HBHEMuonOfflineAnalyzer::Loop() { h_HotCell[cut]->Fill(hcal_cellHot->at(ml)); if (mergeDepth_) { double en1(0), en2(0), energyFill(0), chargeS(0), chargeBG(0); + double enh(0), enc(0); for (int dep=0; depat(ml) : hcal_edepth1->at(ml)); en2 += ((useCorrect_) ? hcal_edepthHotCorrect1->at(ml) : hcal_edepthHot1->at(ml)); + enh += (hcal_edepthHot1->at(ml)); + enc += (hcal_edepthHotCorrect1->at(ml)); energyFill += (hcal_activeHotL1->at(ml)); chargeS += (hcal_cdepthHot1->at(ml)); chargeBG += (hcal_cdepthHotBG1->at(ml)); @@ -859,6 +863,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { } else if (dep == 1) { en1 += ((useCorrect_) ? hcal_edepthCorrect2->at(ml) : hcal_edepth2->at(ml)); en2 += ((useCorrect_) ? hcal_edepthHotCorrect2->at(ml) : hcal_edepthHot2->at(ml)); + enh += (hcal_edepthHot2->at(ml)); + enc += (hcal_edepthHotCorrect2->at(ml)); energyFill += (hcal_activeHotL2->at(ml)); chargeS += (hcal_cdepthHot2->at(ml)); chargeBG += (hcal_cdepthHotBG2->at(ml)); @@ -866,6 +872,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { } else if (dep == 2) { en1 += ((useCorrect_) ? hcal_edepthCorrect3->at(ml) : hcal_edepth3->at(ml)); en2 += ((useCorrect_) ? hcal_edepthHotCorrect3->at(ml) : hcal_edepthHot3->at(ml)); + enh += (hcal_edepthHot3->at(ml)); + enc += (hcal_edepthHotCorrect3->at(ml)); energyFill += (hcal_activeHotL3->at(ml)); chargeS += (hcal_cdepthHot3->at(ml)); chargeBG += (hcal_cdepthHotBG3->at(ml)); @@ -873,6 +881,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { } else if (dep == 3) { en1 += ((useCorrect_) ? hcal_edepthCorrect4->at(ml) : hcal_edepth4->at(ml)); en2 += ((useCorrect_) ? hcal_edepthHotCorrect4->at(ml) : hcal_edepthHot4->at(ml)); + enh += (hcal_edepthHot4->at(ml)); + enc += (hcal_edepthHotCorrect4->at(ml)); energyFill = hcal_activeHotL4->at(ml); chargeS += (hcal_cdepthHot4->at(ml)); chargeBG += (hcal_cdepthHotBG4->at(ml)); @@ -881,6 +891,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (hcal_edepthCorrect5->size() > ml) { en1 += ((useCorrect_) ? hcal_edepthCorrect5->at(ml) : hcal_edepth5->at(ml)); en2 += ((useCorrect_) ? hcal_edepthHotCorrect5->at(ml) : hcal_edepthHot5->at(ml)); + enh += (hcal_edepthHot5->at(ml)); + enc += (hcal_edepthHotCorrect5->at(ml)); energyFill += (hcal_activeHotL5->at(ml)); chargeS += (hcal_cdepthHot5->at(ml)); chargeBG += (hcal_cdepthHotBG5->at(ml)); @@ -890,6 +902,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (hcal_edepthCorrect6->size() > ml) { en1 += ((useCorrect_) ? hcal_edepthCorrect6->at(ml) : hcal_edepth6->at(ml)); en2 += ((useCorrect_) ? hcal_edepthHotCorrect6->at(ml) : hcal_edepthHot6->at(ml)); + enh += (hcal_edepthHot6->at(ml)); + enc += (hcal_edepthHotCorrect6->at(ml)); energyFill += (hcal_activeHotL6->at(ml)); chargeS += (hcal_cdepthHot6->at(ml)); chargeBG += (hcal_cdepthHotBG6->at(ml)); @@ -899,6 +913,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (hcal_edepthCorrect7->size() > ml) { en1 += ((useCorrect_) ? hcal_edepthCorrect7->at(ml) : hcal_edepth7->at(ml)); en2 += ((useCorrect_) ? hcal_edepthHotCorrect7->at(ml) : hcal_edepthHot7->at(ml)); + enh += (hcal_edepthHot7->at(ml)); + enc += (hcal_edepthHotCorrect7->at(ml)); energyFill += (hcal_activeHotL7->at(ml)); chargeS += (hcal_cdepthHot7->at(ml)); chargeBG += (hcal_cdepthHotBG7->at(ml)); @@ -911,9 +927,9 @@ void HBHEMuonOfflineAnalyzer::Loop() { std::cout << "Matched Id " << matchedId->at(ml) << " Hot " << hcal_cellHot->at(ml) << " eta " << etaHcal << ":" << eta << " phi " << phiHcal << ":" << PHI - << " Index " << ind << " E " << en2 << ":" - << energyFill << " Charge " << chargeS << ":" - << chargeBG << std::endl; + << " Index " << ind << " E " << en1 << ":" << en2 + << ":" << enh << ":" << enc << " L " << energyFill + << " Charge " << chargeS << ":" << chargeBG <at(ml))) continue; if (hcal_cellHot->at(ml)==1) { if (energyFill > 0) { @@ -924,7 +940,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { h_charge_signal[cut][ind]->Fill(chargeS); h_charge_bg[cut][ind]->Fill(chargeBG); - t_ene.push_back(en2); + t_ene.push_back(enh); + t_enec.push_back(enc); t_charge.push_back(chargeS); t_actln.push_back(energyFill); t_depth.push_back(0); @@ -938,28 +955,37 @@ void HBHEMuonOfflineAnalyzer::Loop() { if(debug_) std::cout<<"dep:"<at(ml) : hcal_edepth1->at(ml); en2 = (useCorrect_) ? hcal_edepthHotCorrect1->at(ml) : hcal_edepthHot1->at(ml); + enh = (hcal_edepthHot1->at(ml)); + enc = (hcal_edepthHotCorrect1->at(ml)); energyFill = hcal_activeHotL1->at(ml); chargeS = (hcal_cdepthHot1->at(ml)); chargeBG = (hcal_cdepthHotBG1->at(ml)); } else if (dep == 1) { en1 = (useCorrect_) ? hcal_edepthCorrect2->at(ml) : hcal_edepth2->at(ml); en2 = (useCorrect_) ? hcal_edepthHotCorrect2->at(ml) : hcal_edepthHot2->at(ml); + enh = (hcal_edepthHot2->at(ml)); + enc = (hcal_edepthHotCorrect2->at(ml)); energyFill = hcal_activeHotL2->at(ml); chargeS = (hcal_cdepthHot2->at(ml)); chargeBG = (hcal_cdepthHotBG2->at(ml)); } else if (dep == 2) { en1 = (useCorrect_) ? hcal_edepthCorrect3->at(ml) : hcal_edepth3->at(ml); en2 = (useCorrect_) ? hcal_edepthHotCorrect3->at(ml) : hcal_edepthHot3->at(ml); + enh = (hcal_edepthHot3->at(ml)); + enc = (hcal_edepthHotCorrect3->at(ml)); energyFill = hcal_activeHotL3->at(ml); chargeS = (hcal_cdepthHot3->at(ml)); chargeBG = (hcal_cdepthHotBG3->at(ml)); } else if (dep == 3) { en1 = (useCorrect_) ? hcal_edepthCorrect4->at(ml) : hcal_edepth4->at(ml); en2 = (useCorrect_) ? hcal_edepthHotCorrect4->at(ml) : hcal_edepthHot4->at(ml); + enh = (hcal_edepthHot4->at(ml)); + enc = (hcal_edepthHotCorrect4->at(ml)); energyFill = hcal_activeHotL4->at(ml); chargeS = (hcal_cdepthHot4->at(ml)); chargeBG = (hcal_cdepthHotBG4->at(ml)); @@ -967,6 +993,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (hcal_edepthCorrect5->size() > ml) { en1 = (useCorrect_) ? hcal_edepthCorrect5->at(ml) : hcal_edepth5->at(ml); en2 = (useCorrect_) ? hcal_edepthHotCorrect5->at(ml) : hcal_edepthHot5->at(ml); + enh = (hcal_edepthHot5->at(ml)); + enc = (hcal_edepthHotCorrect5->at(ml)); energyFill = hcal_activeHotL5->at(ml); chargeS = (hcal_cdepthHot5->at(ml)); chargeBG = (hcal_cdepthHotBG5->at(ml)); @@ -975,6 +1003,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (hcal_edepthCorrect6->size() > ml) { en1 = (useCorrect_) ? hcal_edepthCorrect6->at(ml) : hcal_edepth6->at(ml); en2 = (useCorrect_) ? hcal_edepthHotCorrect6->at(ml) : hcal_edepthHot6->at(ml); + enh = (hcal_edepthHot6->at(ml)); + enc = (hcal_edepthHotCorrect6->at(ml)); energyFill = hcal_activeHotL6->at(ml); chargeS = (hcal_cdepthHot6->at(ml)); chargeBG = (hcal_cdepthHotBG6->at(ml)); @@ -983,6 +1013,8 @@ void HBHEMuonOfflineAnalyzer::Loop() { if (hcal_edepthCorrect7->size() > ml) { en1 = (useCorrect_) ? hcal_edepthCorrect7->at(ml) : hcal_edepth7->at(ml); en2 = (useCorrect_) ? hcal_edepthHotCorrect7->at(ml) : hcal_edepthHot7->at(ml); + enh = (hcal_edepthHot7->at(ml)); + enc = (hcal_edepthHotCorrect7->at(ml)); energyFill = hcal_activeHotL7->at(ml); chargeS = (hcal_cdepthHot7->at(ml)); chargeBG = (hcal_cdepthHotBG7->at(ml)); @@ -1004,8 +1036,9 @@ void HBHEMuonOfflineAnalyzer::Loop() { << hcal_cellHot->at(ml) << " eta " << etaHcal << ":" << eta << " phi " << phiHcal << ":" << PHI << " depth " << dep << " Index " << ind << " E " - << en2 << ":" << energyFill << " Charge " - << chargeS << ":" << chargeBG << std::endl; + << en1 << ":" << en2 << ":" << enh << ":" << enc + << " L " << energyFill << " Charge " << chargeS + << ":" << chargeBG << std::endl; if (!(matchedId->at(ml))) continue; if (ok1) { if (debug_) std::cout<<"enter ok1"<Fill(p_of_muon->at(ml)); h_charge_signal[cut][ind]->Fill(chargeS); h_charge_bg[cut][ind]->Fill(chargeBG); - t_ene.push_back(en2); + t_ene.push_back(enh); + t_enec.push_back(enc); t_charge.push_back(chargeS); t_actln.push_back(energyFill); // added depth vector AmanKalsi @@ -1026,6 +1060,7 @@ void HBHEMuonOfflineAnalyzer::Loop() { fillTree = true; } else { t_ene.push_back(-999.0); + t_enec.push_back(-999.0); t_charge.push_back(-999.0); t_actln.push_back(-999.0); t_depth.push_back(-999.0); @@ -1079,6 +1114,7 @@ void HBHEMuonOfflineAnalyzer::BookHistograms(const char* fname) { outtree_->Branch("t_nvtx", &t_nvtx); outtree_->Branch("t_p", &t_p); outtree_->Branch("t_ene", &t_ene); + outtree_->Branch("t_enec", &t_enec); outtree_->Branch("t_charge", &t_charge); outtree_->Branch("t_actln", &t_actln); outtree_->Branch("t_depth", &t_depth);