From fe1c621eb674d30e704f039107d5cc093ff57d26 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Tue, 19 Jul 2022 08:06:07 +0200 Subject: [PATCH] Make use of MessageLogger wherever possible in Calibration/IsolatedParticles and HcalIsolatedTrackReco packages Code check --- .../plugins/ECALRegFEDSelector.cc | 5 +- .../plugins/SubdetFEDSelector.cc | 5 +- .../interface/eHCALMatrixExtra.h | 2 + .../interface/eHCALMatrixExtra.icc | 130 +- .../macros/HBHEMuonOfflineSimAnalyzer.C | 1222 +++++----- .../IsolatedParticles/macros/PlotTracks.C | 2054 ++++++++++++----- .../IsolatedParticles/macros/PlotVecGeom.C | 627 +++-- .../IsolatedParticles/macros/StudyHLT.C | 541 +++-- .../plugins/IsolatedTracksNxN.cc | 2 +- 9 files changed, 2922 insertions(+), 1666 deletions(-) diff --git a/Calibration/HcalIsolatedTrackReco/plugins/ECALRegFEDSelector.cc b/Calibration/HcalIsolatedTrackReco/plugins/ECALRegFEDSelector.cc index 732115f407d84..27991d927c317 100644 --- a/Calibration/HcalIsolatedTrackReco/plugins/ECALRegFEDSelector.cc +++ b/Calibration/HcalIsolatedTrackReco/plugins/ECALRegFEDSelector.cc @@ -5,6 +5,7 @@ #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/one/EDProducer.h" #include "FWCore/Framework/interface/Event.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" @@ -109,8 +110,8 @@ void ECALRegFEDSelector::produce(edm::Event& iEvent, const edm::EventSetup& iSet // this fed has data -- lets copy it FEDRawData& fedDataProd = producedData->FEDData(j); if (fedDataProd.size() != 0) { - // std::cout << " More than one FEDRawDataCollection with data in FED "; - // std::cout << j << " Skipping the 2nd\n"; + edm::LogVerbatim("HcalCalib") << " More than one FEDRawDataCollection with data in FED " << j + << " Skipping the 2nd"; continue; } fedDataProd.resize(size); diff --git a/Calibration/HcalIsolatedTrackReco/plugins/SubdetFEDSelector.cc b/Calibration/HcalIsolatedTrackReco/plugins/SubdetFEDSelector.cc index 907bc0db21e17..2569467652b13 100644 --- a/Calibration/HcalIsolatedTrackReco/plugins/SubdetFEDSelector.cc +++ b/Calibration/HcalIsolatedTrackReco/plugins/SubdetFEDSelector.cc @@ -2,6 +2,7 @@ #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/global/EDProducer.h" #include "FWCore/Framework/interface/Event.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" @@ -184,8 +185,8 @@ void SubdetFEDSelector::produce(edm::StreamID, edm::Event& iEvent, const edm::Ev // this fed has data -- lets copy it FEDRawData& fedDataProd = producedData->FEDData(j); if (fedDataProd.size() != 0) { - // std::cout << " More than one FEDRawDataCollection with data in FED "; - // std::cout << j << " Skipping the 2nd\n"; + edm::LogVerbatim("HcalCalib") << " More than one FEDRawDataCollection with data in FED " << j + << " Skipping the 2nd"; continue; } fedDataProd.resize(size); diff --git a/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h b/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h index d0e5f6e3a2285..5bdfec36cce29 100644 --- a/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h +++ b/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h @@ -20,10 +20,12 @@ Created: August 2009 // system include files #include #include +#include #include // user include files #include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" #include "DataFormats/Common/interface/Handle.h" #include "DataFormats/Candidate/interface/Candidate.h" #include "DataFormats/DetId/interface/DetId.h" diff --git a/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.icc b/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.icc index 7c2444d694789..c5c4c2bea2b99 100644 --- a/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.icc +++ b/Calibration/IsolatedParticles/interface/eHCALMatrixExtra.icc @@ -42,8 +42,8 @@ namespace spr { HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1); DetId det(hcid.rawId()); if (debug) - std::cout << "Inside eHCALmatrixTotal " << ieta << "X" << iphi << " Inclusion of HO Flag " << includeHO - << std::endl; + edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrixTotal " << ieta << "X" << iphi << " Inclusion of HO Flag " + << includeHO; spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug); @@ -69,8 +69,8 @@ namespace spr { itrym = itry; } if (debug) - std::cout << "eHCALMatrix::Total energy " << energy << " Max " << energySum << "in trial # " << itrym - << std::endl; + edm::LogVerbatim("IsoTrack") << "eHCALMatrix::Total energy " << energy << " Max " << energySum << "in trial # " + << itrym; } return std::pair(energySum, itrym); @@ -93,12 +93,12 @@ namespace spr { HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1); DetId det(hcid.rawId()); if (debug) - std::cout << "Inside energyHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 << " Inclusion of HO Flag " - << includeHO << std::endl; + edm::LogVerbatim("IsoTrack") << "Inside energyHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1 + << " Inclusion of HO Flag " << includeHO; double energy = 0; if (ieta > 2) { - std::cout << " no matrix > 5x5 is coded ! " << std::endl; + edm::LogVerbatim("IsoTrack") << " no matrix > 5x5 is coded ! "; return energy; } // used dets during building the matrix @@ -113,23 +113,25 @@ namespace spr { // go to east from central tower vNeighboursDetId.clear(); if (debug) - std::cout << "hitHCALMatrix:: det " << (HcalDetId)det << std::endl; + edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: det " << (HcalDetId)det; vNeighboursDetId = topology->east(det); if (debug) { - std::cout << "Neighbour:: East"; + std::ostringstream st1; + st1 << "Neighbour:: East"; for (unsigned int i = 0; i < vNeighboursDetId.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); if (ieta == 2) { for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) { std::vector vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]); if (debug) { - std::cout << "Neighbour:: East"; + std::ostringstream st1; + st1 << "Neighbour:: East"; for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetIdc[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetIdc[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); } @@ -139,20 +141,22 @@ namespace spr { vNeighboursDetId.clear(); vNeighboursDetId = topology->west(det); if (debug) { - std::cout << "Neighbour:: West"; + std::ostringstream st1; + st1 << "Neighbour:: West"; for (unsigned int i = 0; i < vNeighboursDetId.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); if (ieta == 2) { for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) { std::vector vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]); if (debug) { - std::cout << "Neighbour:: West"; + std::ostringstream st1; + st1 << "Neighbour:: West"; for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetIdc[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetIdc[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); } @@ -163,13 +167,14 @@ namespace spr { DetId detnorth = det; for (int inorth = 0; inorth < iphi; inorth++) { if (debug) - std::cout << "hitHCALMatrix:: detNorth " << (HcalDetId)detnorth << std::endl; + edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: detNorth " << (HcalDetId)detnorth; std::vector NorthDetId = topology->north(detnorth); if (debug) { - std::cout << "Neighbour:: North"; + std::ostringstream st1; + st1 << "Neighbour:: North"; for (unsigned int i = 0; i < NorthDetId.size(); i++) - std::cout << " " << (HcalDetId)NorthDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)NorthDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } if (!(NorthDetId.empty())) { energy += energyHCAL(NorthDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); @@ -178,20 +183,22 @@ namespace spr { vNeighboursDetId.clear(); vNeighboursDetId = topology->east(NorthDetId[0]); if (debug) { - std::cout << "Neighbour:: East"; + std::ostringstream st1; + st1 << "Neighbour:: East"; for (unsigned int i = 0; i < vNeighboursDetId.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); if (ieta == 2) { for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) { std::vector vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]); if (debug) { - std::cout << "Neighbour:: East"; + std::ostringstream st1; + st1 << "Neighbour:: East"; for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetIdc[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetIdc[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); @@ -202,20 +209,22 @@ namespace spr { vNeighboursDetId.clear(); vNeighboursDetId = topology->west(NorthDetId[0]); if (debug) { - std::cout << "Neighbour:: West"; + std::ostringstream st1; + st1 << "Neighbour:: West"; for (unsigned int i = 0; i < vNeighboursDetId.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); if (ieta == 2) { for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) { std::vector vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]); if (debug) { - std::cout << "Neighbour:: West"; + std::ostringstream st1; + st1 << "Neighbour:: West"; for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetIdc[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetIdc[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); @@ -232,16 +241,17 @@ namespace spr { DetId detsouth = det; for (int isouth = 0; isouth < iphi; isouth++) { if (debug) - std::cout << "hitHCALMatrix:: detSouth " << (HcalDetId)detsouth << std::endl; + edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: detSouth " << (HcalDetId)detsouth; vNeighboursDetId.clear(); std::vector SouthDetId = topology->south(detsouth); if (!(SouthDetId.empty())) { energy += energyHCAL(SouthDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); if (debug) { - std::cout << "Neighbour:: South"; + std::ostringstream st1; + st1 << "Neighbour:: South"; for (unsigned int i = 0; i < SouthDetId.size(); i++) - std::cout << " " << (HcalDetId)SouthDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)SouthDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } // go to east @@ -250,20 +260,22 @@ namespace spr { vNeighboursDetId.clear(); vNeighboursDetId = topology->east(SouthDetId[0]); if (debug) { - std::cout << "Neighbour:: East"; + std::ostringstream st1; + st1 << "Neighbour:: East"; for (unsigned int i = 0; i < vNeighboursDetId.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); if (ieta == 2) { for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) { std::vector vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]); if (debug) { - std::cout << "Neighbour:: East"; + std::ostringstream st1; + st1 << "Neighbour:: East"; for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetIdc[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetIdc[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); @@ -274,20 +286,22 @@ namespace spr { vNeighboursDetId.clear(); vNeighboursDetId = topology->west(SouthDetId[0]); if (debug) { - std::cout << "Neighbour:: West"; + std::ostringstream st1; + st1 << "Neighbour:: West"; for (unsigned int i = 0; i < vNeighboursDetId.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetId[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetId[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); if (ieta == 2) { for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) { std::vector vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]); if (debug) { - std::cout << "Neighbour:: West"; + std::ostringstream st1; + st1 << "Neighbour:: West"; for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++) - std::cout << " " << (HcalDetId)vNeighboursDetIdc[i]; - std::cout << std::endl; + st1 << " " << (HcalDetId)vNeighboursDetIdc[i]; + edm::LogVerbatim("IsoTrack") << st1.str(); } energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug); @@ -375,7 +389,7 @@ namespace spr { } } if (debug) - std::cout << "energyHCAL:: Energy " << energySum << " from " << kHit << " hits" << std::endl; + edm::LogVerbatim("IsoTrack") << "energyHCAL:: Energy " << energySum << " from " << kHit << " hits"; return energySum; } @@ -400,13 +414,13 @@ namespace spr { khit++; ok = true; if (debug) - std::cout << "energyDetIdHCALneighbours:: Hit " << khit << " " << (HcalDetId)vdets[i] << " energy " - << hit[ihit]->energy() << std::endl; + edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours:: Hit " << khit << " " << (HcalDetId)vdets[i] + << " energy " << hit[ihit]->energy(); } } if (debug) - std::cout << "energyDetIdHCALneighbours::Detector " << i << " " << (HcalDetId)vdets[i] << " energy " << energy - << std::endl; + edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours::Detector " << i << " " << (HcalDetId)vdets[i] + << " energy " << energy; if (ok) { int subdet = ((HcalDetId)(vdets[i].rawId())).subdet(); double eThr = hbThr; @@ -422,7 +436,7 @@ namespace spr { } if (debug) - std::cout << "energyDetIdHCALneighbours::EnergyDetId buffer size " << energyDetId.size() << std::endl; + edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours::EnergyDetId buffer size " << energyDetId.size(); return energyDetId; } diff --git a/Calibration/IsolatedParticles/macros/HBHEMuonOfflineSimAnalyzer.C b/Calibration/IsolatedParticles/macros/HBHEMuonOfflineSimAnalyzer.C index d533a42c0cefc..e2f839f4e0d75 100644 --- a/Calibration/IsolatedParticles/macros/HBHEMuonOfflineSimAnalyzer.C +++ b/Calibration/IsolatedParticles/macros/HBHEMuonOfflineSimAnalyzer.C @@ -1,3 +1,22 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// HBHEMuonOfflineSimAnalyzer h1(infile, outfile, mode, maxDHB, maxDHE); +// h1.Loop() +// +// Offline analysis for MC files +// +// infile const char* Name of the input file +// outfile const char* Name of the output file +// (dyll_PU20_25_output.root) +// 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 (2) +// maxDHB int Maximum number of depths for HB (4) +// maxDHE int Maximum number of depths for HE (7) +// +/////////////////////////////////////////////////////////////////////////////// + #include #include #include @@ -16,178 +35,178 @@ #include class HBHEMuonOfflineSimAnalyzer { - -public : - TTree *fChain; //!pointer to the analyzed TTree or TChain - Int_t fCurrent; - - UInt_t Run_No; - UInt_t Event_No; - UInt_t LumiNumber; - UInt_t BXNumber; - std::vector *pt_of_muon; - std::vector *eta_of_muon; - std::vector *phi_of_muon; - std::vector *p_of_muon; - std::vector *ecal_3x3; - std::vector *ecal_detID; - std::vector *hcal_1x1; - std::vector *matchedId; - std::vector *hcal_detID; - std::vector *hcal_cellHot; - std::vector *activeLength; - std::vector *hcal_edepth1; - std::vector *hcal_edepth2; - std::vector *hcal_edepth3; - std::vector *hcal_edepth4; - std::vector *hcal_activeL1; - std::vector *hcal_activeL2; - std::vector *hcal_activeL3; - std::vector *hcal_activeL4; - std::vector *activeLengthHot; - std::vector *hcal_edepthHot1; - std::vector *hcal_edepthHot2; - std::vector *hcal_edepthHot3; - std::vector *hcal_edepthHot4; - std::vector *hcal_activeHotL1; - std::vector *hcal_activeHotL2; - std::vector *hcal_activeHotL3; - std::vector *hcal_activeHotL4; - std::vector *hcal_edepth5; - std::vector *hcal_activeL5; - std::vector *hcal_edepthHot5; - std::vector *hcal_activeHotL5; - std::vector *hcal_edepth6; - std::vector *hcal_activeL6; - std::vector *hcal_edepthHot6; - std::vector *hcal_activeHotL6; - std::vector *hcal_edepth7; - std::vector *hcal_activeL7; - std::vector *hcal_edepthHot7; - std::vector *hcal_activeHotL7; - - TBranch *b_Run_No; //! - TBranch *b_Event_No; //! - TBranch *b_LumiNumber; //! - TBranch *b_BXNumber; //! - TBranch *b_pt_of_muon; //! - TBranch *b_eta_of_muon; //! - TBranch *b_phi_of_muon; //! - TBranch *b_p_of_muon; //! - TBranch *b_ecal_3x3; //! - TBranch *b_ecal_detID; //! - TBranch *b_hcal_1x1; //! - TBranch *b_hcal_detID; //! - TBranch *b_hcal_cellHot; //! - TBranch *b_activeLength; //! - TBranch *b_hcal_edepth1; //! - TBranch *b_hcal_edepth2; //! - TBranch *b_hcal_edepth3; //! - TBranch *b_hcal_edepth4; //! - TBranch *b_hcal_activeL1; //! - TBranch *b_hcal_activeL2; //! - TBranch *b_hcal_activeL3; //! - TBranch *b_hcal_activeL4; //! - TBranch *b_activeLengthHot; //! - TBranch *b_hcal_edepthHot1; //! - TBranch *b_hcal_edepthHot2; //! - TBranch *b_hcal_edepthHot3; //! - TBranch *b_hcal_edepthHot4; //! - TBranch *b_hcal_activeHotL1; //! - TBranch *b_hcal_activeHotL2; //! - TBranch *b_hcal_activeHotL3; //! - TBranch *b_hcal_activeHotL4; //! - TBranch *b_hcal_edepth5; //! - TBranch *b_hcal_activeL5; //! - TBranch *b_hcal_edepthHot5; //! - TBranch *b_hcal_activeHotL5; //! - TBranch *b_hcal_edepth6; //! - TBranch *b_hcal_activeL6; //! - TBranch *b_hcal_edepthHot6; //! - TBranch *b_hcal_activeHotL6; //! - TBranch *b_hcal_edepth7; //! - TBranch *b_hcal_activeL7; //! - TBranch *b_hcal_edepthHot7; //! - TBranch *b_hcal_activeHotL7; //! - TBranch *b_matchedId; //! - - HBHEMuonOfflineSimAnalyzer(const char *infile, const char *outfile="dyll_PU20_25_output_10.root", const int mode=0, const int maxDHB=4, const int maxDHE=7); - // mode of LHC is kept 1 for 2017 scenario as no change in depth segmentation - // mode of LHC is 0 for 2021 +private: + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; + + UInt_t Run_No; + UInt_t Event_No; + UInt_t LumiNumber; + UInt_t BXNumber; + double pt_of_muon; + double eta_of_muon; + double phi_of_muon; + double p_of_muon; + double ecal_3x3; + unsigned int ecal_detID; + double hcal_1x1; + double matchedId; + unsigned int hcal_detID; + unsigned int hcal_cellHot; + double activeLength; + double hcal_edepth1; + double hcal_edepth2; + double hcal_edepth3; + double hcal_edepth4; + double hcal_activeL1; + double hcal_activeL2; + double hcal_activeL3; + double hcal_activeL4; + double activeLengthHot; + double hcal_edepthHot1; + double hcal_edepthHot2; + double hcal_edepthHot3; + double hcal_edepthHot4; + double hcal_activeHotL1; + double hcal_activeHotL2; + double hcal_activeHotL3; + double hcal_activeHotL4; + double hcal_edepth5; + double hcal_activeL5; + double hcal_edepthHot5; + double hcal_activeHotL5; + double hcal_edepth6; + double hcal_activeL6; + double hcal_edepthHot6; + double hcal_activeHotL6; + double hcal_edepth7; + double hcal_activeL7; + double hcal_edepthHot7; + double hcal_activeHotL7; + + TBranch *b_Run_No; //! + TBranch *b_Event_No; //! + TBranch *b_LumiNumber; //! + TBranch *b_BXNumber; //! + TBranch *b_pt_of_muon; //! + TBranch *b_eta_of_muon; //! + TBranch *b_phi_of_muon; //! + TBranch *b_p_of_muon; //! + TBranch *b_ecal_3x3; //! + TBranch *b_ecal_detID; //! + TBranch *b_hcal_1x1; //! + TBranch *b_hcal_detID; //! + TBranch *b_hcal_cellHot; //! + TBranch *b_activeLength; //! + TBranch *b_hcal_edepth1; //! + TBranch *b_hcal_edepth2; //! + TBranch *b_hcal_edepth3; //! + TBranch *b_hcal_edepth4; //! + TBranch *b_hcal_activeL1; //! + TBranch *b_hcal_activeL2; //! + TBranch *b_hcal_activeL3; //! + TBranch *b_hcal_activeL4; //! + TBranch *b_activeLengthHot; //! + TBranch *b_hcal_edepthHot1; //! + TBranch *b_hcal_edepthHot2; //! + TBranch *b_hcal_edepthHot3; //! + TBranch *b_hcal_edepthHot4; //! + TBranch *b_hcal_activeHotL1; //! + TBranch *b_hcal_activeHotL2; //! + TBranch *b_hcal_activeHotL3; //! + TBranch *b_hcal_activeHotL4; //! + TBranch *b_hcal_edepth5; //! + TBranch *b_hcal_activeL5; //! + TBranch *b_hcal_edepthHot5; //! + TBranch *b_hcal_activeHotL5; //! + TBranch *b_hcal_edepth6; //! + TBranch *b_hcal_activeL6; //! + TBranch *b_hcal_edepthHot6; //! + TBranch *b_hcal_activeHotL6; //! + TBranch *b_hcal_edepth7; //! + TBranch *b_hcal_activeL7; //! + TBranch *b_hcal_edepthHot7; //! + TBranch *b_hcal_activeHotL7; //! + TBranch *b_matchedId; //! + +public: + HBHEMuonOfflineSimAnalyzer(const char *infile, + const char *outfile = "dyll_PU20_25_output.root", + const int mode = 0, + const int maxDHB = 4, + const int maxDHE = 7); virtual ~HBHEMuonOfflineSimAnalyzer(); - virtual Int_t Cut(Long64_t entry); - virtual Int_t GetEntry(Long64_t entry); + 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); + virtual void Init(TTree *tree); + virtual void Loop(); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); std::vector firedTriggers; - void BookHistograms(const char* ); + void BookHistograms(const char *); void WriteHistograms(); - bool LooseMuon(unsigned int ml); - bool tightMuon(unsigned int ml); - bool SoftMuon(unsigned int ml); + bool LooseMuon(); + bool tightMuon(); + bool SoftMuon(); void etaPhiHcal(unsigned int detId, int &eta, int &phi, int &depth); - void etaPhiEcal(unsigned int detId, int& type, int& zside, - int& etaX, int& phiY, int& plane, int& strip); - void calculateP(double pt, double eta, double& pM); + void etaPhiEcal(unsigned int detId, int &type, int &zside, int &etaX, int &phiY, int &plane, int &strip); + void calculateP(double pt, double eta, double &pM); void close(); - int NDepthBins(int ieta); - int NPhiBins(int ieta); - + int NDepthBins(int ieta, int iphi); + int NPhiBins(int ieta); + private: - static const bool debug_=false; - static const int maxDep=7; - static const int maxEta=29; - static const int maxPhi=72; + static const bool debug_ = false; + static const int maxDep = 7; + static const int maxEta = 29; + static const int maxPhi = 72; //3x16x72x2 + 5x4x72x2 + 5x9x36x2 - static const int maxHist=20000;//13032; - int modeLHC, maxDepthHB_, maxDepthHE_, maxDepth_; - int nHist, nDepths[maxEta], nDepthsPhi[maxEta],indxEta[maxEta][maxDep][maxPhi]; + static const int maxHist = 20000; //13032; + int modeLHC_, maxDepthHB_, maxDepthHE_, maxDepth_; + int nHist, nDepths[maxEta], nDepthsPhi[maxEta], indxEta[maxEta][maxDep][maxPhi]; TFile *output_file; - - TH1D *h_Pt_Muon[3], *h_Eta_Muon[3], *h_Phi_Muon[3], *h_P_Muon[3]; - TH1D *h_PF_Muon[3], *h_GlobTrack_Chi[3], *h_Global_Muon_Hits[3]; - TH1D *h_MatchedStations[3], *h_Tight_TransImpactparameter[3]; - TH1D *h_Tight_LongitudinalImpactparameter[3], *h_InnerTrackPixelHits[3]; - TH1D *h_TrackerLayer[3], *h_IsolationR04[3] , *h_Global_Muon[3]; - TH1D *h_LongImpactParameter[3], *h_LongImpactParameterBin1[3], *h_LongImpactParameterBin2[3]; - - TH1D *h_TransImpactParameter[3], *h_TransImpactParameterBin1[3], *h_TransImpactParameterBin2[3]; - TH1D *h_Hot_MuonEnergy_hcal_ClosestCell[3][maxHist] , *h_Hot_MuonEnergy_hcal_HotCell[3][maxHist] , *h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[3][maxHist], *h_HotCell_MuonEnergy_phi[3][maxHist]; - TH2D *h_2D_Bin1[3], *h_2D_Bin2[3]; - TH1D *h_ecal_energy[3], *h_hcal_energy[3], *h_3x3_ecal[3], *h_1x1_hcal[3]; - TH1D *h_MuonHittingEcal[3], *h_HotCell[3], *h_MuonEnergy_hcal[3][maxHist]; - TH1D *h_Hot_MuonEnergy_hcal[3][maxHist]; - TH2D *hcal_ietaVsEnergy[3]; + + TH1D *h_Pt_Muon[3], *h_Eta_Muon[3], *h_Phi_Muon[3], *h_P_Muon[3]; + TH1D *h_PF_Muon[3], *h_GlobTrack_Chi[3], *h_Global_Muon_Hits[3]; + TH1D *h_MatchedStations[3], *h_Tight_TransImpactparameter[3]; + TH1D *h_Tight_LongitudinalImpactparameter[3], *h_InnerTrackPixelHits[3]; + TH1D *h_TrackerLayer[3], *h_IsolationR04[3], *h_Global_Muon[3]; + TH1D *h_LongImpactParameter[3], *h_LongImpactParameterBin1[3], *h_LongImpactParameterBin2[3]; + + TH1D *h_TransImpactParameter[3], *h_TransImpactParameterBin1[3], *h_TransImpactParameterBin2[3]; + TH1D *h_Hot_MuonEnergy_hcal_ClosestCell[3][maxHist], *h_Hot_MuonEnergy_hcal_HotCell[3][maxHist], + *h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[3][maxHist], *h_HotCell_MuonEnergy_phi[3][maxHist]; + TH2D *h_2D_Bin1[3], *h_2D_Bin2[3]; + TH1D *h_ecal_energy[3], *h_hcal_energy[3], *h_3x3_ecal[3], *h_1x1_hcal[3]; + TH1D *h_MuonHittingEcal[3], *h_HotCell[3], *h_MuonEnergy_hcal[3][maxHist]; + TH1D *h_Hot_MuonEnergy_hcal[3][maxHist]; + TH2D *hcal_ietaVsEnergy[3]; TProfile *h_EtaX_hcal[3], *h_PhiY_hcal[3], *h_EtaX_ecal[3], *h_PhiY_ecal[3]; TProfile *h_Eta_ecal[3], *h_Phi_ecal[3]; TProfile *h_MuonEnergy_eta[3][maxDep], *h_MuonEnergy_phi[3][maxDep], *h_MuonEnergy_muon_eta[3][maxDep]; - TProfile *h_Hot_MuonEnergy_eta[3][maxDep], *h_Hot_MuonEnergy_phi[3][maxDep], *h_Hot_MuonEnergy_muon_eta[3][maxDep]; - TProfile *h_IsoHot_MuonEnergy_eta[3][maxDep], *h_IsoHot_MuonEnergy_phi[3][maxDep], *h_IsoHot_MuonEnergy_muon_eta[3][maxDep]; - TProfile *h_IsoWithoutHot_MuonEnergy_eta[3][maxDep], *h_IsoWithoutHot_MuonEnergy_phi[3][maxDep], *h_IsoWithoutHot_MuonEnergy_muon_eta[3][maxDep]; - TProfile *h_HotWithoutIso_MuonEnergy_eta[3][maxDep], *h_HotWithoutIso_MuonEnergy_phi[3][maxDep], *h_HotWithoutIso_MuonEnergy_muon_eta[3][maxDep]; - + TProfile *h_Hot_MuonEnergy_eta[3][maxDep], *h_Hot_MuonEnergy_phi[3][maxDep], *h_Hot_MuonEnergy_muon_eta[3][maxDep]; + TProfile *h_IsoHot_MuonEnergy_eta[3][maxDep], *h_IsoHot_MuonEnergy_phi[3][maxDep], + *h_IsoHot_MuonEnergy_muon_eta[3][maxDep]; + TProfile *h_IsoWithoutHot_MuonEnergy_eta[3][maxDep], *h_IsoWithoutHot_MuonEnergy_phi[3][maxDep], + *h_IsoWithoutHot_MuonEnergy_muon_eta[3][maxDep]; + TProfile *h_HotWithoutIso_MuonEnergy_eta[3][maxDep], *h_HotWithoutIso_MuonEnergy_phi[3][maxDep], + *h_HotWithoutIso_MuonEnergy_muon_eta[3][maxDep]; }; -HBHEMuonOfflineSimAnalyzer::HBHEMuonOfflineSimAnalyzer(const char* infile, - const char* outFileName, - const int mode, - const int maxDHB, - const int maxDHE) { - modeLHC = mode; +HBHEMuonOfflineSimAnalyzer::HBHEMuonOfflineSimAnalyzer( + const char *infile, const char *outFileName, const int mode, const int maxDHB, const int maxDHE) { + modeLHC_ = mode; maxDepthHB_ = maxDHB; maxDepthHE_ = maxDHE; - maxDepth_ = (maxDepthHB_ > maxDepthHE_) ? maxDepthHB_ : maxDepthHE_; + maxDepth_ = (maxDepthHB_ > maxDepthHE_) ? maxDepthHB_ : maxDepthHE_; // if parameter tree is not specified (or zero), connect the file // used to generate this class and read the Tree. - // std::cout<<"maxDepth_"<Get("HcalHBHEMuonAnalyzer"); + TFile *f = new TFile(infile); + TDirectory *dir = (TDirectory *)f->Get("HcalHBHEMuonAnalyzer"); TTree *tree(0); - dir->GetObject("TREE",tree); + dir->GetObject("TREE", tree); Init(tree); //Now book histograms @@ -195,11 +214,12 @@ HBHEMuonOfflineSimAnalyzer::HBHEMuonOfflineSimAnalyzer(const char* infile, } HBHEMuonOfflineSimAnalyzer::~HBHEMuonOfflineSimAnalyzer() { - if (!fChain) return; + if (!fChain) + return; delete fChain->GetCurrentFile(); } -Int_t HBHEMuonOfflineSimAnalyzer::Cut(Long64_t ) { +Int_t HBHEMuonOfflineSimAnalyzer::Cut(Long64_t) { // This function may be called from Loop. // returns 1 if entry is accepted. // returns -1 otherwise. @@ -208,15 +228,18 @@ Int_t HBHEMuonOfflineSimAnalyzer::Cut(Long64_t ) { Int_t HBHEMuonOfflineSimAnalyzer::GetEntry(Long64_t entry) { // Read contents of entry. - if (!fChain) return 0; + if (!fChain) + return 0; return fChain->GetEntry(entry); } Long64_t HBHEMuonOfflineSimAnalyzer::LoadTree(Long64_t entry) { // Set the environment to read one entry - if (!fChain) return -5; + if (!fChain) + return -5; Long64_t centry = fChain->LoadTree(entry); - if (centry < 0) return centry; + if (centry < 0) + return centry; if (fChain->GetTreeNumber() != fCurrent) { fCurrent = fChain->GetTreeNumber(); Notify(); @@ -232,7 +255,7 @@ void HBHEMuonOfflineSimAnalyzer::Init(TTree *tree) { // 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 pt_of_muon = 0; eta_of_muon = 0; @@ -274,7 +297,8 @@ void HBHEMuonOfflineSimAnalyzer::Init(TTree *tree) { hcal_edepthHot7 = 0; hcal_activeHotL7 = 0; matchedId = 0; - if (!tree) return; + if (!tree) + return; fChain = tree; fCurrent = -1; fChain->SetMakeClass(1); @@ -327,130 +351,145 @@ void HBHEMuonOfflineSimAnalyzer::Init(TTree *tree) { } void HBHEMuonOfflineSimAnalyzer::Loop() { - //declarations - if (fChain == 0) return; - + if (fChain == 0) + return; + Long64_t nentries = fChain->GetEntriesFast(); - - if (debug_) std::cout << "nevent = " << nentries << std::endl; + + if (debug_) + std::cout << "nevent = " << nentries << std::endl; Long64_t nbytes = 0, nb = 0; - for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; - - for (unsigned int ml = 0;ml< pt_of_muon->size();ml++) { - - if (debug_) std::cout << "ecal_det_id " << ecal_detID->at(ml) << std::endl; - - if (debug_) std::cout << "hcal_det_id " << std::hex << hcal_detID->at(ml) << std::dec; - - int etaHcal, phiHcal, depthHcal; - etaPhiHcal(hcal_detID->at(ml),etaHcal,phiHcal,depthHcal); - - int eta = (etaHcal > 0) ? etaHcal-1 : -(1+etaHcal); - int nDepth = NDepthBins(eta+1); - int nPhi = NPhiBins(eta+1); - - double phiYHcal = (phiHcal-0.5); - if (debug_) std::cout<<"phiHcal"<Fill(p_of_muon->at(ml)); - 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)); - - double energyFill; - for (int dep=0; dep 36) ? (phiHcal-1) : (phiHcal-1)/2; - double en1(-9999), en2(-9999); - if (dep == 0) { - en1 = hcal_edepth1->at(ml); - en2 = hcal_edepthHot1->at(ml); - energyFill = (hcal_activeHotL1->at(ml) > 0) ? hcal_activeHotL1->at(ml) : 999; - - } else if (dep == 1) { - en1 = hcal_edepth2->at(ml); - en2 = hcal_edepthHot2->at(ml); - energyFill = (hcal_activeHotL2->at(ml) > 0) ? hcal_activeHotL2->at(ml) : 999; - if (debug_) std::cout<<"problem here.. lets see if it got printed"<at(ml); - en2 = hcal_edepthHot3->at(ml); - energyFill = (hcal_activeHotL3->at(ml) > 0) ? hcal_activeHotL3->at(ml) : 999; - } else if (dep == 3) { - en1 = hcal_edepth4->at(ml); - en2 = hcal_edepthHot4->at(ml); - if (debug_) std::cout<<"Hello in 4"<at(ml) > 0) ? hcal_activeHotL4->at(ml) : 999; - } else if (dep == 4) { - en1 = hcal_edepth5->at(ml); - en2 = hcal_edepthHot5->at(ml); - energyFill = (hcal_activeHotL5->at(ml) > 0) ? hcal_activeHotL5->at(ml) : 999; - } else if (dep == 5) { - if (debug_) std::cout << "Energy in depth 6 " << hcal_edepth6->size() << ":" << hcal_edepthHot6->size() << std::endl; - en1 = (hcal_edepth6->size() > ml) ? hcal_edepth6->at(ml) : 0; - en2 = (hcal_edepthHot6->size() >ml) ? hcal_edepthHot6->at(ml) : 0; - energyFill = (hcal_activeHotL6->at(ml) > 0) ? hcal_activeHotL6->at(ml) : 999; - } else if (dep == 6) { - if (debug_) std::cout << "Energy in depth 7 " << hcal_edepth7->size() << ":" << hcal_edepthHot7->size() << std::endl; - en1 = (hcal_edepth7->size() > ml) ? hcal_edepth7->at(ml) : 0; - en2 = (hcal_edepthHot7->size() >ml) ? hcal_edepthHot7->at(ml) : 0; - energyFill = (hcal_activeHotL7->at(ml) > 0) ? hcal_activeHotL7->at(ml) : 999; - } - - if (debug_) std::cout<<" Debug2"< -9999); - bool ok2 = (en2 > -9999); - - if (debug_) std::cout<<"Before Index"< 0) ? indxEta[eta][dep][PHI] : 1+indxEta[eta][dep][PHI]; - - if (debug_) { - std::cout<<"ieta"<at(ml))) continue; - if (ok1) { - if (debug_) std::cout<<"enter ok1"<at(ml)==1) { - if (en2 > 0) { - h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[cut][ind]->Fill(en2/energyFill); - } - if (debug_) std::cout<<"enter hot cell"<at(ml)!=1) { - } - } - - if (debug_) std::cout<<"ETA \t"<GetEntry(jentry); + nbytes += nb; + + if (debug_) { + std::cout << "ecal_det_id " << ecal_detID << std::endl; + std::cout << "hcal_det_id " << std::hex << hcal_detID << std::dec; + } + int etaHcal, phiHcal, depthHcal; + etaPhiHcal(hcal_detID, etaHcal, phiHcal, depthHcal); + + int eta = (etaHcal > 0) ? (etaHcal - 1) : -(1 + etaHcal); + int nDepth = NDepthBins(eta + 1, phiHcal); + int nPhi = NPhiBins(eta + 1); + + double phiYHcal = (phiHcal - 0.5); + if (debug_) + std::cout << "phiHcal" << phiHcal << " phiYHcal" << phiYHcal << std::endl; + + for (int cut = 0; cut < 3; ++cut) { + bool select(false); + if (cut == 0) + select = tightMuon(); + else if (cut == 1) + select = SoftMuon(); + else + select = LooseMuon(); + + if (select) { + // h_P_Muon[cut]->Fill(p_of_muon); + h_P_Muon[cut]->Fill(p_of_muon); + h_Pt_Muon[cut]->Fill(pt_of_muon); + h_Eta_Muon[cut]->Fill(eta_of_muon); + + double energyFill; + for (int dep = 0; dep < nDepth; ++dep) { + if (debug_) { + std::cout << "why on 15/2 only" << std::endl; + std::cout << "dep:" << dep << std::endl; + } + int PHI = (nPhi > 36) ? (phiHcal - 1) : (phiHcal - 1) / 2; + double en1(-9999), en2(-9999); + if (dep == 0) { + en1 = hcal_edepth1; + en2 = hcal_edepthHot1; + energyFill = (hcal_activeHotL1 > 0) ? hcal_activeHotL1 : 999; + } else if (dep == 1) { + en1 = hcal_edepth2; + en2 = hcal_edepthHot2; + energyFill = (hcal_activeHotL2 > 0) ? hcal_activeHotL2 : 999; + if (debug_) + std::cout << "problem here.. lets see if it got printed\n"; + } else if (dep == 2) { + en1 = hcal_edepth3; + en2 = hcal_edepthHot3; + energyFill = (hcal_activeHotL3 > 0) ? hcal_activeHotL3 : 999; + } else if (dep == 3) { + en1 = hcal_edepth4; + en2 = hcal_edepthHot4; + if (debug_) + std::cout << "Hello in 4" << std::endl; + energyFill = (hcal_activeHotL4 > 0) ? hcal_activeHotL4 : 999; + } else if (dep == 4) { + en1 = hcal_edepth5; + en2 = hcal_edepthHot5; + energyFill = (hcal_activeHotL5 > 0) ? hcal_activeHotL5 : 999; + } else if (dep == 5) { + if (debug_) + std::cout << "Energy in depth 6 " << maxDepth_ << ":" << hcal_edepth6 << ":" << hcal_edepthHot6 + << std::endl; + en1 = (maxDepth_ > 5) ? hcal_edepth6 : 0; + en2 = (maxDepth_ > 5) ? hcal_edepthHot6 : 0; + energyFill = (hcal_activeHotL6 > 0) ? hcal_activeHotL6 : 999; + } else if (dep == 6) { + if (debug_) + std::cout << "Energy in depth 7 " << maxDepth_ << ":" << hcal_edepth7 << ":" << hcal_edepthHot7 + << std::endl; + en1 = (maxDepth_ > 6) ? hcal_edepth7 : 0; + en2 = (maxDepth_ > 6) ? hcal_edepthHot7 : 0; + energyFill = (hcal_activeHotL7 > 0) ? hcal_activeHotL7 : 999; + } + + if (debug_) { + std::cout << " Debug2" << std::endl; + std::cout << "ok1" << en1 << std::endl; + std::cout << "ok2" << en2 << std::endl; + } + bool ok1 = (en1 > -9999); + bool ok2 = (en2 > -9999); + + if (debug_) + std::cout << "Before Index" << std::endl; + + int ind = (etaHcal > 0) ? indxEta[eta][dep][PHI] : 1 + indxEta[eta][dep][PHI]; + + if (debug_) { + std::cout << "ieta " << eta << "depth " << dep << "indxEta[eta][dep]:" << indxEta[eta][dep][PHI] + << std::endl; + std::cout << "index showing eta,depth:" << ind << std::endl; + std::cout << "etaHcal: " << etaHcal << " eta " << eta << " dep " << dep << " indx " << ind << std::endl; + } + if (!(matchedId)) + continue; + if (ok1) { + if (debug_) + std::cout << "enter ok1" << std::endl; + if (hcal_cellHot == 1) { + if (en2 > 0) { + h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[cut][ind]->Fill(en2 / energyFill); + } + if (debug_) + std::cout << "enter hot cell" << std::endl; + } + } + + if (ok2) { + if (debug_) + std::cout << "enter ok2" << std::endl; + if (hcal_cellHot != 1) { + } + } + + if (debug_) + std::cout << "ETA \t" << eta << "DEPTH \t" << dep << std::endl; + } } } } @@ -470,234 +509,244 @@ Bool_t HBHEMuonOfflineSimAnalyzer::Notify() { void HBHEMuonOfflineSimAnalyzer::Show(Long64_t entry) { // Print contents of entry. // If entry is not specified, print current entry - if (!fChain) return; + if (!fChain) + return; fChain->Show(entry); } -void HBHEMuonOfflineSimAnalyzer::BookHistograms(const char* fname) { - output_file = TFile::Open(fname,"RECREATE"); +void HBHEMuonOfflineSimAnalyzer::BookHistograms(const char *fname) { + output_file = TFile::Open(fname, "RECREATE"); //output_file->cd(); - std::string type[]={"tight","soft","loose"}; + std::string type[] = {"tight", "soft", "loose"}; char name[128], title[500]; - - std::cout<<"BookHistograms"<cd(); - for (int i=0; i<3; ++i) { - sprintf (name, "h_Pt_Muon_%s", type[i].c_str()); - sprintf (title, "p_{T} of %s muons (GeV)", type[i].c_str()); - h_Pt_Muon[i] = new TH1D(name, title,100,0,200); - - sprintf (name, "h_Eta_Muon_%s", type[i].c_str()); - sprintf (title, "#eta of %s muons", type[i].c_str()); - h_Eta_Muon[i] = new TH1D(name, title,50,-2.5,2.5); - - sprintf (name, "h_Phi_Muon_%s", type[i].c_str()); - sprintf (title, "#phi of %s muons", type[i].c_str()); - h_Phi_Muon[i] = new TH1D(name, title,100,-3.1415926,3.1415926); - - sprintf (name, "h_P_Muon_%s", type[i].c_str()); - sprintf (title, "p of %s muons (GeV)", type[i].c_str()); - h_P_Muon[i] = new TH1D(name, title,100,0,200); - - sprintf (name, "h_PF_Muon_%s", type[i].c_str()); - sprintf (title, "PF %s muons (GeV)", type[i].c_str()); - h_PF_Muon[i] = new TH1D(name, title,2,0,2); - - sprintf (name, "h_Global_Muon_Chi2_%s", type[i].c_str()); - sprintf (title, "Chi2 Global %s muons (GeV)", type[i].c_str()); - h_GlobTrack_Chi[i] = new TH1D(name, title,15,0,15); - - sprintf (name, "h_Global_Muon_Hits_%s", type[i].c_str()); - sprintf (title, "Global Hits %s muons (GeV)", type[i].c_str()); - h_Global_Muon_Hits[i] = new TH1D(name, title,10,0,10) ; - - sprintf (name, "h_Matched_Stations_%s", type[i].c_str()); - sprintf (title, "Matched Stations %s muons (GeV)", type[i].c_str()); - h_MatchedStations[i] = new TH1D(name, title,10,0,10); - - sprintf (name, "h_Transverse_ImpactParameter_%s", type[i].c_str()); - sprintf (title, "Transverse_ImpactParameter of %s muons (GeV)", type[i].c_str()); - h_Tight_TransImpactparameter[i] = new TH1D(name, title,50,0,10); - - sprintf (name, "h_Longitudinal_ImpactParameter_%s", type[i].c_str()); - sprintf (title, "Longitudinal_ImpactParameter of %s muons (GeV)", type[i].c_str()); - h_Tight_LongitudinalImpactparameter[i] = new TH1D(name, title,20,0,10); - - sprintf (name, "h_InnerTrack_PixelHits_%s", type[i].c_str()); - sprintf (title, "InnerTrack_PixelHits of %s muons (GeV)", type[i].c_str()); - h_InnerTrackPixelHits[i]= new TH1D(name, title,20,0,20); - - sprintf (name, "h_TrackLayers_%s", type[i].c_str()); - sprintf (title, "No. of Tracker Layers of %s muons (GeV)", type[i].c_str()); - h_TrackerLayer[i]= new TH1D(name, title,20,0,20);; - - sprintf (name, "h_IsolationR04_%s", type[i].c_str()); - sprintf (title, "IsolationR04 %s muons (GeV)", type[i].c_str()); - h_IsolationR04[i] = new TH1D(name, title,45,0,5);; - - sprintf (name, "h_Global_Muon_%s", type[i].c_str()); - sprintf (title, "Global %s muons (GeV)", type[i].c_str()); - h_Global_Muon[i]= new TH1D(name, title,2,0,2); - - sprintf (name, "h_TransImpactParameter_%s", type[i].c_str()); - sprintf (title, "TransImpactParameter of %s muons (GeV)", type[i].c_str()); - h_TransImpactParameter[i] = new TH1D(name, title,100,0,0.5); - - sprintf (name, "h_TransImpactParameterBin1_%s", type[i].c_str()); - sprintf (title, "TransImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str()); - h_TransImpactParameterBin1[i] = new TH1D(name, title,100,0,0.5); - - sprintf (name, "h_TransImpactParameterBin2_%s", type[i].c_str()); - sprintf (title, "TransImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str()); - h_TransImpactParameterBin2[i] = new TH1D(name, title,100,0,0.5); + for (int i = 0; i < 3; ++i) { + sprintf(name, "h_Pt_Muon_%s", type[i].c_str()); + sprintf(title, "p_{T} of %s muons (GeV)", type[i].c_str()); + h_Pt_Muon[i] = new TH1D(name, title, 100, 0, 200); + + sprintf(name, "h_Eta_Muon_%s", type[i].c_str()); + sprintf(title, "#eta of %s muons", type[i].c_str()); + h_Eta_Muon[i] = new TH1D(name, title, 50, -2.5, 2.5); + + sprintf(name, "h_Phi_Muon_%s", type[i].c_str()); + sprintf(title, "#phi of %s muons", type[i].c_str()); + h_Phi_Muon[i] = new TH1D(name, title, 100, -3.1415926, 3.1415926); + + sprintf(name, "h_P_Muon_%s", type[i].c_str()); + sprintf(title, "p of %s muons (GeV)", type[i].c_str()); + h_P_Muon[i] = new TH1D(name, title, 100, 0, 200); + + sprintf(name, "h_PF_Muon_%s", type[i].c_str()); + sprintf(title, "PF %s muons (GeV)", type[i].c_str()); + h_PF_Muon[i] = new TH1D(name, title, 2, 0, 2); + + sprintf(name, "h_Global_Muon_Chi2_%s", type[i].c_str()); + sprintf(title, "Chi2 Global %s muons (GeV)", type[i].c_str()); + h_GlobTrack_Chi[i] = new TH1D(name, title, 15, 0, 15); + + sprintf(name, "h_Global_Muon_Hits_%s", type[i].c_str()); + sprintf(title, "Global Hits %s muons (GeV)", type[i].c_str()); + h_Global_Muon_Hits[i] = new TH1D(name, title, 10, 0, 10); + + sprintf(name, "h_Matched_Stations_%s", type[i].c_str()); + sprintf(title, "Matched Stations %s muons (GeV)", type[i].c_str()); + h_MatchedStations[i] = new TH1D(name, title, 10, 0, 10); + + sprintf(name, "h_Transverse_ImpactParameter_%s", type[i].c_str()); + sprintf(title, "Transverse_ImpactParameter of %s muons (GeV)", type[i].c_str()); + h_Tight_TransImpactparameter[i] = new TH1D(name, title, 50, 0, 10); + + sprintf(name, "h_Longitudinal_ImpactParameter_%s", type[i].c_str()); + sprintf(title, "Longitudinal_ImpactParameter of %s muons (GeV)", type[i].c_str()); + h_Tight_LongitudinalImpactparameter[i] = new TH1D(name, title, 20, 0, 10); + + sprintf(name, "h_InnerTrack_PixelHits_%s", type[i].c_str()); + sprintf(title, "InnerTrack_PixelHits of %s muons (GeV)", type[i].c_str()); + h_InnerTrackPixelHits[i] = new TH1D(name, title, 20, 0, 20); + + sprintf(name, "h_TrackLayers_%s", type[i].c_str()); + sprintf(title, "No. of Tracker Layers of %s muons (GeV)", type[i].c_str()); + h_TrackerLayer[i] = new TH1D(name, title, 20, 0, 20); + ; + + sprintf(name, "h_IsolationR04_%s", type[i].c_str()); + sprintf(title, "IsolationR04 %s muons (GeV)", type[i].c_str()); + h_IsolationR04[i] = new TH1D(name, title, 45, 0, 5); + ; + + sprintf(name, "h_Global_Muon_%s", type[i].c_str()); + sprintf(title, "Global %s muons (GeV)", type[i].c_str()); + h_Global_Muon[i] = new TH1D(name, title, 2, 0, 2); + + sprintf(name, "h_TransImpactParameter_%s", type[i].c_str()); + sprintf(title, "TransImpactParameter of %s muons (GeV)", type[i].c_str()); + h_TransImpactParameter[i] = new TH1D(name, title, 100, 0, 0.5); + + sprintf(name, "h_TransImpactParameterBin1_%s", type[i].c_str()); + sprintf(title, "TransImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str()); + h_TransImpactParameterBin1[i] = new TH1D(name, title, 100, 0, 0.5); + + sprintf(name, "h_TransImpactParameterBin2_%s", type[i].c_str()); + sprintf(title, "TransImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str()); + h_TransImpactParameterBin2[i] = new TH1D(name, title, 100, 0, 0.5); // - sprintf (name, "h_LongImpactParameter_%s", type[i].c_str()); - sprintf (title, "LongImpactParameter of %s muons (GeV)", type[i].c_str()); - h_LongImpactParameter[i] = new TH1D(name, title,100,0,30); - - sprintf (name, "h_LongImpactParameterBin1_%s", type[i].c_str()); - sprintf (title, "LongImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str()); - h_LongImpactParameterBin1[i] = new TH1D(name, title,100,0,30); - - sprintf (name, "h_LongImpactParameterBin2_%s", type[i].c_str()); - sprintf (title, "LongImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str()); - h_LongImpactParameterBin2[i] = new TH1D(name, title,100,0,30); - - sprintf (name, "h_2D_Bin1_%s", type[i].c_str()); - sprintf (title, "Trans/Long ImpactParameter of %s muons (GeV) in -1.5 <= #phi< 0.5 ", type[i].c_str()); - h_2D_Bin1[i] = new TH2D(name, title, 100,0,0.5,100,0,30); - - sprintf (name, "h_2D_Bin2_%s", type[i].c_str()); - sprintf (title, "Trans/Long ImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str()); - h_2D_Bin2[i] = new TH2D(name, title, 100,0,0.5,100,0,30); - - sprintf (name, "h_ecal_energy_%s", type[i].c_str()); - sprintf (title, "ECAL energy for %s muons", type[i].c_str()); - h_ecal_energy[i] = new TH1D(name, title,1000,-10.0,90.0); - - sprintf (name, "h_hcal_energy_%s", type[i].c_str()); - sprintf (title, "HCAL energy for %s muons", type[i].c_str()); - h_hcal_energy[i] = new TH1D(name, title,500,-10.0,90.0); - - sprintf (name, "h_3x3_ecal_%s", type[i].c_str()); - sprintf (title, "ECAL energy in 3x3 for %s muons", type[i].c_str()); - h_3x3_ecal[i] = new TH1D(name, title,1000,-10.0,90.0); - - sprintf (name, "h_1x1_hcal_%s", type[i].c_str()); - sprintf (title, "HCAL energy in 1x1 for %s muons", type[i].c_str()); - h_1x1_hcal[i] = new TH1D(name, title,500,-10.0,90.0); - - sprintf (name, "h_EtaX_hcal_%s", type[i].c_str()); - sprintf (title, "HCAL energy as a function of i#eta for %s muons", type[i].c_str()); - h_EtaX_hcal[i] = new TProfile(name, title,60,-30.0,30.0); - - sprintf (name, "h_PhiY_hcal_%s", type[i].c_str()); - sprintf (title, "HCAL energy as a function of i#phi for %s muons", type[i].c_str()); - h_PhiY_hcal[i] = new TProfile(name, title,72,0,72); - - sprintf (name, "h_EtaX_ecal_%s", type[i].c_str()); - sprintf (title, "EB energy as a function of i#eta for %s muons", type[i].c_str()); - h_EtaX_ecal[i] = new TProfile(name, title,170,-85.0,85.0); - - sprintf (name, "h_PhiY_ecal_%s", type[i].c_str()); - sprintf (title, "EB energy as a function of i#phi for %s muons", type[i].c_str()); - h_PhiY_ecal[i] = new TProfile(name, title,360,0,360); - - sprintf (name, "h_Eta_ecal_%s", type[i].c_str()); - sprintf (title, "ECAL energy as a function of #eta for %s muons", type[i].c_str()); - h_Eta_ecal[i] = new TProfile(name, title,100,-2.5,2.5); - - sprintf (name, "h_Phi_ecal_%s", type[i].c_str()); - sprintf (title, "ECAL energy as a function of #phi for %s muons", type[i].c_str()); - h_Phi_ecal[i] = new TProfile(name, title,100,-3.1415926,3.1415926); - - sprintf (name, "h_MuonHittingEcal_%s", type[i].c_str()); - sprintf (title, "%s muons hitting ECAL", type[i].c_str()); - h_MuonHittingEcal[i] = new TH1D(name, title,100,0,5.0); - - sprintf (name, "h_HotCell_%s", type[i].c_str()); - sprintf (title, "Hot cell for %s muons", type[i].c_str()); - h_HotCell[i] = new TH1D(name, title,100,0,2); - - std::cout<<"problem here"<cd(); - for (int eta=0; eta<29; ++eta) { - int nDepth = NDepthBins(eta+1); - int nPhi = NPhiBins(eta+1); - //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); - //d_output_file[i][eta]->cd(); - for (int depth=0; depthSumw2(); - - ih++; - 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()); - h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih] = new TH1D(name, title,4000,0.0,1.0); - h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Sumw2(); - } + sprintf(name, "h_LongImpactParameter_%s", type[i].c_str()); + sprintf(title, "LongImpactParameter of %s muons (GeV)", type[i].c_str()); + h_LongImpactParameter[i] = new TH1D(name, title, 100, 0, 30); + + sprintf(name, "h_LongImpactParameterBin1_%s", type[i].c_str()); + sprintf(title, "LongImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str()); + h_LongImpactParameterBin1[i] = new TH1D(name, title, 100, 0, 30); + + sprintf(name, "h_LongImpactParameterBin2_%s", type[i].c_str()); + sprintf(title, "LongImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str()); + h_LongImpactParameterBin2[i] = new TH1D(name, title, 100, 0, 30); + + sprintf(name, "h_2D_Bin1_%s", type[i].c_str()); + sprintf(title, "Trans/Long ImpactParameter of %s muons (GeV) in -1.5 <= #phi< 0.5 ", type[i].c_str()); + h_2D_Bin1[i] = new TH2D(name, title, 100, 0, 0.5, 100, 0, 30); + + sprintf(name, "h_2D_Bin2_%s", type[i].c_str()); + sprintf(title, "Trans/Long ImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str()); + h_2D_Bin2[i] = new TH2D(name, title, 100, 0, 0.5, 100, 0, 30); + + sprintf(name, "h_ecal_energy_%s", type[i].c_str()); + sprintf(title, "ECAL energy for %s muons", type[i].c_str()); + h_ecal_energy[i] = new TH1D(name, title, 1000, -10.0, 90.0); + + sprintf(name, "h_hcal_energy_%s", type[i].c_str()); + sprintf(title, "HCAL energy for %s muons", type[i].c_str()); + h_hcal_energy[i] = new TH1D(name, title, 500, -10.0, 90.0); + + sprintf(name, "h_3x3_ecal_%s", type[i].c_str()); + sprintf(title, "ECAL energy in 3x3 for %s muons", type[i].c_str()); + h_3x3_ecal[i] = new TH1D(name, title, 1000, -10.0, 90.0); + + sprintf(name, "h_1x1_hcal_%s", type[i].c_str()); + sprintf(title, "HCAL energy in 1x1 for %s muons", type[i].c_str()); + h_1x1_hcal[i] = new TH1D(name, title, 500, -10.0, 90.0); + + sprintf(name, "h_EtaX_hcal_%s", type[i].c_str()); + sprintf(title, "HCAL energy as a function of i#eta for %s muons", type[i].c_str()); + h_EtaX_hcal[i] = new TProfile(name, title, 60, -30.0, 30.0); + + sprintf(name, "h_PhiY_hcal_%s", type[i].c_str()); + sprintf(title, "HCAL energy as a function of i#phi for %s muons", type[i].c_str()); + h_PhiY_hcal[i] = new TProfile(name, title, 72, 0, 72); + + sprintf(name, "h_EtaX_ecal_%s", type[i].c_str()); + sprintf(title, "EB energy as a function of i#eta for %s muons", type[i].c_str()); + h_EtaX_ecal[i] = new TProfile(name, title, 170, -85.0, 85.0); + + sprintf(name, "h_PhiY_ecal_%s", type[i].c_str()); + sprintf(title, "EB energy as a function of i#phi for %s muons", type[i].c_str()); + h_PhiY_ecal[i] = new TProfile(name, title, 360, 0, 360); + + sprintf(name, "h_Eta_ecal_%s", type[i].c_str()); + sprintf(title, "ECAL energy as a function of #eta for %s muons", type[i].c_str()); + h_Eta_ecal[i] = new TProfile(name, title, 100, -2.5, 2.5); + + sprintf(name, "h_Phi_ecal_%s", type[i].c_str()); + sprintf(title, "ECAL energy as a function of #phi for %s muons", type[i].c_str()); + h_Phi_ecal[i] = new TProfile(name, title, 100, -3.1415926, 3.1415926); + + sprintf(name, "h_MuonHittingEcal_%s", type[i].c_str()); + sprintf(title, "%s muons hitting ECAL", type[i].c_str()); + h_MuonHittingEcal[i] = new TH1D(name, title, 100, 0, 5.0); + + sprintf(name, "h_HotCell_%s", type[i].c_str()); + sprintf(title, "Hot cell for %s muons", type[i].c_str()); + h_HotCell[i] = new TH1D(name, title, 100, 0, 2); + + std::cout << "problem here" << std::endl; + for (int eta = 0; eta < 29; ++eta) { + int nDepth = NDepthBins(eta + 1, -1); + int nPhi = NPhiBins(eta + 1); + for (int depth = 0; depth < nDepth; ++depth) { + for (int PHI = 0; PHI < nPhi; ++PHI) { + int PHI0 = (nPhi == 72) ? PHI + 1 : 2 * PHI + 1; + int ih = indxEta[eta][depth][PHI]; + std::cout << "eta:" << eta << " depth:" << depth << " PHI:" << PHI << ":" << PHI0 << " ih:" << ih + << std::endl; + + 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()); + h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih] = new TH1D(name, title, 4000, 0.0, 1.0); + h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Sumw2(); + + ih++; + 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()); + h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih] = new TH1D(name, title, 4000, 0.0, 1.0); + h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Sumw2(); + } } //output_file->cd(); - } + } } //output_file->cd(); } -bool HBHEMuonOfflineSimAnalyzer::LooseMuon(unsigned int ml){ - if (pt_of_muon->at(ml) > 20.){ - if (fabs(eta_of_muon->at(ml)) <= 5.0) { - +bool HBHEMuonOfflineSimAnalyzer::LooseMuon() { + if (pt_of_muon > 20.) { + if (fabs(eta_of_muon) <= 5.0) { return true; } } - + return false; -} +} -bool HBHEMuonOfflineSimAnalyzer::SoftMuon(unsigned int ml){ - if (pt_of_muon->at(ml) > 20.){ - if (fabs(eta_of_muon->at(ml)) <= 5.0) { - +bool HBHEMuonOfflineSimAnalyzer::SoftMuon() { + if (pt_of_muon > 20.) { + if (fabs(eta_of_muon) <= 5.0) { return true; } } - + return false; } -bool HBHEMuonOfflineSimAnalyzer::tightMuon(unsigned int ml) { - if (pt_of_muon->at(ml) > 20.){ - if (fabs(eta_of_muon->at(ml)) <= 5.0) { - +bool HBHEMuonOfflineSimAnalyzer::tightMuon() { + if (pt_of_muon > 20.) { + if (fabs(eta_of_muon) <= 5.0) { return true; } } @@ -707,52 +756,48 @@ bool HBHEMuonOfflineSimAnalyzer::tightMuon(unsigned int ml) { void HBHEMuonOfflineSimAnalyzer::etaPhiHcal(unsigned int detId, int &eta, int &phi, int &depth) { int zside, etaAbs; - if ((detId&0x1000000)==0) { - zside = (detId&0x2000)?(1):(-1); - etaAbs = (detId>>7)&0x3F; - phi = detId&0x7F; - depth = (detId>>14)&0x1F; + if ((detId & 0x1000000) == 0) { + zside = (detId & 0x2000) ? (1) : (-1); + etaAbs = (detId >> 7) & 0x3F; + phi = detId & 0x7F; + depth = (detId >> 14) & 0x1F; } else { - zside = (detId&0x80000)?(1):(-1); - etaAbs = (detId>>10)&0x1FF; - phi = detId&0x3FF; - depth = (detId>>20)&0xF; + zside = (detId & 0x80000) ? (1) : (-1); + etaAbs = (detId >> 10) & 0x1FF; + phi = detId & 0x3FF; + depth = (detId >> 20) & 0xF; } - eta = etaAbs*zside; + eta = etaAbs * zside; } -void HBHEMuonOfflineSimAnalyzer::etaPhiEcal(unsigned int detId, int& type, - int& zside, int& etaX, int& phiY, - int& plane, int& strip) { - - type = ((detId>>25)&0x7); - // std::cout<<"type"<> 25) & 0x7); plane = strip = 0; - if (type==1) { + if (type == 1) { //Ecal Barrel - zside = (detId&0x10000)?(1):(-1); - etaX = (detId>>9)&0x7F; - phiY = detId&0x1FF; - } else if (type==2) { - zside = (detId&0x4000)?(1):(-1); - etaX = (detId>>7)&0x7F; - phiY = (detId&0x7F); - } else if (type==3) { - zside = (detId&0x80000)?(1):(-1); - etaX = (detId>>6)&0x3F; + zside = (detId & 0x10000) ? (1) : (-1); + etaX = (detId >> 9) & 0x7F; + phiY = detId & 0x1FF; + } else if (type == 2) { + zside = (detId & 0x4000) ? (1) : (-1); + etaX = (detId >> 7) & 0x7F; + phiY = (detId & 0x7F); + } else if (type == 3) { + zside = (detId & 0x80000) ? (1) : (-1); + etaX = (detId >> 6) & 0x3F; /** get the sensor iy */ - phiY = (detId>>12)&0x3F; + phiY = (detId >> 12) & 0x3F; /** get the strip */ - plane = ((detId>>18)&0x1)+1; - strip = detId&0x3F; + plane = ((detId >> 18) & 0x1) + 1; + strip = detId & 0x3F; } else { zside = etaX = phiY = 0; } } - -void HBHEMuonOfflineSimAnalyzer::calculateP(double pt, double eta, double& pM) { - pM = (pt*cos(2*(1/atan(exp(eta))))); +void HBHEMuonOfflineSimAnalyzer::calculateP(double pt, double eta, double &pM) { + pM = (pt * cos(2 * (1 / atan(exp(eta))))); } void HBHEMuonOfflineSimAnalyzer::close() { @@ -765,33 +810,29 @@ void HBHEMuonOfflineSimAnalyzer::close() { std::cout << "now doing return" << std::endl; } - void HBHEMuonOfflineSimAnalyzer::WriteHistograms() { - - //output_file->cd(); - std::string type[]={"tight","soft","loose"}; + std::string type[] = {"tight", "soft", "loose"}; char name[128]; - - std::cout<<"BookHistograms"<cd(); - for (int i=0; i<3; ++i) { - + for (int i = 0; i < 3; ++i) { h_Pt_Muon[i]->Write(); h_Eta_Muon[i]->Write(); h_Phi_Muon[i]->Write(); @@ -802,7 +843,7 @@ void HBHEMuonOfflineSimAnalyzer::WriteHistograms() { h_Global_Muon_Hits[i]->Write(); h_MatchedStations[i]->Write(); - h_Tight_TransImpactparameter[i]->Write(); + h_Tight_TransImpactparameter[i]->Write(); h_Tight_LongitudinalImpactparameter[i]->Write(); h_InnerTrackPixelHits[i]->Write(); @@ -810,7 +851,8 @@ void HBHEMuonOfflineSimAnalyzer::WriteHistograms() { h_IsolationR04[i]->Write(); h_Global_Muon[i]->Write(); - h_TransImpactParameter[i]->Write();; + h_TransImpactParameter[i]->Write(); + ; h_TransImpactParameterBin1[i]->Write(); h_TransImpactParameterBin2[i]->Write(); // @@ -819,62 +861,102 @@ void HBHEMuonOfflineSimAnalyzer::WriteHistograms() { h_LongImpactParameterBin2[i]->Write(); h_ecal_energy[i]->Write(); - h_hcal_energy[i]->Write();; + h_hcal_energy[i]->Write(); + ; h_3x3_ecal[i]->Write(); - h_1x1_hcal[i]->Write();; + h_1x1_hcal[i]->Write(); + ; h_EtaX_hcal[i]->Write(); - h_PhiY_hcal[i]->Write();; - - h_EtaX_ecal[i]->Write();; - h_PhiY_ecal[i]->Write();; - h_Eta_ecal[i]->Write();; - h_Phi_ecal[i]->Write();; - h_MuonHittingEcal[i]->Write();; - h_HotCell[i]->Write();; - + h_PhiY_hcal[i]->Write(); + ; + + h_EtaX_ecal[i]->Write(); + ; + h_PhiY_ecal[i]->Write(); + ; + h_Eta_ecal[i]->Write(); + ; + h_Phi_ecal[i]->Write(); + ; + h_MuonHittingEcal[i]->Write(); + ; + h_HotCell[i]->Write(); + ; + output_file->cd(); - for (int eta=0; eta<29; ++eta) { - int nDepth = NDepthBins(eta+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); - //output_file->cd(name); + 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); d_output_file[i][eta]->cd(); - for (int depth=0; depthWrite(); - - ih++; - h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Write(); - } + for (int depth = 0; depth < nDepth; ++depth) { + for (int PHI = 0; PHI < nPhi; ++PHI) { + int ih = indxEta[eta][depth][PHI]; + h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Write(); + ih++; + h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Write(); + } } output_file->cd(); - } + } } output_file->cd(); } -int HBHEMuonOfflineSimAnalyzer::NDepthBins(int eta) { - // int nDepth[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}; - // for 2021 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}; - int nDepth[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}; - - int nbin = nDepth[eta-1]; - if (modeLHC == 0) { - if (eta<=15) { +int HBHEMuonOfflineSimAnalyzer::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}; + // for a test scenario with multi depth segmentation considered during Run 1 + // 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) + // = 5 --> corresponds to Run 4 (post LS3) + int nbin(0); + if (modeLHC_ == 0) { + if (eta <= 15) { nbin = maxDepthHB_; } else if (eta == 16) { nbin = 4; } else { nbin = maxDepthHE_; } + } else if (modeLHC_ == 1) { + nbin = nDepthR1[eta - 1]; + } else if (modeLHC_ == 2) { + nbin = nDepthR2[eta - 1]; + } else if (modeLHC_ == 3) { + nbin = nDepthR3[eta - 1]; + } else if (modeLHC_ == 4) { + 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; } diff --git a/Calibration/IsolatedParticles/macros/PlotTracks.C b/Calibration/IsolatedParticles/macros/PlotTracks.C index fd7eeae3fcebf..ac2d503f4b920 100644 --- a/Calibration/IsolatedParticles/macros/PlotTracks.C +++ b/Calibration/IsolatedParticles/macros/PlotTracks.C @@ -16,10 +16,10 @@ // logy = If y-axis scale shuld by linear/logarithmic [true] // pos = position of the statistics boxes [0] // pv = flag deciding call to plotEnergyPV vs plotEnergy [false] -// savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1] +// savePlot= Plot to be saved: no(-1), eps(0), gif(1), pdf(2), C(3) [-1] // // void plotEnergyAll(std::string fname, std::string hlt, int pv, int data, -// bool varbin, int rebin, bool approve, bool logy, int pos, +// bool varbin, int rebin, bool approve, bool logy, int pos, // int var, int ene, int eta, int savePlot) // Plots energy distributions for a number of variables, eta and energy bins // @@ -36,9 +36,9 @@ // var = variable name (-1 means all variables) [-1] // ene = energy bin (-1 means all energy bins) [-1] // eta = eta bin (-1 means all eta bins) [-1] -// savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1] +// savePlot= Plot to be saved: no(-1), eps(0), gif(1), pdf(2), C(3) [-1] // -// void plotEMeanAll(int data, int models, bool ratio, bool approve, +// void plotEMeanAll(int data, int models, bool ratio, bool approve, // std::string postfix, int savePlot) // // Plots mean energy response as a function of track momentum or the ratio @@ -49,14 +49,14 @@ // wrt a reference // approve= If meant for approval talks [true] // postfix= String to be added in the name of saved file [""] -// savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1] +// savePlot= Plot to be saved: no(-1), eps(0), gif(1), pdf(2), C(3) [-1] // // void plotEMean(std::string fname, std::string hlt, int models, int var, // int eta, int pv, int data, bool ratio, bool approve, // std::string postfix, int savePlot) // // Plots mean energy response as a function of track momentum or the ratio -// MC/Data +// MC/Data // // fname = name of the i/p root file [""] // HLT = type of HLT used (to be given in the figure legend) ["All HLTs"] @@ -72,9 +72,9 @@ // a reference // approve= If meant for approval talks [true] // postfix= String to be added in the name of saved file [""] -// savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1] +// savePlot= Plot to be saved: no(-1), eps(0), gif(1), pdf(2), C(3) [-1] // -// void plotEnergy(std::string fname, std::string HLT, int var, int ien, +// void plotEnergy(std::string fname, std::string HLT, int var, int ien, // int eta, bool varbin, int rebin, bool approve, bool logy, // int pos, int coloff) // @@ -120,10 +120,10 @@ // logy = If y-axis scale shuld by linear/logarithmic [true] // pos = position of the statistics boxes [0] // -// void plotTrack(std::string fname, std::string HLT, int var, bool varbin, +// void plotTrack(std::string fname, std::string HLT, int var, bool varbin, // int rebin, bool approve, bool logy, int pos) // -// Plots kinematic propeties of the track +// Plots kinematic propeties of the track // // fname = name of the i/p root file ["hlt.root"] // HLT = type of HLT used (to be given in the figure legend ["All HLTs"] @@ -163,6 +163,14 @@ // approve= If meant for approval talks [false] // logy = If y-axis scale shuld by linear/logarithmic [true] // pos = position of the statistics boxes [0] +// +// void plotEMeanRatioXAll(std::string fname=, std::string hlt, std::string postfix="F", int savePlot=2); +// +// fname = name of the i/p root file ["pikp/FBE4p3vMixStudyHLT.root"] +// hlt = name of the tag for the file ["10.4p03 FTFP_BERT_EMM"] +// postfix = string fixed to the saved file ["F"] +// savePlot = flag to save the canvas [2] +// //////////////////////////////////////////////////////////////////////////////// #include "TCanvas.h" @@ -171,12 +179,14 @@ #include "TFile.h" #include "TFitResult.h" #include "TGraph.h" +#include "TGraphErrors.h" #include "TGraphAsymmErrors.h" #include "TH1D.h" #include "TH2D.h" #include "THStack.h" #include "TLegend.h" #include "TMath.h" +#include "TPolyLine.h" #include "TProfile.h" #include "TPaveStats.h" #include "TPaveText.h" @@ -189,71 +199,71 @@ #include #include +//const int nmodelm=1, nmodelx=2, nmodels=3; //const int nmodelm=15, nmodelx=16, nmodels=2; -const int nmodelm=3, nmodelx=4, nmodels=2; +//const int nmodelm=3, nmodelx=4, nmodels=3; +//const int nmodelm=4, nmodelx=5, nmodels=3; //const int nmodelm=4, nmodelx=5, nmodels=2; -//const int nmodelm=2, nmodelx=3, nmodels=2; -int styles[7] = {20, 21, 24, 22, 23, 33, 25}; -int colors[7] = {1, 2, 4, 6, 7, 38, 3}; -std::string names[7] = {"All", "Quality", "okEcal", "EcalCharIso", - "HcalCharIso", "EcalNeutIso", "HcalNeutIso"}; -std::string namefull[7]= {"All tracks", "Good quality tracks", - "Tracks reaching ECAL", "Charge isolation in ECAL", - "Charge isolation in HCAL", "Isolated in ECAL", - "Isolated in HCAL"}; -std::string nameEta[4] = {"i#eta 1:6","i#eta 7:12","i#eta 13:16","i#eta 17:22"}; -std::string nameEtas[4]= {"|#eta| < 0.52", "0.52 < |#eta| < 1.04", - "1.04 < |#eta| < 1.39", "1.39 < |#eta| < 2.01"}; -std::string namePV[5] = {"all PV","PV 1:1","PV 2:2","PV 3:5","PV > 5"}; +//const int nmodelm=2, nmodelx=3, nmodels=3; +const int nmodelm = 5, nmodelx = 6, nmodels = 3; +const int etaMax = 3; +//int styles[7] = {20, 23, 22, 24, 25, 21, 33}; +//int colors[7] = {1, 2, 4, 6, 7, 46, 38}; +int styles[7] = {20, 23, 21, 22, 24, 25, 33}; +int colors[7] = {1, 2, 6, 4, 3, 7, 38}; +int lstyle[7] = {1, 2, 3, 4, 5, 6, 7}; +std::string names[7] = {"All", "Quality", "okEcal", "EcalCharIso", "HcalCharIso", "EcalNeutIso", "HcalNeutIso"}; +std::string namefull[7] = {"All tracks", + "Good quality tracks", + "Tracks reaching ECAL", + "Charge isolation in ECAL", + "Charge isolation in HCAL", + "Isolated in ECAL", + "Isolated in HCAL"}; +std::string nameEta[4] = {"i#eta 1:6", "i#eta 7:12", "i#eta 13:16", "i#eta 17:22"}; +std::string nameEtas[4] = {"|#eta| < 0.52", "0.52 < |#eta| < 1.04", "1.04 < |#eta| < 1.39", "1.39 < |#eta| < 2.01"}; +std::string namePV[5] = {"all PV", "PV 1:1", "PV 2:2", "PV 3:5", "PV > 5"}; std::string varname[4] = {"p", "pt", "eta", "phi"}; -std::string vartitle[4]= {"p (GeV/c)", "p_{T} (GeV/c)", "#eta", "#phi"}; -std::string nameC[2] = {"Ecal", "Hcal"}; -std::string nameCF[2] = {"ECAL", "HCAL"}; -std::string varnameC[4]= {"maxNearP", "ediff", "ene1", "ene2"}; -std::string vartitlC[4]= {"Charge isolation energy", - "Neutral isolation energy", - "Energy in smaller cone", - "Energy in larger cone"}; -const int NPT=10; -double mom[NPT]={1.5,2.5,3.5,4.5,5.5,6.5,8.0,10.0,13.0,17.5}; -double dmom[NPT]={0.5,0.5,0.5,0.5,0.5,0.5,1.0,1.0,2.0,2.5}; -std::string varPs[NPT] = {"1:2","2:3","3:4","4:5","5:6","6:7","7:9", - "9:11","11:15","15:20"}; -std::string varPs1[NPT] = {"1","2","3","4","5","6","7","9","11","15"}; -std::string varPPs[NPT] = {"1-2 GeV","2-3 GeV","3-4 GeV","4-5 GeV","5-6 GeV", - "6-7 GeV","7-9 GeV","9-11 GeV","11-15 GeV", - - "15-20 GeV"}; +std::string vartitle[4] = {"p (GeV/c)", "p_{T} (GeV/c)", "#eta", "#phi"}; +std::string nameC[2] = {"Ecal", "Hcal"}; +std::string nameCF[2] = {"ECAL", "HCAL"}; +std::string varnameC[4] = {"maxNearP", "ediff", "ene1", "ene2"}; +std::string vartitlC[4] = { + "Charge isolation energy", "Neutral isolation energy", "Energy in smaller cone", "Energy in larger cone"}; +std::string cmsP = "CMS Preliminary"; +//std::string cmsP = "CMS"; +std::string fileData = "AllDataStudyHLT.root"; +const int NPT = 10; +double mom[NPT] = {1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 8.0, 10.0, 13.0, 17.5}; +double dmom[NPT] = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0, 2.0, 2.5}; +std::string varPs[NPT] = {"1:2", "2:3", "3:4", "4:5", "5:6", "6:7", "7:9", "9:11", "11:15", "15:20"}; +std::string varPs1[NPT] = {"1", "2", "3", "4", "5", "6", "7", "9", "11", "15"}; +std::string varPPs[NPT] = { + "1-2 GeV", "2-3 GeV", "3-4 GeV", "4-5 GeV", "5-6 GeV", "6-7 GeV", "7-9 GeV", "9-11 GeV", "11-15 GeV", "15-20 GeV"}; std::string varEta[4] = {"1:6", "7:12", "13:16", "17:23"}; -std::string varEta1[4]= {"1", "2", "3", "4"}; -std::string varEne[6] = {"E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})", - "E_{11x11}", "H_{5x5}", "(E_{11x11}+H_{5x5})"}; -std::string varEne1[6]= {"E7x7","H3x3","E7x7H_3x3","E11x11","H5x5", "E11x11H5x5"}; -const int nbins=100; -double xbins[nbins+1] = {0.00,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09, - 0.10,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19, - 0.20,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29, - 0.30,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39, - 0.40,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49, - 0.50,0.52,0.54,0.56,0.58,0.60,0.62,0.64,0.66,0.68, - 0.70,0.72,0.74,0.76,0.78,0.80,0.82,0.84,0.86,0.88, - 0.90,0.92,0.94,0.96,0.98,1.00,1.05,1.10,1.15,1.20, - 1.25,1.30,1.35,1.40,1.45,1.50,1.55,1.60,1.65,1.70, - 1.75,1.80,1.85,1.90,1.95,2.00,2.10,2.20,2.30,2.40, - 2.50}; -int ibins[nbins+1] = { 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, - 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, - 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, - 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, - 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, - 61, 63, 65, 67, 69, 71, 73, 75, 77, 79, - 81, 83, 85, 87, 89, 91, 93, 95, 97, 99, - 101,103,105,107,109,111,116,121,126,131, - 136,141,146,151,156,161,166,171,176,181, - 186,191,196,201,206,211,221,231,241,251,261}; - +std::string varEta1[4] = {"1", "2", "3", "4"}; +std::string varEne[6] = {"E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})", "E_{11x11}", "H_{5x5}", "(E_{11x11}+H_{5x5})"}; +std::string varEne1[6] = {"E7x7", "H3x3", "E7x7H3x3", "E11x11", "H5x5", "E11x11H5x5"}; +const int nbins = 100; +double xbins[nbins + 1] = {0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, + 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, + 0.30, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.40, 0.41, 0.42, 0.43, 0.44, + 0.45, 0.46, 0.47, 0.48, 0.49, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64, 0.66, 0.68, + 0.70, 0.72, 0.74, 0.76, 0.78, 0.80, 0.82, 0.84, 0.86, 0.88, 0.90, 0.92, 0.94, 0.96, 0.98, + 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, 1.50, 1.55, 1.60, 1.65, 1.70, + 1.75, 1.80, 1.85, 1.90, 1.95, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50}; +int ibins[nbins + 1] = {11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, + 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, + 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, + 63, 65, 67, 69, 71, 73, 75, 77, 79, 81, 83, 85, 87, 89, 91, 93, 95, + 97, 99, 101, 103, 105, 107, 109, 111, 116, 121, 126, 131, 136, 141, 146, 151, 156, + 161, 166, 171, 176, 181, 186, 191, 196, 201, 206, 211, 221, 231, 241, 251, 261}; +/* std::string files[nmodels]={"ZeroBiasStudyHLT.root","MinimumBiasStudyHLT.root"}; std::string types[nmodels]={"Zero Bias Data","Minimum Bias Data"}; +*/ +std::string files[nmodels] = {"2017C_ZB.root", "2017G_ZB.root", "2017H_ZB.root"}; +std::string types[nmodels] = {"Zero Bias (2017C)", "Zero Bias (2017G)", "Zero Bias (2017H)"}; /* std::string filem[nmodelm]={"pikp/FBE2p2StudyHLT.root", @@ -335,12 +345,166 @@ std::string typem[nmodelm]={"10.3.ref06 FTFP_BERT_EMM (Native)", "10.3.ref06 FTFP_BERT_EMM (VecGeom 4)", "10.3.ref06 FTFP_BERT_EMM (VecGeom Corr)"}; */ +/* std::string filem[nmodelm]={"pikp/FBE4r00vMixStudyHLT.root", "pikp/FBE4r01MixStudyHLT.root", "pikp/FBE4r01vMixStudyHLT.root"}; std::string typem[nmodelm]={"10.4 VecGeom FTFP_BERT_EMM", "10.4.ref01 FTFP_BERT_EMM", "10.4.ref01 VecGeom FTFP_BERT_EMM"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE4MixStudyHLT10d.root", + "pikp/FBE4MixStudyHLT10s.root", + "pikp/FBE4MixStudyHLT10t.root", + "pikp/FBE4MixStudyHLT10st.root"}; +std::string typem[nmodelm]={"10.4 FTFP_BERT_EMM (Default)", + "10.4 FTFP_BERT_EMM (New StepIntegrator)", + "10.4 FTFP_BERT_EMM (Smart Track)", + "10.4 FTFP_BERT_EMM (New Stepper + Smart Track)"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE2p2StudyHLT.root", + "pikp/FBE4MixStudyHLT1025.root", + "pikp/FBE5c01vMixStudyHLT.root", + "pikp/FBE5c02vMixStudyHLT.root", + "pikp/QFB5c02vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.2.p02 FTFP_BERT_EMM", + "10.4 FTFP_BERT_EMM", + "10.5.cand01 FTFP_BERT_EMM", + "10.5.cand02 FTFP_BERT_EMM", + "10.5.cand02 QGSP_FTFP_BERT_EML"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE5r00vMixStudyHLT.root", + "pikp/FBE5r05vMixStudyHLT.root", + "pikp/FBE5r07vMixStudyHLT.root", + "pikp/FBE5r08vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4.p03 FTFP_BERT_EMM", + "10.5 FTFP_BERT_EMM", + "10.5.ref05 FTFP_BERT_EMM", + "10.5.ref07 FTFP_BERT_EMM", + "10.5.ref08 FTFP_BERT_EMM"}; +*/ +/* +std::string filem[nmodelm]={"pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB5r00vMixStudyHLT.root", + "pikp/QFB5r05vMixStudyHLT.root", + "pikp/QFB5r07vMixStudyHLT.root", + "pikp/QFB5r08vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4.p03 QGSP_FTFP_BERT_EML", + "10.5 QGSP_FTFP_BERT_EML", + "10.5.ref05 QGSP_FTFP_BERT_EML", + "10.5.ref07 QGSP_FTFP_BERT_EML", + "10.5.ref08 QGSP_FTFP_BERT_EML"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE2p2StudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p00vMixStudyHLT.root", + "pikp/FBE6r01vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.2.p02 FTFP_BERT_EMM", + "10.4.p03 FTFP_BERT_EMM", + "10.6.ref00 FTFP_BERT_EMM + Birk", + "10.6.ref01 FTFP_BERT_EMM + Birk"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p01AvMixStudyHLT.root", + "pikp/FBE6p01BvMixStudyHLT.root", + "pikp/FBE6p01CvMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4.p03 FTFP_BERT_EMM", + "10.6.p01 FTFP_BERT_EMM + Birk C1+C3", + "10.6.p01 FTFP_BERT_EMM + Birk C1", + "10.6.p01 FTFP_BERT_EMM + Birk C1 + #pi"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p01CvMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4.p03 FTFP_BERT_EMM", + "10.6.p01 FTFP_BERT_EMM + Birk C1 + #pi"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE5r08vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4.p03 FTFP_BERT_EMM", + "10.5.ref08 FTFP_BERT_EMM"}; +*/ +/* +std::string filem[nmodelm]={"pikp/QFBE2p2StudyHLT.root", + "pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB6p00vMixStudyHLT.root", + "pikp/QFB6r01vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.2.p02 QGSP_FTFP_BERT_EMM", + "10.4.p03 QGSP_FTFP_BERT_EML", + "10.6.ref00 QGSP_FTFP_BERT_EML + Birk", + "10.6.ref01 QGSP_FTFP_BERT_EML + Birk"}; +*/ +/* +std::string filem[nmodelm]={"pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB6p01AvMixStudyHLT.root", + "pikp/QFB6p01BvMixStudyHLT.root", + "pikp/QFB6p01CvMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4.p03 QGSP_FTFP_BERT_EML", + "10.6.p01 QGSP_FTFP_BERT_EML + Birk C1+C3", + "10.6.p01 QGSP_FTFP_BERT_EML + Birk C1", + "10.6.p01 QGSP_FTFP_BERT_EML + Birk C1 + #pi"}; +*/ +std::string filem[nmodelm] = {"pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE7p02vMixStudyHLT.root", + "pikp/FBE110p01vMixStudyHLT.root", + "pikp/FBE110r05MixStudyHLT.root", + "pikp/FBE110r06vMixStudyHLT.root"}; +std::string typem[nmodelm] = {"10.4.p03 FTFP_BERT_EMM", + "10.7.p02 FTFP_BERT_EMM + Birk C1 + #pi", + "11.0.p01 FTFP_BERT_EMM + Birk C1 + #pi", + "11.0.ref05 FTFP_BERT_EMM + Birk C1 + #pi", + "11.0.ref06 FTFP_BERT_EMM + Birk C1 + #pi"}; +/* +std::string filem[nmodelm]={"pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB7p02vMixStudyHLT.root", + "pikp/QFB110p01vMixStudyHLT.root", + "pikp/QFB110r05MixStudyHLT.root", + "pikp/QFB110r06vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4.p03 QGSP_FTFP_BERT_EML", + "10.7.p02 QGSP_FTFP_BERT_EML + Birk C1", + "11.0.p01 QGSP_FTFP_BERT_EML + Birk C1", + "11.0.ref05 QGSP_FTFP_BERT_EML + Birk C1", + "11.0.ref06 QGSP_FTFP_BERT_EML + Birk C1"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE7r00vMixStudyHLT.root", + "pikp/FBE7p01MixStudyHLT.root", + "pikp/FBE7p02vMixStudyHLT.root", + "pikp/FBE7p03vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.7 FTFP_BERT_EMM + Birk C1 + #pi", + "10.7.p01 FTFP_BERT_EMM + Birk C1 + #pi", + "10.7.p02 FTFP_BERT_EMM + Birk C1 + #pi", + "10.7.p03 FTFP_BERT_EMM + Birk C1 + #pi"}; +*/ +/* +std::string filem[nmodelm]={"pikp/QFB7r00vMixStudyHLT.root", + "pikp/QFB7p01MixStudyHLT.root", + "pikp/QFB7p02vMixStudyHLT.root", + "pikp/QFB7p03vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.7 QGSP_BERT_EML + Birk C1 + #pi", + "10.7.p01 QGSP_BERT_EML + Birk C1 + #pi", + "10.7.p02 QGSP_BERT_EML + Birk C1 + #pi", + "10.7.p03 QGSP_BERT_EML + Birk C1 + #pi"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p2vMixStudyHLT.root", + "pikp/FBE7p1vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4p03 FTFP_BERT_EMM", + "10.6p02 FTFP_BERT_EMM", + "10.7p01 FTFP_BERT_EMM"}; +*/ +/* +std::string filem[nmodelm]={"pikp/FBE4p3vMixStudyHLT.root"}; +std::string typem[nmodelm]={"10.4p03 FTFP_BERT_EMM"}; +*/ /* std::string filex[nmodelx]={"AllDataStudyHLT.root", "pikp/FBE2p2StudyHLT.root", @@ -435,6 +599,7 @@ std::string typex[nmodelx]={"Data (2016B)", "10.3.ref06 FTFP_BERT_EMM (VecGeom 4)", "10.3.ref06 FTFP_BERT_EMM (VecGeom Corr)"}; */ +/* std::string filex[nmodelx]={"AllDataStudyHLT.root", "pikp/FBE4r00vMixStudyHLT.root", "pikp/FBE4r01MixStudyHLT.root", @@ -443,6 +608,213 @@ std::string typex[nmodelx]={"Data (2016B)", "10.4 VecGeom FTFP_BERT_EMM", "10.4.ref01 FTFP_BERT_EMM", "10.4.ref01 VecGeom FTFP_BERT_EMM"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4MixStudyHLT10d.root", + "pikp/FBE4MixStudyHLT10s.root", + "pikp/FBE4MixStudyHLT10t.root", + "pikp/FBE4MixStudyHLT10st.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4 FTFP_BERT_EMM (Default)", + "10.4 FTFP_BERT_EMM (New StepIntegrator)", + "10.4 FTFP_BERT_EMM (Smart Track)", + "10.4 FTFP_BERT_EMM (New Stepper + Smart Track)"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE2p2StudyHLT.root", + "pikp/FBE4MixStudyHLT1025.root", + "pikp/FBE5c01vMixStudyHLT.root", + "pikp/FBE5c02vMixStudyHLT.root", + "pikp/QFB5c02vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.2.p02 FTFP_BERT_EMM", + "10.4 FTFP_BERT_EMM", + "10.5.cand01 FTFP_BERT_EMM", + "10.5.cand02 FTFP_BERT_EMM", + "10.5.cand02 QGSP_FTFP_BERT_EML"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE5r00vMixStudyHLT.root", + "pikp/FBE5r05vMixStudyHLT.root", + "pikp/FBE5r07vMixStudyHLT.root", + "pikp/FBE5r08vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 FTFP_BERT_EMM", + "10.5 FTFP_BERT_EMM", + "10.5.ref05 FTFP_BERT_EMM", + "10.5.ref07 FTFP_BERT_EMM", + "10.5.ref08 FTFP_BERT_EMM"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB5r00vMixStudyHLT.root", + "pikp/QFB5r05vMixStudyHLT.root", + "pikp/QFB5r07vMixStudyHLT.root", + "pikp/QFB5r08vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 QGSP_FTFP_BERT_EML", + "10.5 QGSP_FTFP_BERT_EML", + "10.5.ref05 QGSP_FTFP_BERT_EML", + "10.5.ref07 QGSP_FTFP_BERT_EML", + "10.5.ref08 QGSP_FTFP_BERT_EML"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE5r10vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 FTFP_BERT_EMM", + "10.5.ref10 FTFP_BERT_EMM"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE2p2StudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p00vMixStudyHLT.root", + "pikp/FBE6r01vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.2.p02 FTFP_BERT_EMM", + "10.4.p03 FTFP_BERT_EMM", + "10.6.ref00 FTFP_BERT_EMM + Birk", + "10.6.ref01 FTFP_BERT_EMM + Birk"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p01AvMixStudyHLT.root", + "pikp/FBE6p01BvMixStudyHLT.root", + "pikp/FBE6p01CvMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 FTFP_BERT_EMM", + "10.6.p01 FTFP_BERT_EMM + Birk C1+C3", + "10.6.p01 FTFP_BERT_EMM + Birk C1", + "10.6.p01 FTFP_BERT_EMM + Birk C1 + #pi"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p01CvMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 FTFP_BERT_EMM", + "10.6.p01 FTFP_BERT_EMM + Birk C1 + #pi"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB5r00vMixStudyHLT.root", + "pikp/QFB6p00vMixStudyHLT.root", + "pikp/QFB6r01vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 QGSP_FTFP_BERT_EML", + "10.5 QGSP_FTFP_BERT_EML", + "10.6.ref00 QGSP_FTFP_BERT_EML + Birk", + "10.6.ref01 QGSP_FTFP_BERT_EML + Birk"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB6p01AvMixStudyHLT.root", + "pikp/QFB6p01BvMixStudyHLT.root", + "pikp/QFB6p01CvMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 QGSP_FTFP_BERT_EML", + "10.6.p01 QGSP_FTFP_BERT_EML + Birk C1+C3", + "10.6.p01 QGSP_FTFP_BERT_EML + Birk C1", + "10.6.p01 QGSP_FTFP_BERT_EML + Birk C1 + #pi"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE7p02vMixStudyHLT.root", + "pikp/FBE110p01vMixStudyHLT.root", + "pikp/FBE110r05MixStudyHLT.root", + "pikp/FBE110r06vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 FTFP_BERT_EMM", + "10.7.p02 FTFP_BERT_EMM", + "11.0.p01 FTFP_BERT_EMM", + "11.0.ref05 FTFP_BERT_EMM", + "11.0.ref06 FTFP_BERT_EMM"}; +*/ +std::string filex[nmodelx] = {"AllDataStudyHLT.root", + "pikp/QFB4p3vMixStudyHLT.root", + "pikp/QFB7p02vMixStudyHLT.root", + "pikp/QFB110p01vMixStudyHLT.root", + "pikp/QFB110r05MixStudyHLT.root", + "pikp/QFB110r06vMixStudyHLT.root"}; +std::string typex[nmodelx] = {"Data (2016B)", + "10.4.p03 QGSP_FTFP_BERT_EML", + "10.7.p02 QGSP_FTFP_BERT_EML", + "11.0.p01 QGSP_FTFP_BERT_EML", + "11.0.ref05 QGSP_FTFP_BERT_EML", + "11.0.ref06 QGSP_FTFP_BERT_EML"}; +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/QFB110p01vMixStudyHLT.root", + "pikp/QFB110r01vMixStudyHLT.root", + "pikp/QFB110r02vMixStudyHLT.root", + "pikp/QFB110r03MixStudyHLT.root", + "pikp/QFB110r04MixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "11.0.p01 QGSP_BERT_EMM + Birk C1 + #pi", + "11.0.ref01 QGSP_BERT_EMM + Birk C1 + #pi", + "11.0.ref02 QGSP_BERT_EMM + Birk C1 + #pi", + "11.0.ref03 QGSP_BERT_EMM + Birk C1 + #pi", + "11.0.ref04 QGSP_BERT_EMM + Birk C1 + #pi"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/QFB7r00vMixStudyHLT.root", + "pikp/QFB7p01vMixStudyHLT.root", + "pikp/QFB7p02vMixStudyHLT.root", + "pikp/QFB7p03vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.7 QGSP_BERT_EML + Birk C1 + #pi", + "10.7.p01 QGSP_BERT_EML + Birk C1 + #pi", + "10.7.p02 QGSP_BERT_EML + Birk C1 + #pi", + "10.7.p03 QGSP_BERT_EML + Birk C1 + #pi"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE1250p30vMixStudyHLT.root", + "pikp/FBE1250p31vMixStudyHLT.root", + "pikp/FBE1250p32vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.7.p02 FTFP_BERT_EMM Default", + "10.7.p02 FTFP_BERT_EMM Emax = 0.02", + "10.7.p02 FTFP_BERT_EMM Emax = 0.04"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE6p2MixStudyHLT.root", + "pikp/FBE7p1MixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "Geant4.10.4p03", + "Geant4.10.6p02", + "Geant4.10.7p01"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/FBE4p3vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "Geant4.10.4p03"}; +*/ +/* +std::string filex[nmodelx]={"AllDataStudyHLT.root", + "pikp/QFB4vp3vMixStudyHLT.root", + "pikp/QFB6vp2vMixStudyHLT.root", + "pikp/QFB7vb1vMixStudyHLT.root"}; +std::string typex[nmodelx]={"Data (2016B)", + "10.4.p03 QGSP_FTFP_BERT_EML", + "10.6.p02 QGSP_FTFP_BERT_EML", + "10.7.beta QGSP_FTFP_BERT_EML"}; +*/ /* std::string files[nmodelx]={"StudyHLT_ZeroBias_1PV.root","StudyHLT_PixelTrack_1PV.root","StudyHLT_1PV.root"}; std::string types[nmodels]={"Zero Bias HLT","Pixel Track HLT","All HLTs"}; @@ -454,458 +826,842 @@ std::string typex[nmodelx]={"Data (2016B)", std::string typex[nmodelx]={"Data (2016B)","10.0.p02 QGSP_FTFP_BERT_EML","10.2.p02 FTFP_BERT_EMM"}; */ -void plotAll(std::string fname="hlt.root", std::string HLT="All HLTs", int var=-1, int ien=-1, int eta=-1, bool varbin=false, int rebin=1, bool approve=true, bool logy=true, int pos=0, bool pv=false, int savePlot=-1); -void plotEnergyAll(std::string fname="", std::string hlt="All HLTs", int models=15, int pv=0, int data=4, bool varbin=false, int rebin=5, bool approve=true, bool logy=true, int pos=0, int var=-1, int ene=-1, int eta=-1, int savePlot=-1); -void plotEMeanAll(int data=4, int models=15, bool ratio=false, bool approve=true, std::string postfix="", int savePlot=-1); -void plotEMean(std::string fname="", std::string hlt="All HLTs", int models=15, int var=0, int eta=0, int pv=0, int dataMC=1, bool raio=false, bool approve=true, std::string postfix="",int savePlot=-1); -TCanvas* plotEMeanDraw(std::vector fnames, std::vector hlts, int var, int eta, int pv=0, bool approve=false, std::string dtype="Data", int coloff=0); -TCanvas* plotEMeanRatioDraw(std::vector fnames, std::vector hlts, int var, int eta, int pv=0, bool approve=false, std::string dtype="Data", int coloff=0); -TCanvas* plotEnergies(std::vector fnames, std::vector hlts, int var=0, int ien=0, int eta=0, int pv=0, bool varbin=false, int rebin=1, bool approve=false, std::string dtype="Data", bool logy=true, int pos=0, int coloff=0); -TCanvas* plotEnergy(std::string fname="hlt.root", std::string HLT="All HLTs", int var=0, int ien=0, int eta=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0, int coloff=0); -void plotEMeanPVAll(std::string fname="StudyHLT_ZeroBias_1PV.root", std::string HLT="Zero Bias", int var=-1, int eta=-1, bool approve=true); -TCanvas* plotEMeanDrawPV(std::string fname="StudyHLT_ZeroBias_1PV.root", std::string HLT="Zero Bias", int var=0, int eta=0, bool approve=true); -TCanvas* plotEnergyPV(std::string fnamee="StudyHLT_HLTZeroBias.root", std::string HLT="Zero Bias", int var=0, int ien=0, int eta=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0); -TCanvas* plotTrack(std::string fname="hlt.root", std::string HLT="All HLTs", int var=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0); -TCanvas* plotIsolation(std::string fname="hlt.root", std::string HLT="All HLTs", int var=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0); -TCanvas* plotHLT(std::string fname="hlt.root", std::string HLT="All HLTs", int run=-1, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0); -TCanvas* plotHisto(char* cname, std::string HLT, TObjArray& histArr, std::vector& labels, std::vector& color, char* name, double ymx0, bool logy, int pos, double yloff, double yhoff, double xmax=-1, bool varbin=false, int rebin=1, bool approve=false); -void getHistStats(TH1D *h, int& entries, int& integral, double& mean, double& meanE, double& rms, double& rmsE, int& uflow, int& oflow); -TFitResultPtr getHistFitStats(TH1F *h, const char* formula, double xlow, double xup, unsigned int& nPar, double* par, double* epar); -void setHistAttr(TH1F *h, int icol, int lwid=1, int ltype=1); -double getWeightedMean (int npt, int Start, std::vector& mean, std::vector& emean); +void plotAll(std::string fname = "", + std::string HLT = "", + int var = -1, + int ien = -1, + int eta = -1, + bool varbin = false, + int rebin = 1, + bool approve = true, + bool logy = true, + int pos = 0, + bool pv = false, + int savePlot = -1); +void plotEnergyAll(std::string fname = "", + std::string hlt = "All HLTs", + int models = 15, + int pv = 0, + int data = 4, + bool varbin = false, + int rebin = 5, + bool approve = true, + bool logy = true, + int pos = 0, + int var = -1, + int ene = -1, + int eta = -1, + int savePlot = -1); +void plotEMeanAll( + int data = 4, int models = 63, bool ratio = true, bool approve = true, std::string postfix = "F", int savePlot = 2); +void plotEMean(std::string fname = "", + std::string hlt = "All HLTs", + int models = 15, + int var = 0, + int eta = 0, + int pv = 0, + int dataMC = 1, + bool raio = false, + bool approve = true, + std::string postfix = "", + int savePlot = -1); +TCanvas* plotEMeanDraw(std::vector fnames, + std::vector hlts, + int var, + int eta, + int pv = 0, + bool approve = false, + std::string dtype = "Data", + int coloff = 0); +TCanvas* plotEMeanRatioDraw(std::vector fnames, + std::vector hlts, + int var, + int eta, + int pv = 0, + bool approve = false, + std::string dtype = "Data", + int coloff = 0); +void plotEMeanRatioXAll(std::string fname = "pikp/FBE4p3vMixStudyHLT.root", + std::string hlt = "10.4p03 FTFP_BERT_EMM", + std::string postfix = "F", + int savePlot = 2); +TCanvas* plotEMeanRatioDrawX(std::string fname = "pikp/FBE4p3vMixStudyHLT.root", + std::string hlt = "10.4p03 FTFP_BERT_EMM", + int var = 0, + int pv = 0); +TCanvas* plotEnergies(std::vector fnames, + std::vector hlts, + int var = 0, + int ien = 0, + int eta = 0, + int pv = 0, + bool varbin = false, + int rebin = 1, + bool approve = false, + std::string dtype = "Data", + bool logy = true, + int pos = 0, + int coloff = 0); +TCanvas* plotEnergy(std::string fname = "hlt.root", + std::string HLT = "All HLTs", + int var = 0, + int ien = 0, + int eta = 0, + bool varbin = false, + int rebin = 1, + bool approve = false, + bool logy = true, + int pos = 0, + int coloff = 0); +void plotEMeanPVAll(std::string fname = "StudyHLT_ZeroBias_1PV.root", + std::string HLT = "Zero Bias", + int var = -1, + int eta = -1, + bool approve = true); +TCanvas* plotEMeanDrawPV(std::string fname = "StudyHLT_ZeroBias_1PV.root", + std::string HLT = "Zero Bias", + int var = 0, + int eta = 0, + bool approve = true); +TCanvas* plotEnergyPV(std::string fnamee = "StudyHLT_HLTZeroBias.root", + std::string HLT = "Zero Bias", + int var = 0, + int ien = 0, + int eta = 0, + bool varbin = false, + int rebin = 1, + bool approve = false, + bool logy = true, + int pos = 0); +TCanvas* plotTrack(std::string fname = "hlt.root", + std::string HLT = "All HLTs", + int var = 0, + bool varbin = false, + int rebin = 1, + bool approve = false, + bool logy = true, + int pos = 0); +TCanvas* plotIsolation(std::string fname = "hlt.root", + std::string HLT = "All HLTs", + int var = 0, + bool varbin = false, + int rebin = 1, + bool approve = false, + bool logy = true, + int pos = 0); +TCanvas* plotHLT(std::string fname = "hlt.root", + std::string HLT = "All HLTs", + int run = -1, + bool varbin = false, + int rebin = 1, + bool approve = false, + bool logy = true, + int pos = 0); +TCanvas* plotHisto(char* cname, + std::string HLT, + TObjArray& histArr, + std::vector& labels, + std::vector& color, + char* name, + double ymx0, + bool logy, + int pos, + double yloff, + double yhoff, + double xmax = -1, + bool varbin = false, + int rebin = 1, + bool approve = false); +void getHistStats(TH1D* h, + int& entries, + int& integral, + double& mean, + double& meanE, + double& rms, + double& rmsE, + int& uflow, + int& oflow); +TFitResultPtr getHistFitStats( + TH1F* h, const char* formula, double xlow, double xup, unsigned int& nPar, double* par, double* epar); +void setHistAttr(TH1F* h, int icol, int lwid = 1, int ltype = 1); +double getWeightedMean(int npt, int Start, std::vector& mean, std::vector& emean); TH1D* rebin(TH1D* histin, int); -void plotAll(std::string fname, std::string HLT, int var, int ien, int eta, - bool varbin, int rebin, bool approve, bool logy, int pos, bool pv, - int savePlot) { - int varmin(0), varmax(5), enemin(0), enemax(9), etamin(0), etamax(3); - if (var >= 0) varmin = varmax = var; - if (ien >= 0) enemin = enemax = ien; - if (eta >= 0) etamin = etamax = eta; -// std::cout << "Var " << varmin << ":" << varmax << " Ene " << enemin << ":" -// << enemax << " Eta " << etamin << ":" << etamax << std::endl; - for (int var=varmin; var<=varmax; ++var) { - for (int ene=enemin; ene<=enemax; ++ene) { - for (int eta=etamin; eta<=etamax; ++eta) { - TCanvas *c(0); - if (pv) { - c = plotEnergyPV(fname, HLT, var, ene, eta, varbin, rebin, approve, logy, pos); - } else { - c = plotEnergy(fname, HLT, var, ene, eta, varbin, rebin, approve, logy, pos); - } - if (c != 0 && savePlot >= 0 && savePlot < 3) { - std::string ext[3] = {"eps", "gif", "pdf"}; - char name[200]; - sprintf (name, "%s.%s", c->GetName(), ext[savePlot].c_str()); - c->Print(name); - } +void plotAll(std::string fname, + std::string HLT, + int var, + int ien, + int eta, + bool varbin, + int rebin, + bool approve, + bool logy, + int pos, + bool pv, + int savePlot) { + int varmin(0), varmax(5), enemin(0), enemax(9), etamin(0), etamax(etaMax); + if (var >= 0) + varmin = varmax = var; + if (ien >= 0) + enemin = enemax = ien; + if (eta >= 0) + etamin = etamax = eta; + for (int var = varmin; var <= varmax; ++var) { + for (int ene = enemin; ene <= enemax; ++ene) { + for (int eta = etamin; eta <= etamax; ++eta) { + TCanvas* c(0); + if (pv) { + c = plotEnergyPV(fname, HLT, var, ene, eta, varbin, rebin, approve, logy, pos); + } else { + c = plotEnergy(fname, HLT, var, ene, eta, varbin, rebin, approve, logy, pos); + } + if (c != 0 && savePlot >= 0 && savePlot <= 3) { + std::string ext[4] = {"eps", "gif", "pdf", "C"}; + char name[200]; + sprintf(name, "%s.%s", c->GetName(), ext[savePlot].c_str()); + c->Print(name); + } } } } } -void plotEnergyAll(std::string fname, std::string hlt, int models, int pv, - int data, bool varbin, int rebin, bool approve, bool logy, - int pos, int var, int ene, int eta, int savePlot) { - +void plotEnergyAll(std::string fname, + std::string hlt, + int models, + int pv, + int data, + bool varbin, + int rebin, + bool approve, + bool logy, + int pos, + int var, + int ene, + int eta, + int savePlot) { std::vector fnames, hlts; std::string dtype = (data == 1) ? "Data" : "MC"; int modeluse(models); int coloff = (data == 4) ? 0 : 1; if (fname == "") { if (data == 1) { - for (int i=0; i= varmin && var <= varmax) varmin = varmax = var; - if (ene >= enemin && ene <= enemax) enemin = enemax = ene; - if (eta >= etamin && eta <= etamax) etamin = etamax = eta; -// std::cout << "Var " << varmin << ":" << varmax << " Ene " << enemin << ":" -// << enemax << " Eta " << etamin << ":" << etamax << std::endl; - for (int var=varmin; var<=varmax; ++var) { - for (int ene=enemin; ene<=enemax; ++ene) { - for (int eta=etamin; eta<=etamax; ++eta) { - TCanvas *c = plotEnergies(fnames, hlts, var, ene, eta, pv, varbin, - rebin, approve, dtype, logy, pos, coloff); - if (c != 0 && savePlot >= 0 && savePlot < 3) { - std::string ext[3] = {"eps", "gif", "pdf"}; - char name[200]; - sprintf (name, "%s.%s", c->GetName(), ext[savePlot].c_str()); - c->Print(name); - } + int varmin(0), varmax(5), enemin(0), enemax(9), etamin(0), etamax(etaMax); + if (var >= varmin && var <= varmax) + varmin = varmax = var; + if (ene >= enemin && ene <= enemax) + enemin = enemax = ene; + if (eta >= etamin && eta <= etamax) + etamin = etamax = eta; + for (int var = varmin; var <= varmax; ++var) { + for (int ene = enemin; ene <= enemax; ++ene) { + for (int eta = etamin; eta <= etamax; ++eta) { + TCanvas* c = plotEnergies(fnames, hlts, var, ene, eta, pv, varbin, rebin, approve, dtype, logy, pos, coloff); + if (c != 0 && savePlot >= 0 && savePlot <= 3) { + std::string ext[4] = {"eps", "gif", "pdf", "C"}; + char name[200]; + sprintf(name, "%s.%s", c->GetName(), ext[savePlot].c_str()); + c->Print(name); + } } } } } -void plotEMeanAll(int data, int models, bool ratio, bool approve, - std::string postfix, int savePlot){ - int varmin(0), varmax(5), pvmin(0), pvmax(0), etamin(0), etamax(3); -// std::cout << "Var " << varmin << ":" << varmax << " PV " << pvmin << ":" -// << pvmax << " Eta " << etamin << ":" << etamax << std::endl; - for (int var=varmin; var<=varmax; ++var) { - for (int eta=etamin; eta<=etamax; ++eta) { - for (int pv=pvmin; pv<=pvmax; ++pv) { - plotEMean("", "", models, var, eta, pv, data, ratio, approve, postfix, - savePlot); +void plotEMeanAll(int data, int models, bool ratio, bool approve, std::string postfix, int savePlot) { + int varmin(0), varmax(5), pvmin(0), pvmax(0), etamin(0), etamax(etaMax); + for (int var = varmin; var <= varmax; ++var) { + for (int eta = etamin; eta <= etamax; ++eta) { + for (int pv = pvmin; pv <= pvmax; ++pv) { + plotEMean("", "", models, var, eta, pv, data, ratio, approve, postfix, savePlot); } } } } -void plotEMean(std::string fname, std::string hlt, int models, int var, int eta, - int pv, int data, bool ratio, bool approve, std::string postfix, - int savePlot) { - +void plotEMean(std::string fname, + std::string hlt, + int models, + int var, + int eta, + int pv, + int data, + bool ratio, + bool approve, + std::string postfix, + int savePlot) { std::vector fnames, hlts; std::string dtype = (data == 1) ? "Data" : "MC"; int modeluse(models); int coloff = (data == 4 || data == 3) ? 0 : 1; if (fname == "") { if (data == 1) { - for (int i=0; i= 0) varmin = varmax = var; - if (eta >= 0) etamin = etamax = eta; - if (pv >= 0) pvmin = pvmax = pv; -// std::cout << "Var " << varmin << ":" << varmax << " Eta " << etamin << ":" -// << etamax << std::endl; - for (int var=varmin; var<=varmax; ++var) { - for (int eta=etamin; eta<=etamax; ++eta) { - for (int pv=pvmin; pv<=pvmax; ++pv) { - TCanvas* c(0); - if (ratio) { - c = plotEMeanRatioDraw(fnames, hlts, var, eta, pv, approve, dtype, coloff); - } else { - c = plotEMeanDraw(fnames, hlts, var, eta, pv, approve, dtype, coloff); - } - if (c != 0 && savePlot >= 0 && savePlot < 3) { - std::string ext[3] = {"eps", "gif", "pdf"}; - char name[200]; - sprintf (name, "%s%s.%s", c->GetName(), postfix.c_str(), - ext[savePlot].c_str()); - c->Print(name); - } + int varmin(0), varmax(5), etamin(0), etamax(etaMax), pvmin(0), pvmax(4); + if (var >= 0) + varmin = varmax = var; + if (eta >= 0) + etamin = etamax = eta; + if (pv >= 0) + pvmin = pvmax = pv; + for (int var = varmin; var <= varmax; ++var) { + for (int eta = etamin; eta <= etamax; ++eta) { + for (int pv = pvmin; pv <= pvmax; ++pv) { + TCanvas* c = ((ratio) ? plotEMeanRatioDraw(fnames, hlts, var, eta, pv, approve, dtype, coloff) + : plotEMeanDraw(fnames, hlts, var, eta, pv, approve, dtype, coloff)); + if (c != 0 && savePlot >= 0 && savePlot <= 3) { + std::string ext[4] = {"eps", "gif", "pdf", "C"}; + char name[200]; + sprintf(name, "%s%s.%s", c->GetName(), postfix.c_str(), ext[savePlot].c_str()); + c->Print(name); + } } } } } -TCanvas* plotEMeanDraw(std::vector fnames, - std::vector hlts, int var, int eta, int pv, - bool approve, std::string dtype, int coloff) { - +TCanvas* plotEMeanDraw(std::vector fnames, + std::vector hlts, + int var, + int eta, + int pv, + bool approve, + std::string dtype, + int coloff) { bool debug(false); std::vector graphs; - TLegend* legend = new TLegend(0.25, 0.80, 0.975, 0.95); - legend->SetBorderSize(1); legend->SetFillColor(kWhite); + double yminx = (fnames.size() < 3) ? 0.85 : 0.75; + TLegend* legend = new TLegend(0.60, yminx, 0.975, 0.95); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); legend->SetMargin(0.2); - for (unsigned k=0; kFindObjectAny(name); + for (int i = 0; i < NPT; ++i) { + char name[100]; + sprintf(name, "h_energy_%d_%d_%d_%d", pv + 3, i, eta, var); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { - mean[i] = histo->GetMean(); - dmean[i]= histo->GetMeanError(); + mean[i] = histo->GetMean(); + dmean[i] = histo->GetMeanError(); } else { - mean[i] = -100.; - dmean[i]= 0; + mean[i] = -100.; + dmean[i] = 0; } } if (debug) { std::cout << "Get mean for " << NPT << " points" << std::endl; - for (int i=0; iSetMarkerStyle(styles[coloff+k]); - graph->SetMarkerColor(colors[coloff+k]); - graph->SetMarkerSize(1.6); - graph->SetLineColor(colors[coloff+k]); + TGraphAsymmErrors* graph = new TGraphAsymmErrors(NPT, mom, mean, dmom, dmom, dmean, dmean); + graph->SetMarkerStyle(styles[coloff + k]); + graph->SetMarkerColor(colors[coloff + k]); + graph->SetMarkerSize(1.2); + graph->SetLineColor(colors[coloff + k]); graph->SetLineWidth(2); graphs.push_back(graph); legend->AddEntry(graph, hlts[k].c_str(), "lp"); - if (debug) std::cout << "Complete " << hlts[k] << std::endl; + if (debug) + std::cout << "Complete " << hlts[k] << std::endl; file->Close(); } char cname[100], name[200]; - sprintf (cname, "c_%s_%d_%d_%s", varEne1[var].c_str(), eta, pv, dtype.c_str()); - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptTitle(kFALSE); gStyle->SetPadBorderMode(0); + sprintf(cname, "c_%s_%d_%d_%s", varEne1[var].c_str(), eta, pv, dtype.c_str()); gStyle->SetCanvasBorderMode(0); - TCanvas *canvas = new TCanvas(cname, cname, 500, 400); - gStyle->SetOptStat(0); gPad->SetTopMargin(0.05); - gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(kFALSE); + gStyle->SetPadBorderMode(0); + gStyle->SetCanvasBorderMode(0); + TCanvas* canvas = new TCanvas(cname, cname, 500, 400); + gStyle->SetOptStat(0); + gPad->SetTopMargin(0.05); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); gPad->SetBottomMargin(0.20); - TH1F *vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5); - vFrame->GetYaxis()->SetRangeUser(0.0,1.6); + TH1F* vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5); + vFrame->GetYaxis()->SetRangeUser(0.0, 1.6); vFrame->GetXaxis()->SetLabelSize(0.06); vFrame->GetYaxis()->SetLabelSize(0.05); vFrame->GetXaxis()->SetTitleSize(0.06); vFrame->GetYaxis()->SetTitleSize(0.06); vFrame->GetYaxis()->SetTitleOffset(0.9); - vFrame->GetXaxis()->SetRangeUser(1.0,20.0); + vFrame->GetXaxis()->SetRangeUser(1.0, 20.0); if (approve) { - sprintf (name, "Mean of %s/p_{Track}", varEne[var].c_str()); + sprintf(name, "Mean of %s/p_{Track}", varEne[var].c_str()); } else { - sprintf (name, "<%s/p_{Track}>", varEne[var].c_str()); + sprintf(name, "<%s/p_{Track}>", varEne[var].c_str()); } - vFrame->GetYaxis()->SetTitle(name); - sprintf (name, "p_{Track} (GeV/c)"); + vFrame->GetYaxis()->SetTitle(name); + sprintf(name, "p_{Track} (GeV/c)"); vFrame->GetXaxis()->SetTitle(name); - for (unsigned int ii=0; iiDraw("P"); + for (unsigned int ii = 0; ii < graphs.size(); ++ii) + graphs[ii]->Draw("P"); legend->Draw(); - TLine *line = new TLine(1.0, 1.0, 20.0, 1.0); - line->SetLineStyle(2); line->SetLineWidth(2); line->SetLineColor(kRed); + TLine* line = new TLine(1.0, 1.0, 20.0, 1.0); + line->SetLineStyle(2); + line->SetLineWidth(2); + line->SetLineColor(kRed); line->Draw(); TPaveText* text = new TPaveText(0.25, 0.74, 0.55, 0.79, "brNDC"); if (approve) { - sprintf (name, "(%s)", nameEtas[eta].c_str()); + sprintf(name, "(%s)", nameEtas[eta].c_str()); } else { - sprintf (name, "(%s, %s)", nameEta[eta].c_str(), namePV[pv].c_str()); + sprintf(name, "(%s, %s)", nameEta[eta].c_str(), namePV[pv].c_str()); } - if (debug) std::cout << "Name " << name << " |" << std::endl; + if (debug) + std::cout << "Name " << name << " |" << std::endl; text->AddText(name); text->Draw("same"); - TPaveText* text2 = new TPaveText(0.55, 0.74, 0.97, 0.79, "brNDC"); - sprintf (name, "CMS Preliminary"); + TPaveText* text2 = new TPaveText(0.55, yminx - 0.06, 0.97, yminx - 0.01, "brNDC"); + sprintf(name, cmsP.c_str()); text2->AddText(name); text2->Draw("same"); return canvas; } -TCanvas* plotEMeanRatioDraw(std::vector fnames, - std::vector hlts, int var, int eta, - int pv, bool approve, std::string dtype, - int coloff) { - +TCanvas* plotEMeanRatioDraw(std::vector fnames, + std::vector hlts, + int var, + int eta, + int pv, + bool approve, + std::string dtype, + int coloff) { bool debug(false); std::vector graphs; - TLegend* legend = new TLegend(0.25, 0.80, 0.975, 0.95); - legend->SetBorderSize(1); legend->SetFillColor(kWhite); + double yminx = (fnames.size() < 3) ? 0.88 : 0.80; + TLegend* legend = new TLegend(0.60, yminx, 0.975, 0.95); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); legend->SetMargin(0.2); - double mean0[NPT], dmean0[NPT]; - for (unsigned k=0; kFindObjectAny(name); + for (int i = 0; i < NPT; ++i) { + char name[100]; + sprintf(name, "h_energy_%d_%d_%d_%d", pv + 3, i, eta, var); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { - mean[i] = histo->GetMean(); - dmean[i]= histo->GetMeanError(); + mean[i] = histo->GetMean(); + dmean[i] = histo->GetMeanError(); } else { - mean[i] = -100.; - dmean[i]= 0; + mean[i] = -100.; + dmean[i] = 0; } } if (debug) { std::cout << "Get mean for " << NPT << " points" << std::endl; - for (int i=0; i 0) ? dmean[i] / mean[i] : 0.0; + dmean1[NPT2 - i - 1] = -dmean1[i]; } + pmom1[NPT + 1] = pmom1[NPT] = mom[NPT - 1] + dmom[NPT - 1]; + dmean1[NPT] = dmean1[NPT - 1]; + dmean1[NPT + 1] = -dmean1[NPT]; + pmom1[NPT2] = pmom1[0]; + dmean1[NPT2] = dmean1[0]; + if (debug) + for (int k = 0; k < NPT2 + 1; ++k) + std::cout << "[" << k << "] " << pmom1[k] << " " << dmean1[k] << "\n"; } else { double sumNum(0), sumDen(0), sumNum1(0), sumDen1(0); - for (int i=0; i 0 && dmean0[i] > 0) { - double er1 = dmean[i]/mean[i]; - double er2 = dmean0[i]/mean0[i]; - mean[i] = mean[i]/mean0[i]; - dmean[i]= mean[i]*sqrt(er1*er1+er2*er2); - double temp1 = (mean[i]>1.0) ? 1.0/mean[i] : mean[i]; - double temp2 = (mean[i]>1.0) ? dmean[i]/(mean[i]*mean[i]) : dmean[i]; - if (i > 0) { - sumNum += (fabs(1-temp1)/(temp2*temp2)); - sumDen += (1.0/(temp2*temp2)); - } - sumNum1 += (fabs(1-temp1)/(temp2*temp2)); - sumDen1 += (1.0/(temp2*temp2)); - } else { - mean[i] = -100.; - dmean[i]= 0; - } + for (int i = 0; i < NPT; ++i) { + if (dmean[i] > 0 && dmean0[i] > 0) { + double er1 = dmean[i] / mean[i]; + double er2 = dmean0[i] / mean0[i]; + mean[i] = mean[i] / mean0[i]; + dmean[i] = mean[i] * sqrt(er1 * er1 + er2 * er2); + double temp1 = (mean[i] > 1.0) ? 1.0 / mean[i] : mean[i]; + double temp2 = (mean[i] > 1.0) ? dmean[i] / (mean[i] * mean[i]) : dmean[i]; + if (i > 0) { + sumNum += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen += (1.0 / (temp2 * temp2)); + } + sumNum1 += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen1 += (1.0 / (temp2 * temp2)); + } else { + mean[i] = -100.; + dmean[i] = 0; + } } - sumNum = (sumDen>0) ? (sumNum/sumDen) : 0; - sumDen = (sumDen>0) ? 1.0/sqrt(sumDen) : 0; - sumNum1 = (sumDen1>0) ? (sumNum1/sumDen1) : 0; - sumDen1 = (sumDen1>0) ? 1.0/sqrt(sumDen1) : 0; - std::cout << "Get Ratio of mean for " << NPT << " points: Mean " << sumNum << " +- " << sumDen << " (" << sumNum1 << " +- " << sumDen1 << ") Input: " << fnames[k] << " var|eta|pv " << var << ":" << eta << ":" << pv << std::endl; + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + sumNum1 = (sumDen1 > 0) ? (sumNum1 / sumDen1) : 0; + sumDen1 = (sumDen1 > 0) ? 1.0 / sqrt(sumDen1) : 0; + std::cout << "Get Ratio of mean for " << NPT << " points: Mean " << sumNum << " +- " << sumDen << " (" << sumNum1 + << " +- " << sumDen1 << ") Input: " << fnames[k] << " var|eta|pv " << var << ":" << eta << ":" << pv + << std::endl; if (debug) { - std::cout << "Get Ratio of mean for " << NPT << " points: Mean " - << sumNum << " +- " << sumDen << std::endl; - for (int i=0; iSetMarkerStyle(styles[coloff+k]); - graph->SetMarkerColor(colors[coloff+k]); - graph->SetMarkerSize(1.6); - graph->SetLineColor(colors[coloff+k]); + TGraphAsymmErrors* graph = new TGraphAsymmErrors(NPT, mom, mean, dmom, dmom, dmean, dmean); + graph->SetMarkerStyle(styles[coloff + k]); + graph->SetMarkerColor(colors[coloff + k]); + graph->SetMarkerSize(1.2); + graph->SetLineColor(colors[coloff + k]); graph->SetLineWidth(2); graphs.push_back(graph); char text[100]; if (approve) { - sprintf (text,"%s", hlts[k].c_str()); + sprintf(text, "%s", hlts[k].c_str()); } else { - sprintf (text,"%5.3f #pm %5.3f %s", sumNum, sumDen, hlts[k].c_str()); + sprintf(text, "%5.3f #pm %5.3f %s", sumNum, sumDen, hlts[k].c_str()); } legend->AddEntry(graph, text, "lp"); - if (debug) std::cout << "Complete " << hlts[k] << std::endl; + if (debug) + std::cout << "Complete " << hlts[k] << std::endl; } file->Close(); } + TGraphErrors* graph = new TGraphErrors(NPT, mom, ones, dmean1); + graph->SetFillColor(4); + graph->SetFillStyle(3005); + TPolyLine* pline = new TPolyLine(NPT2 + 1, pmom1, dmean1); + pline->SetFillColor(5); + pline->SetLineColor(6); + pline->SetLineWidth(2); char cname[100], name[200]; - sprintf (cname, "cR_%s_%d_%d_%s", varEne1[var].c_str(), eta, pv, dtype.c_str()); - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptTitle(kFALSE); gStyle->SetPadBorderMode(0); + sprintf(cname, "cR_%s_%d_%d_%s", varEne1[var].c_str(), eta, pv, dtype.c_str()); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(kFALSE); + gStyle->SetPadBorderMode(0); gStyle->SetCanvasBorderMode(0); - TCanvas *canvas = new TCanvas(cname, cname, 500, 400); - gStyle->SetOptStat(0); gPad->SetTopMargin(0.05); - gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025); + TCanvas* canvas = new TCanvas(cname, cname, 500, 400); + gStyle->SetOptStat(0); + gPad->SetTopMargin(0.05); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); gPad->SetBottomMargin(0.20); - TH1F *vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5); - vFrame->GetYaxis()->SetRangeUser(0.4,1.5); + TH1F* vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5); + vFrame->GetYaxis()->SetRangeUser(0.4, 1.5); vFrame->GetXaxis()->SetLabelSize(0.06); vFrame->GetYaxis()->SetLabelSize(0.05); vFrame->GetXaxis()->SetTitleSize(0.06); vFrame->GetYaxis()->SetTitleSize(0.045); vFrame->GetYaxis()->SetTitleOffset(1.2); - vFrame->GetXaxis()->SetRangeUser(1.0,20.0); + vFrame->GetXaxis()->SetRangeUser(1.0, 20.0); if (approve) { - sprintf (name, "%s/Data for mean of %s/p_{Track}", dtype.c_str(), varEne[var].c_str()); + sprintf(name, "%s/Data for mean of %s/p_{Track}", dtype.c_str(), varEne[var].c_str()); } else { - sprintf (name, "#frac{%s}{%s} for #frac{%s}{p_{Track}}", dtype.c_str(), hlts[0].c_str(), varEne[var].c_str()); + sprintf(name, "#frac{%s}{%s} for #frac{%s}{p_{Track}}", dtype.c_str(), hlts[0].c_str(), varEne[var].c_str()); } - vFrame->GetYaxis()->SetTitle(name); - sprintf (name, "p_{Track} (GeV/c)"); + vFrame->GetYaxis()->SetTitle(name); + sprintf(name, "p_{Track} (GeV/c)"); vFrame->GetXaxis()->SetTitle(name); - for (unsigned int ii=0; iiDraw("P"); + for (unsigned int ii = 0; ii < graphs.size(); ++ii) + graphs[ii]->Draw("P"); + graph->Draw("3"); legend->Draw(); - TLine *line = new TLine(1.0, 1.0, 20.0, 1.0); - line->SetLineStyle(2); line->SetLineWidth(2); line->SetLineColor(kRed); + TLine* line = new TLine(1.0, 1.0, 20.0, 1.0); + line->SetLineStyle(2); + line->SetLineWidth(2); + line->SetLineColor(kRed); line->Draw(); + pline->Draw("f L SAME"); TPaveText* text = new TPaveText(0.55, 0.40, 0.95, 0.45, "brNDC"); if (approve) { - sprintf (name, "(%s)", nameEtas[eta].c_str()); + sprintf(name, "(%s)", nameEtas[eta].c_str()); } else { - sprintf (name, "(%s, %s)", nameEta[eta].c_str(), namePV[pv].c_str()); + sprintf(name, "(%s, %s)", nameEta[eta].c_str(), namePV[pv].c_str()); } - if (debug) std::cout << "Name " << name << " |" << std::endl; + if (debug) + std::cout << "Name " << name << " |" << std::endl; text->AddText(name); text->Draw("same"); TPaveText* text2 = new TPaveText(0.55, 0.45, 0.95, 0.50, "brNDC"); - sprintf (name, "CMS Preliminary"); + sprintf(name, cmsP.c_str()); text2->AddText(name); text2->Draw("same"); return canvas; } -TCanvas* plotEnergies(std::vector fnames, - std::vector hlts, int var, int ien, int eta, - int pv, bool varbin, int rebin, bool approve, - std::string dtype, bool logy, int pos, int coloff) { +void plotEMeanRatioXAll(std::string fname, std::string hlt, std::string postfix, int savePlot) { + int varmin(0), varmax(5), pvmin(0), pvmax(0); + for (int var = varmin; var <= varmax; ++var) { + for (int pv = pvmin; pv <= pvmax; ++pv) { + TCanvas* c = plotEMeanRatioDrawX(fname, hlt, var, pv); + if (savePlot >= 0 && savePlot <= 3) { + std::string ext[4] = {"eps", "gif", "pdf", "C"}; + char name[200]; + sprintf(name, "%s%s.%s", c->GetName(), postfix.c_str(), ext[savePlot].c_str()); + c->Print(name); + } + } + } +} - // bool debug(false); - TLegend* legend = new TLegend(0.55, 0.70, 0.95, 0.85); - legend->SetBorderSize(1); legend->SetFillColor(kWhite); +TCanvas* plotEMeanRatioDrawX(std::string fname, std::string hlt, int var, int pv) { + bool debug(false); + std::vector graphs; + double yminx = 0.80; + TLegend* legend = new TLegend(0.60, yminx, 0.975, 0.95); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); + legend->SetMargin(0.2); + double mean0[NPT], dmean0[NPT], mean[NPT], dmean[NPT]; + char name[200], cname[100]; + for (int eta = 0; eta <= etaMax; ++eta) { + // First Data + TFile* file = TFile::Open(fileData.c_str()); + for (int i = 0; i < NPT; ++i) { + sprintf(name, "h_energy_%d_%d_%d_%d", pv + 3, i, eta, var); + TH1D* histo = (TH1D*)file->FindObjectAny(name); + if (histo) { + mean0[i] = histo->GetMean(); + dmean0[i] = histo->GetMeanError(); + } else { + mean0[i] = -100.; + dmean0[i] = 0; + } + } + file->Close(); + // Now MC + file = TFile::Open(fname.c_str()); + for (int i = 0; i < NPT; ++i) { + sprintf(name, "h_energy_%d_%d_%d_%d", pv + 3, i, eta, var); + TH1D* histo = (TH1D*)file->FindObjectAny(name); + if (histo) { + mean[i] = histo->GetMean(); + dmean[i] = histo->GetMeanError(); + } else { + mean[i] = -100.; + dmean[i] = 0; + } + } + file->Close(); + if (debug) { + std::cout << "Get mean for " << NPT << " points" << std::endl; + for (int i = 0; i < NPT; ++i) + std::cout << "[" << i << "]" + << " Momentum " << mom[i] << " +- " << dmom[i] << " Data:Mean " << mean0[i] << " +- " << dmean0[i] + << " MC:Mean " << mean[i] << " +- " << dmean[i] << std::endl; + } + double sumNum(0), sumDen(0), sumNum1(0), sumDen1(0); + for (int i = 0; i < NPT; ++i) { + if (dmean[i] > 0 && dmean0[i] > 0) { + double er1 = dmean[i] / mean[i]; + double er2 = dmean0[i] / mean0[i]; + mean[i] = mean[i] / mean0[i]; + dmean[i] = mean[i] * sqrt(er1 * er1 + er2 * er2); + double temp1 = (mean[i] > 1.0) ? 1.0 / mean[i] : mean[i]; + double temp2 = (mean[i] > 1.0) ? dmean[i] / (mean[i] * mean[i]) : dmean[i]; + if (i > 0) { + sumNum += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen += (1.0 / (temp2 * temp2)); + } + sumNum1 += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen1 += (1.0 / (temp2 * temp2)); + } else { + mean[i] = -100.; + dmean[i] = 0; + } + } + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + sumNum1 = (sumDen1 > 0) ? (sumNum1 / sumDen1) : 0; + sumDen1 = (sumDen1 > 0) ? 1.0 / sqrt(sumDen1) : 0; + std::cout << "Get Ratio of mean for " << NPT << " points: Mean " << sumNum << " +- " << sumDen << " (" << sumNum1 + << " +- " << sumDen1 << ") Input: " << fname << " var|eta|pv " << var << ":" << eta << ":" << pv + << std::endl; + if (debug) { + std::cout << "Get Ratio of mean for " << NPT << " points: Mean " << sumNum << " +- " << sumDen << std::endl; + for (int i = 0; i < NPT; ++i) + std::cout << "[" << i << "]" + << " Momentum " << mom[i] << " +- " << dmom[i] << " Mean " << mean[i] << " +- " << dmean[i] + << std::endl; + } + TGraphAsymmErrors* graph = new TGraphAsymmErrors(NPT, mom, mean, dmom, dmom, dmean, dmean); + graph->SetMarkerStyle(styles[eta]); + graph->SetMarkerColor(colors[eta]); + graph->SetMarkerSize(1.2); + graph->SetLineColor(colors[eta]); + graph->SetLineWidth(2); + graphs.push_back(graph); + sprintf(cname, "%s", nameEtas[eta].c_str()); + legend->AddEntry(graph, cname, "lp"); + } + + sprintf(cname, "cR_%s_%d", varEne1[var].c_str(), pv); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(kFALSE); + gStyle->SetPadBorderMode(0); + gStyle->SetCanvasBorderMode(0); + TCanvas* canvas = new TCanvas(cname, cname, 750, 600); + gStyle->SetOptStat(0); + gPad->SetTopMargin(0.05); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); + gPad->SetBottomMargin(0.20); + TH1F* vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5); + vFrame->GetYaxis()->SetRangeUser(0.4, 1.5); + vFrame->GetXaxis()->SetLabelSize(0.06); + vFrame->GetYaxis()->SetLabelSize(0.05); + vFrame->GetXaxis()->SetTitleSize(0.06); + vFrame->GetYaxis()->SetTitleSize(0.045); + vFrame->GetYaxis()->SetTitleOffset(1.2); + vFrame->GetXaxis()->SetRangeUser(1.0, 20.0); + sprintf(name, "MC/Data for mean of %s/p_{Track}", varEne[var].c_str()); + vFrame->GetYaxis()->SetTitle(name); + sprintf(name, "p_{Track} (GeV/c)"); + vFrame->GetXaxis()->SetTitle(name); + for (unsigned int ii = 0; ii < graphs.size(); ++ii) + graphs[ii]->Draw("P"); + legend->Draw(); + TLine* line = new TLine(1.0, 1.0, 20.0, 1.0); + line->SetLineStyle(2); + line->SetLineWidth(2); + line->SetLineColor(kRed); + line->Draw(); + TPaveText* text = new TPaveText(0.60, 0.40, 0.92, 0.45, "brNDC"); + sprintf(name, "%s", hlt.c_str()); + text->AddText(name); + text->Draw("same"); + TPaveText* text2 = new TPaveText(0.55, 0.45, 0.95, 0.50, "brNDC"); + sprintf(name, "%s", cmsP.c_str()); + text2->AddText(name); + text2->Draw("same"); + return canvas; +} + +TCanvas* plotEnergies(std::vector fnames, + std::vector hlts, + int var, + int ien, + int eta, + int pv, + bool varbin, + int rebin, + bool approve, + std::string dtype, + bool logy, + int pos, + int coloff) { + TLegend* legend = new TLegend(0.55, 0.70, 0.95, 0.85); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); legend->SetMargin(0.4); - TObjArray histArr; - char name[100]; + TObjArray histArr; + char name[100]; std::vector labels; - std::vector color; - double ymx0(0), ent0(0); - sprintf (name, "h_energy_%d_%d_%d_%d", pv+3, ien, eta, var); - for (unsigned k=0; kFindObjectAny(name); + std::vector color; + double ymx0(0), ent0(0); + sprintf(name, "h_energy_%d_%d_%d_%d", pv + 3, ien, eta, var); + for (unsigned int k = 0; k < fnames.size(); ++k) { + TFile* file = TFile::Open(fnames[k].c_str()); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { if (k == 0) { - ent0 = histo->GetEntries(); + ent0 = histo->GetEntries(); } else { - double scale = (histo->GetEntries() > 0) ? ent0/histo->GetEntries() : 1; - histo->Scale(scale); + double scale = (histo->GetEntries() > 0) ? ent0 / histo->GetEntries() : 1; + histo->Scale(scale); } histArr.AddLast(histo); labels.push_back(hlts[k]); - color.push_back(colors[coloff+k]); + color.push_back(colors[coloff + k]); int ibin = histo->GetMaximumBin(); - if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin); + if (histo->GetBinContent(ibin) > ymx0) + ymx0 = histo->GetBinContent(ibin); } } TCanvas* c(0); - if (histArr.GetEntries()>0) { - sprintf (name, "p=%s, %s", varPPs[ien].c_str(), nameEtas[eta].c_str()); + if (histArr.GetEntries() > 0) { + sprintf(name, "p=%s, %s", varPPs[ien].c_str(), nameEtas[eta].c_str()); /* if (approve) { sprintf (name, "p=%s, %s", varPPs[ien].c_str(), nameEtas[eta].c_str()); @@ -915,256 +1671,283 @@ TCanvas* plotEnergies(std::vector fnames, */ std::string clabel(name); char cname[50]; - sprintf (cname, "c_%s_%d_%d_%s", varEne1[var].c_str(), ien, eta, dtype.c_str()); - sprintf ( name, "%s/p", varEne[var].c_str()); - c = plotHisto(cname, clabel, histArr, labels, color, name, ymx0, logy, pos, - 0.10, 0.05, 2.5, varbin, rebin, approve); - } + sprintf(cname, "c_%s_%d_%d_%s", varEne1[var].c_str(), ien, eta, dtype.c_str()); + sprintf(name, "%s/p", varEne[var].c_str()); + c = plotHisto( + cname, clabel, histArr, labels, color, name, ymx0, logy, pos, 0.10, 0.05, 2.5, varbin, rebin, approve); + } return c; } -TCanvas* plotEnergy(std::string fname, std::string HLT, int var, int ien, - int eta, bool varbin, int rebin, bool approve, bool logy, - int pos, int coloff) { - - TFile *file = TFile::Open(fname.c_str()); - char name[100]; - TObjArray histArr; +TCanvas* plotEnergy(std::string fname, + std::string HLT, + int var, + int ien, + int eta, + bool varbin, + int rebin, + bool approve, + bool logy, + int pos, + int coloff) { + TFile* file = TFile::Open(fname.c_str()); + char name[100]; + TObjArray histArr; std::vector labels; - std::vector color; - double ymx0(0); - for (int i=0; i<4; ++i) { - sprintf (name, "h_energy_%d_%d_%d_%d", i, ien, eta, var); - TH1D *histo = (TH1D*) file->FindObjectAny(name); + std::vector color; + double ymx0(0); + for (int i = 0; i < 4; ++i) { + sprintf(name, "h_energy_%d_%d_%d_%d", i, ien, eta, var); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { histArr.AddLast(histo); - sprintf (name, "p=%s, #eta=%s %s", varPs[ien].c_str(), varEta[eta].c_str(), namefull[i+3].c_str()); + sprintf(name, "p=%s, #eta=%s %s", varPs[ien].c_str(), varEta[eta].c_str(), namefull[i + 3].c_str()); labels.push_back(name); - color.push_back(colors[coloff+i]); + color.push_back(colors[coloff + i]); int ibin = histo->GetMaximumBin(); - if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin); + if (histo->GetBinContent(ibin) > ymx0) + ymx0 = histo->GetBinContent(ibin); } } TCanvas* c(0); - if (histArr.GetEntries()>0) { + if (histArr.GetEntries() > 0) { char cname[50]; - sprintf (cname, "c_%s_%d_%d", varEne1[var].c_str(), ien, eta); - sprintf ( name, "%s/p", varEne[var].c_str()); + sprintf(cname, "c_%s_%d_%d", varEne1[var].c_str(), ien, eta); + sprintf(name, "%s/p", varEne[var].c_str()); c = plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, 0.10, 0.05, 2.5, varbin, rebin, approve); - } + } return c; } -void plotEMeanPVAll(std::string fname, std::string HLT, int var, int eta, - bool approve) { - - int varmin(0), varmax(5), etamin(0), etamax(3); - if (var >= 0) varmin = varmax = var; - if (eta >= 0) etamin = etamax = eta; - for (int var=varmin; var<=varmax; ++var) { - for (int eta=etamin; eta<=etamax; ++eta) { +void plotEMeanPVAll(std::string fname, std::string HLT, int var, int eta, bool approve) { + int varmin(0), varmax(5), etamin(0), etamax(etaMax); + if (var >= 0) + varmin = varmax = var; + if (eta >= 0) + etamin = etamax = eta; + for (int var = varmin; var <= varmax; ++var) { + for (int eta = etamin; eta <= etamax; ++eta) { plotEMeanDrawPV(fname, HLT, var, eta, approve); } } } -TCanvas* plotEMeanDrawPV(std::string fname, std::string HLT, int var, int eta, - bool approve) { - +TCanvas* plotEMeanDrawPV(std::string fname, std::string HLT, int var, int eta, bool approve) { bool debug(false); std::vector graphs; - TLegend* legend = new TLegend(0.575, 0.80, 0.975, 0.95); - legend->SetBorderSize(1); legend->SetFillColor(kWhite); + TLegend* legend = new TLegend(0.575, 0.80, 0.975, 0.95); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); legend->SetMargin(0.4); - TFile *file = TFile::Open(fname.c_str()); - const int nPVBin=4; - int pvBins[nPVBin+1] = {1, 2, 3, 5, 100}; - for (int k=0; kFindObjectAny(name); + for (int i = 0; i < NPT; ++i) { + sprintf(name, "h_energy_%d_%d_%d_%d", k, i, eta, var); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { - mean[i] = histo->GetMean(); - dmean[i]= histo->GetMeanError(); + mean[i] = histo->GetMean(); + dmean[i] = histo->GetMeanError(); } else { - mean[i] = -100.; - dmean[i]= 0; + mean[i] = -100.; + dmean[i] = 0; } } if (debug) { std::cout << "Get mean for " << NPT << " points" << std::endl; - for (int i=0; iSetMarkerStyle(styles[k]); graph->SetMarkerColor(colors[k]); - graph->SetMarkerSize(1.6); + graph->SetMarkerSize(1.2); graph->SetLineColor(colors[k]); graph->SetLineWidth(2); graphs.push_back(graph); - sprintf (name, "PV=%d:%d", pvBins[k], pvBins[k+1]-1); + sprintf(name, "PV=%d:%d", pvBins[k], pvBins[k + 1] - 1); legend->AddEntry(graph, name, "lp"); - if (debug) std::cout << "Complete " << name << std::endl; + if (debug) + std::cout << "Complete " << name << std::endl; } file->Close(); char cname[100], name[200]; - sprintf (cname, "c_%s_PV_%d", varEne1[var].c_str(), eta); - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptTitle(kFALSE); gStyle->SetPadBorderMode(0); + sprintf(cname, "c_%s_PV_%d", varEne1[var].c_str(), eta); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(kFALSE); + gStyle->SetPadBorderMode(0); gStyle->SetCanvasBorderMode(0); - TCanvas *canvas = new TCanvas(cname, cname, 500, 400); - gStyle->SetOptStat(0); gPad->SetTopMargin(0.05); - gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025); + TCanvas* canvas = new TCanvas(cname, cname, 500, 400); + gStyle->SetOptStat(0); + gPad->SetTopMargin(0.05); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); gPad->SetBottomMargin(0.20); - TH1F *vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5); - vFrame->GetYaxis()->SetRangeUser(0.0,1.5); + TH1F* vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5); + vFrame->GetYaxis()->SetRangeUser(0.0, 1.5); vFrame->GetXaxis()->SetLabelSize(0.06); vFrame->GetYaxis()->SetLabelSize(0.05); vFrame->GetXaxis()->SetTitleSize(0.06); vFrame->GetYaxis()->SetTitleSize(0.06); vFrame->GetYaxis()->SetTitleOffset(0.9); - vFrame->GetXaxis()->SetRangeUser(1.0,20.0); + vFrame->GetXaxis()->SetRangeUser(1.0, 20.0); if (approve) { - sprintf (name, "Mean of %s/p_{Track}", varEne[var].c_str()); + sprintf(name, "Mean of %s/p_{Track}", varEne[var].c_str()); } else { - sprintf (name, "<%s/p_{Track}>", varEne[var].c_str()); + sprintf(name, "<%s/p_{Track}>", varEne[var].c_str()); } - vFrame->GetYaxis()->SetTitle(name); - sprintf (name, "p_{Track} (GeV/c)"); + vFrame->GetYaxis()->SetTitle(name); + sprintf(name, "p_{Track} (GeV/c)"); vFrame->GetXaxis()->SetTitle(name); - for (unsigned int ii=0; iiDraw("P"); + for (unsigned int ii = 0; ii < graphs.size(); ++ii) + graphs[ii]->Draw("P"); legend->Draw(); - TLine *line = new TLine(1.0, 1.0, 20.0, 1.0); - line->SetLineStyle(2); line->SetLineWidth(2); line->SetLineColor(kRed); + TLine* line = new TLine(1.0, 1.0, 20.0, 1.0); + line->SetLineStyle(2); + line->SetLineWidth(2); + line->SetLineColor(kRed); line->Draw(); TPaveText* text = new TPaveText(0.575, 0.75, 0.97, 0.79, "brNDC"); - sprintf (name, "%s (%s)", HLT.c_str(), nameEtas[eta].c_str()); - if (debug) std::cout << "Name " << name << " |" << std::endl; + sprintf(name, "%s (%s)", HLT.c_str(), nameEtas[eta].c_str()); + if (debug) + std::cout << "Name " << name << " |" << std::endl; text->AddText(name); text->Draw("same"); TPaveText* text2 = new TPaveText(0.575, 0.71, 0.97, 0.75, "brNDC"); - sprintf (name, "CMS Preliminary"); + sprintf(name, cmsP.c_str()); text2->AddText(name); text2->Draw("same"); return canvas; } -TCanvas* plotEnergyPV(std::string fname, std::string HLT, int var, int ien, - int eta, bool varbin, int rebin, bool approve, bool logy, - int pos) { - - const int nPVBin=4; - int pvBins[nPVBin+1] = {1, 2, 3, 5, 100}; - TFile *file = TFile::Open(fname.c_str()); - char name[100]; - TObjArray histArr; +TCanvas* plotEnergyPV(std::string fname, + std::string HLT, + int var, + int ien, + int eta, + bool varbin, + int rebin, + bool approve, + bool logy, + int pos) { + const int nPVBin = 4; + int pvBins[nPVBin + 1] = {1, 2, 3, 5, 100}; + TFile* file = TFile::Open(fname.c_str()); + char name[100]; + TObjArray histArr; std::vector labels; - std::vector color; - double ymx0(0); - for (int i=4; iFindObjectAny(name); + std::vector color; + double ymx0(0); + for (int i = 4; i < nPVBin + 4; ++i) { + sprintf(name, "h_energy_%d_%d_%d_%d", i, ien, eta, var); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { histArr.AddLast(histo); - sprintf (name, "p=%s, #eta=%s, PV=%d:%d (%s)", varPs[ien].c_str(), varEta[eta].c_str(), pvBins[i-4], pvBins[i-3]-1, namefull[6].c_str()); + sprintf(name, + "p=%s, #eta=%s, PV=%d:%d (%s)", + varPs[ien].c_str(), + varEta[eta].c_str(), + pvBins[i - 4], + pvBins[i - 3] - 1, + namefull[6].c_str()); labels.push_back(name); - color.push_back(colors[i-4]); + color.push_back(colors[i - 4]); int ibin = histo->GetMaximumBin(); - if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin); + if (histo->GetBinContent(ibin) > ymx0) + ymx0 = histo->GetBinContent(ibin); } } - TCanvas *c(0); - if (histArr.GetEntries()>0) { + TCanvas* c(0); + if (histArr.GetEntries() > 0) { char cname[50]; - sprintf (cname, "c_%s_%d_%d", varEne1[var].c_str(), ien, eta); - sprintf ( name, "%s/p", varEne[var].c_str()); - c = plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, - 0.10, 0.05, 2.5, varbin, rebin, approve); + sprintf(cname, "c_%s_%d_%d", varEne1[var].c_str(), ien, eta); + sprintf(name, "%s/p", varEne[var].c_str()); + c = plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, 0.10, 0.05, 2.5, varbin, rebin, approve); } return c; } -TCanvas* plotTrack(std::string fname, std::string HLT, int var, bool varbin, - int rebin, bool approve, bool logy, int pos) { - - TFile *file = TFile::Open(fname.c_str()); - char name[100]; - TObjArray histArr; +TCanvas* plotTrack( + std::string fname, std::string HLT, int var, bool varbin, int rebin, bool approve, bool logy, int pos) { + TFile* file = TFile::Open(fname.c_str()); + char name[100]; + TObjArray histArr; std::vector labels; - std::vector color; - double ymx0(0); - for (int i=0; i<7; ++i) { - sprintf (name, "h_%s_%s", varname[var].c_str(), names[i].c_str()); - TH1D *histo = (TH1D*) file->FindObjectAny(name); + std::vector color; + double ymx0(0); + for (int i = 0; i < 7; ++i) { + sprintf(name, "h_%s_%s", varname[var].c_str(), names[i].c_str()); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { histArr.AddLast(histo); labels.push_back(namefull[i]); color.push_back(colors[i]); int ibin = histo->GetMaximumBin(); - if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin); + if (histo->GetBinContent(ibin) > ymx0) + ymx0 = histo->GetBinContent(ibin); } } - if (histArr.GetEntries()>0) { + if (histArr.GetEntries() > 0) { char cname[50]; - sprintf (cname, "c_%s", varname[var].c_str()); - sprintf ( name, "%s", vartitle[var].c_str()); - return plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, - 0.10, 0.05, -1, varbin, rebin, approve); + sprintf(cname, "c_%s", varname[var].c_str()); + sprintf(name, "%s", vartitle[var].c_str()); + return plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, 0.10, 0.05, -1, varbin, rebin, approve); } else { return 0; } } -TCanvas* plotIsolation(std::string fname, std::string HLT, int var, bool varbin, - int rebin, bool approve, bool logy, int pos) { - - TFile *file = TFile::Open(fname.c_str()); - char name[100]; - TObjArray histArr; +TCanvas* plotIsolation( + std::string fname, std::string HLT, int var, bool varbin, int rebin, bool approve, bool logy, int pos) { + TFile* file = TFile::Open(fname.c_str()); + char name[100]; + TObjArray histArr; std::vector labels; - std::vector color; - double ymx0(0); - for (int i=0; i<2; ++i) { - sprintf (name, "h_%s_%s", varnameC[var].c_str(), nameC[i].c_str()); - TH1D *histo = (TH1D*) file->FindObjectAny(name); + std::vector color; + double ymx0(0); + for (int i = 0; i < 2; ++i) { + sprintf(name, "h_%s_%s", varnameC[var].c_str(), nameC[i].c_str()); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { histArr.AddLast(histo); labels.push_back(nameCF[i]); color.push_back(colors[i]); int ibin = histo->GetMaximumBin(); - if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin); + if (histo->GetBinContent(ibin) > ymx0) + ymx0 = histo->GetBinContent(ibin); } } - if (histArr.GetEntries()>0) { + if (histArr.GetEntries() > 0) { char cname[50]; - sprintf (cname, "c_%s", varnameC[var].c_str()); - sprintf ( name, "%s (GeV)", vartitlC[var].c_str()); - return plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, - 0.10, 0.05, -1, varbin, rebin, approve); + sprintf(cname, "c_%s", varnameC[var].c_str()); + sprintf(name, "%s (GeV)", vartitlC[var].c_str()); + return plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, 0.10, 0.05, -1, varbin, rebin, approve); } else { return 0; } } -TCanvas* plotHLT(std::string fname, std::string HLT, int run, bool varbin, - int rebin, bool approve, bool logy, int pos) { - - TFile *file = TFile::Open(fname.c_str()); - char name[100]; - TObjArray histArr; +TCanvas* plotHLT(std::string fname, std::string HLT, int run, bool varbin, int rebin, bool approve, bool logy, int pos) { + TFile* file = TFile::Open(fname.c_str()); + char name[100]; + TObjArray histArr; std::vector labels; - std::vector color; - double ymx0(0); - if (run > 0) sprintf (name, "h_HLTAccepts_%d", run); - else sprintf (name, "h_HLTAccept"); - TH1D *histo = (TH1D*) file->FindObjectAny(name); + std::vector color; + double ymx0(0); + if (run > 0) + sprintf(name, "h_HLTAccepts_%d", run); + else + sprintf(name, "h_HLTAccept"); + TH1D* histo = (TH1D*)file->FindObjectAny(name); if (histo) { histArr.AddLast(histo); labels.push_back(HLT); @@ -1172,93 +1955,117 @@ TCanvas* plotHLT(std::string fname, std::string HLT, int run, bool varbin, int ibin = histo->GetMaximumBin(); ymx0 = histo->GetBinContent(ibin); } - if (histArr.GetEntries()>0) { + if (histArr.GetEntries() > 0) { char cname[50], hname[50]; - if (run > 0) {sprintf (cname,"c_HLT_%d",run); sprintf (name,"Run %d",run); - } else {sprintf (cname,"c_HLTs"); sprintf (name, "All runs");} - sprintf (hname, " "); - return plotHisto(cname, "", histArr, labels, color, hname, ymx0, logy, pos, - 0.40, 0.01, -1, varbin, rebin, approve); + if (run > 0) { + sprintf(cname, "c_HLT_%d", run); + sprintf(name, "Run %d", run); + } else { + sprintf(cname, "c_HLTs"); + sprintf(name, "All runs"); + } + sprintf(hname, " "); + return plotHisto(cname, "", histArr, labels, color, hname, ymx0, logy, pos, 0.40, 0.01, -1, varbin, rebin, approve); } else { return 0; } } -TCanvas* plotHisto(char* cname, std::string HLT, TObjArray& histArr, - std::vector& labels, std::vector& color, - char* name, double ymx0, bool logy, int pos, double yloff, - double yhoff, double xmax, bool varbin, int rebin0, - bool approve) { - +TCanvas* plotHisto(char* cname, + std::string HLT, + TObjArray& histArr, + std::vector& labels, + std::vector& color, + char* name, + double ymx0, + bool logy, + int pos, + double yloff, + double yhoff, + double xmax, + bool varbin, + int rebin0, + bool approve) { int nentry = histArr.GetEntries(); double ymax = 10.; - for (int i=0; i<10; ++i) { - if (ymx0 < ymax) break; + for (int i = 0; i < 10; ++i) { + if (ymx0 < ymax) + break; ymax *= 10.; } - double ystep = ymax*0.1; - for (int i=0; i<9; ++i) { - if (ymax-ystep < ymx0) break; + double ystep = ymax * 0.1; + for (int i = 0; i < 9; ++i) { + if (ymax - ystep < ymx0) + break; ymax -= ystep; } double ymin(0); - if (logy) ymin = 1; -// std::cout << "ymin " << ymin << " ymax " << ymx0 << ":" << ymax << std::endl; + if (logy) + ymin = 1; - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptTitle(kFALSE); gStyle->SetPadBorderMode(0); gStyle->SetCanvasBorderMode(0); - if (approve) gStyle->SetOptStat(0); - else gStyle->SetOptStat(1110); - TCanvas *canvas = new TCanvas(cname, cname, 500, 500); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(kFALSE); + gStyle->SetPadBorderMode(0); + gStyle->SetCanvasBorderMode(0); + if (approve) + gStyle->SetOptStat(0); + else + gStyle->SetOptStat(1110); + TCanvas* canvas = new TCanvas(cname, cname, 500, 500); gPad->SetTopMargin(yhoff); - gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); gPad->SetBottomMargin(yloff); - if (logy) canvas->SetLogy(); + if (logy) + canvas->SetLogy(); double height = 0.08; - double dx = (nentry > 2) ? 0.50 : 0.30; - double dy = (nentry > 2) ? 0.12 : 0.04; - double dy2 = 0.035; - double xmin1 = (pos > 1) ? 0.375 : 0.75-dx; + double dx = (nentry > 2) ? 0.50 : 0.30; + double dy = (nentry > 2) ? 0.12 : 0.04; + double dy2 = 0.035; + double xmin1 = (pos > 1) ? 0.375 : 0.75 - dx; double xmin2 = (pos > 1) ? 0.12 : 0.65; double xmin3 = (pos > 1) ? 0.15 : 0.75; - double ymin1 = (pos%2 == 0) ? (1.0-yhoff-dy) : (yloff+0.02); - double ymin2 = (pos%2 == 0) ? (0.96-yhoff-nentry*height) : (yloff+0.025+nentry*height); - double ymin3 = (pos%2 == 0) ? (ymin2-dy2) : (ymin2+dy2); - double dx2 = (approve) ? dx : 0.32; - double dx3 = (approve) ? dx : 0.22; + double ymin1 = (pos % 2 == 0) ? (1.0 - yhoff - dy) : (yloff + 0.02); + double ymin2 = (pos % 2 == 0) ? (0.96 - yhoff - nentry * height) : (yloff + 0.025 + nentry * height); + double ymin3 = (pos % 2 == 0) ? (ymin2 - dy2) : (ymin2 + dy2); + double dx2 = (approve) ? dx : 0.32; + double dx3 = (approve) ? dx : 0.22; if (approve) { - xmin1 = xmin2 = xmin3 = 0.975-dx; - ymin1 = (1.0-yhoff-dy); - ymin2 = ymin1-dy2-0.01; - ymin3 = ymin2-dy2; + xmin1 = xmin2 = xmin3 = 0.975 - dx; + ymin1 = (1.0 - yhoff - dy); + ymin2 = ymin1 - dy2 - 0.01; + ymin3 = ymin2 - dy2; } -// std::cout << nentry << " " << pos << " " << dx << " " << dy << " " << xmin1 << " " << xmin2 << std::endl; - TLegend *legend = new TLegend(xmin1, ymin1, xmin1+dx, ymin1+dy); + TLegend* legend = new TLegend(xmin1, ymin1, xmin1 + dx, ymin1 + dy); TPaveText *text, *text2; if (varbin || rebin0 == 1) { - text = new TPaveText(xmin2, ymin2, xmin2+dx2, ymin2+dy2, "brNDC"); - text2 = new TPaveText(xmin3, ymin3, xmin3+dx3, ymin3+dy2, "brNDC"); + text = new TPaveText(xmin2, ymin2, xmin2 + dx2, ymin2 + dy2, "brNDC"); + text2 = new TPaveText(xmin3, ymin3, xmin3 + dx3, ymin3 + dy2, "brNDC"); } else { - text = new TPaveText(0.10, 0.95, dx2+0.10, dy2+0.95, "brNDC"); - text2 = new TPaveText(dx2+0.10, 0.95, dx2+dx3+0.10, dy2+0.95, "brNDC"); + text = new TPaveText(0.10, 0.95, dx2 + 0.10, dy2 + 0.95, "brNDC"); + text2 = new TPaveText(dx2 + 0.10, 0.95, dx2 + dx3 + 0.10, dy2 + 0.95, "brNDC"); } - legend->SetBorderSize(1); legend->SetFillColor(kWhite); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); char texts[200]; - sprintf (texts, "CMS Preliminary"); + sprintf(texts, cmsP.c_str()); text2->AddText(texts); - THStack *Hs = new THStack("hs2"," "); - for (int i=0; iRebin(rebin0); + THStack* Hs = new THStack("hs2", " "); + for (int i = 0; i < nentry; i++) { + TH1D* h = (varbin) ? rebin((TH1D*)histArr[i], i) : (TH1D*)((TH1D*)histArr[i])->Rebin(rebin0); h->SetLineColor(color[i]); + h->SetLineStyle(lstyle[i]); h->SetLineWidth(2); - h->SetMarkerSize(0.8); - double ymax0 = (varbin) ? ymax : rebin0*ymax; - h->GetYaxis()->SetRangeUser(ymin,ymax0); - if (xmax > 0 && (!varbin)) h->GetXaxis()->SetRangeUser(0,xmax); + h->SetMarkerSize(1.0); + double ymax0 = (varbin) ? ymax : rebin0 * ymax; + h->GetYaxis()->SetRangeUser(ymin, ymax0); + if (xmax > 0 && (!varbin)) + h->GetXaxis()->SetRangeUser(0, xmax); Hs->Add(h, "hist sames"); - legend->AddEntry(h,labels[i].c_str(),"l"); + legend->AddEntry(h, labels[i].c_str(), "l"); } Hs->Draw("nostack"); canvas->Update(); @@ -1269,29 +2076,34 @@ TCanvas* plotHisto(char* cname, std::string HLT, TObjArray& histArr, Hs->GetHistogram()->GetYaxis()->SetTitle("Tracks/0.01"); } else { Hs->GetHistogram()->GetYaxis()->SetTitle("Tracks"); - if (xmax > 0) Hs->GetHistogram()->GetXaxis()->SetRangeUser(0,xmax); + if (xmax > 0) + Hs->GetHistogram()->GetXaxis()->SetRangeUser(0, xmax); } canvas->Modified(); - + canvas->Update(); if (!approve) { - for (int i=0; iGetListOfFunctions()->FindObject("stats"); - if (st1 != NULL) { - if (pos%2 == 0) { - st1->SetY1NDC(1.0-yhoff-(i+1)*height); st1->SetY2NDC(1.0-yhoff-i*height); - } else { - st1->SetY1NDC(yloff+0.02+i*height); st1->SetY2NDC(yloff+0.02+(i+1)*height); - } - if (pos > 1) { - st1->SetX1NDC(0.15); st1->SetX2NDC(.375); - } else { - st1->SetX1NDC(0.75); st1->SetX2NDC(.975); - } - st1->SetTextColor(colors[i]); - } + TPaveStats* st1 = (TPaveStats*)h->GetListOfFunctions()->FindObject("stats"); + if (st1 != NULL) { + if (pos % 2 == 0) { + st1->SetY1NDC(1.0 - yhoff - (i + 1) * height); + st1->SetY2NDC(1.0 - yhoff - i * height); + } else { + st1->SetY1NDC(yloff + 0.02 + i * height); + st1->SetY2NDC(yloff + 0.02 + (i + 1) * height); + } + if (pos > 1) { + st1->SetX1NDC(0.15); + st1->SetX2NDC(.375); + } else { + st1->SetX1NDC(0.75); + st1->SetX2NDC(.975); + } + st1->SetTextColor(colors[i]); + } } } } @@ -1301,87 +2113,89 @@ TCanvas* plotHisto(char* cname, std::string HLT, TObjArray& histArr, text->Draw("same"); } text2->Draw("same"); - + return canvas; } -void getHistStats(TH1D *h, int& entries, int& integral, double& mean, double& meanE, double& rms, double& rmsE, int& uflow, int& oflow) { - entries = h->GetEntries(); +void getHistStats(TH1D* h, + int& entries, + int& integral, + double& mean, + double& meanE, + double& rms, + double& rmsE, + int& uflow, + int& oflow) { + entries = h->GetEntries(); integral = h->Integral(); - mean = h->GetMean(); - meanE = h->GetMeanError(); - rms = h->GetRMS(); - rmsE = h->GetRMSError(); - uflow = h->GetBinContent(0); - oflow = h->GetBinContent( h->GetNbinsX()+1 ); + mean = h->GetMean(); + meanE = h->GetMeanError(); + rms = h->GetRMS(); + rmsE = h->GetRMSError(); + uflow = h->GetBinContent(0); + oflow = h->GetBinContent(h->GetNbinsX() + 1); } -TFitResultPtr getHistFitStats(TH1F *h, const char* formula, double xlow, - double xup, unsigned int& nPar, double* par, - double* epar) { - - TFitResultPtr fit = h->Fit(formula, "+qRB0", "", xlow, xup); +TFitResultPtr getHistFitStats( + TH1F* h, const char* formula, double xlow, double xup, unsigned int& nPar, double* par, double* epar) { + TFitResultPtr fit = h->Fit(formula, "+qRB0", "", xlow, xup); nPar = fit->NPar(); const std::vector errors = fit->Errors(); - for (unsigned int i=0; iValue(i); + for (unsigned int i = 0; i < nPar; i++) { + par[i] = fit->Value(i); epar[i] = errors[i]; } return fit; } -void setHistAttr(TH1F *h, int icol, int lwid, int ltype) { +void setHistAttr(TH1F* h, int icol, int lwid, int ltype) { h->SetLineColor(icol); h->SetLineStyle(ltype); h->SetLineWidth(lwid); - TF1 *f = h->GetFunction("gaus"); + TF1* f = h->GetFunction("gaus"); if (!f->IsZombie()) { f->SetLineColor(icol); f->SetLineStyle(2); } } -double getWeightedMean (int npt, int Start, std::vector& mean, std::vector& emean) { - double sumDen=0, sumNum=0; - for (int i=Start; i& mean, std::vector& emean) { + double sumDen = 0, sumNum = 0; + for (int i = Start; i < npt; i++) { + if (mean[i] == 0.0 || emean[i] == 0.0) { sumNum += 0; sumDen += 0; - } else { - sumNum += mean[i]/(emean[i]*emean[i]); - sumDen += 1.0/(emean[i]*emean[i]); + } else { + sumNum += mean[i] / (emean[i] * emean[i]); + sumDen += 1.0 / (emean[i] * emean[i]); } } - double WeightedMean = sumNum/sumDen; + double WeightedMean = sumNum / sumDen; return WeightedMean; } TH1D* rebin(TH1D* histin, int indx) { - std::string nameIn(histin->GetName()); char name[200]; - sprintf (name, "%sRebin%d", nameIn.c_str(), indx); - TH1D* hist = new TH1D(name,histin->GetXaxis()->GetTitle(),nbins,xbins); + sprintf(name, "%sRebin%d", nameIn.c_str(), indx); + TH1D* hist = new TH1D(name, histin->GetXaxis()->GetTitle(), nbins, xbins); std::vector cont; - for (int i=1; i<=histin->GetNbinsX(); ++i) { + for (int i = 1; i <= histin->GetNbinsX(); ++i) { double value = histin->GetBinContent(i); cont.push_back(value); -// std::cout << "Input Bin [" << i << "] = " << value << std::endl; } - for (int i=0; i0) totl /= kount; -// std::cout << "Output Bin [" << i+1 << "] Count " << kount << " Value " << totl << std::endl; - hist->SetBinContent(i+1,totl); + if (kount > 0) + totl /= kount; + hist->SetBinContent(i + 1, totl); } return hist; } diff --git a/Calibration/IsolatedParticles/macros/PlotVecGeom.C b/Calibration/IsolatedParticles/macros/PlotVecGeom.C index aa316ecc7835d..759560a1db37f 100644 --- a/Calibration/IsolatedParticles/macros/PlotVecGeom.C +++ b/Calibration/IsolatedParticles/macros/PlotVecGeom.C @@ -1,3 +1,31 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// plotEffiAll(approve, ratio, savePlot); +// Plots reconstructionefficiency as a function of p, pt, eta, phi +// +// approve (bool) To include "CMS Preliminary" or not (true) +// ratio (bool) If the ratio with the first input to be plotted (true) +// savePlot (int) Saving the plot or not: (no:gif:jpg:pdf -1:0:1:2) (2) +// +// plotCompareAll(cdir1, cdir2, cvers, cfil1, cfil2, ctype1, ctype2, +// postfix, ratio, logy, save, norm); +// Plots comparison plots of calorimetric quatities from native/vecgeom +// versions +// +// cdir1 (std::string) Name of the native directory ("10.7.r06.g4") +// cdir2 (std::string) Name of the vecgeom directory ("10.7.r06.vg") +// cvers (std::string) Name of the input version ("10.7.r06 MinBias") +// cfil1 (std::string) Name of the first input file ("minbias.root") +// cfil2 (std::string) Name of the second input file ("minbias.root") +// ctype1 (std::string) Name of the first type ("Native") +// ctype2 (std::string) Name of the second type ("VecGeom v1.1.16") +// postfix (std::string) Tag to be given to the canvas name ("MBR") +// ratio (bool) If the ratio to be plotted (false) +// logy (bool) If the scale is log or linear (true) +// save (bool) If the canvas is to be saved (true) +// norm (bool) If the plots to be normalized (true) +// +/////////////////////////////////////////////////////////////////////////////// #include "TCanvas.h" #include "TDirectory.h" #include "TF1.h" @@ -22,35 +50,80 @@ #include #include +TH1D* getEffi(TFile* file, std::string varname, unsigned int ifl); +TCanvas* plotEffi(int type, bool approve); +TCanvas* plotEffiRatio(int type, bool approve); +void plotEffiAll(bool approve = true, bool ratio = true, int savePlot = 2); +void plotCompare(const char* infile1, + const char* text1, + const char* infile2, + const char* text2, + int type1 = -1, + int type2 = -1, + int type3 = -1, + bool ratio = false, + bool logy = true, + std::string postfix = "", + bool save = false, + bool normalize = false); +void plotCompareAll(std::string cdir1 = "11.0.r06.g4", + std::string cdir2 = "11.0.r06.vg", + std::string cvers = "11.0.r06 MinBias", + std::string cfil1 = "minbias.root", + std::string cfil2 = "minbias.root", + std::string ctype1 = "Native", + std::string ctype2 = "VecGeom v1.1.20", + std::string postfix = "MBR", + bool ratio = false, + bool logy = true, + bool save = true, + bool norm = true); -int markStyle[7] = {20, 21, 24, 22, 23, 33, 25}; -int colors[7] = {1, 2, 4, 6, 7, 38, 3}; -int lineStyle[7] = {1, 2, 3, 4, 1, 2, 3}; +int markStyle[7] = {23, 22, 24, 25, 21, 33, 20}; +int colors[7] = {2, 4, 6, 7, 46, 3, 1}; +int lineStyle[7] = {1, 2, 3, 4, 1, 2, 3}; -const unsigned int nmodels=2; -//std::string filem[nmodels]={"pikp/FBE4r00MixStudyHLT.root","pikp/FBE4r00vMixStudyHLT.root"}; -std::string filem[nmodels]={"pikp/FBE4cMixStudyHLT10.root","pikp/FBE4vcMixStudyHLT10.root"}; -std::string typem[nmodels]={"10.4 FTFP_BERT_EMM (Native)","10.4 FTFP_BERT_EMM (VecGeom v0.5)"}; +/* +const unsigned int nmodels = 4; +std::string filem[nmodels] = {"pikp/FBE7r00vMixStudyHLT.root", + "pikp/FBE7p01vMixStudyHLT.root", + "pikp/FBE7p02vMixStudyHLT.root", + "pikp/FBE7p03vMixStudyHLT.root"}; +std::string typem[nmodels] = {"10.7 FTFP_BERT_EMM", + "10.7.p01 FTFP_BERT_EMM", + "10.7.p02 FTFP_BERT_EMM", + "10.7.p03 FTFP_BERT_EMM"}; +*/ +const unsigned int nmodels = 5; +std::string filem[nmodels] = {"pikp/FBE4p3vMixStudyHLT.root", + "pikp/FBE7p02vMixStudyHLT.root", + "pikp/FBE110p01vMixStudyHLT.root", + "pikp/FBE110r04MixStudyHLT.root", + "pikp/FBE110r06vMixStudyHLT.root"}; +std::string typem[nmodels] = {"10.4.p03 FTFP_BERT_EMM", + "10.7.p01 FTFP_BERT_EMM", + "11.0.p01 FTFP_BERT_EMM", + "11.0.ref04 FTFP_BERT_EMM", + "11.0.ref06 FTFP_BERT_EMM"}; TH1D* getEffi(TFile* file, std::string varname, unsigned int ifl) { - char name[100]; sprintf(name, "h_%s_All_0", varname.c_str()); - TH1D *hist1 = (TH1D*) file->FindObjectAny(name); + TH1D* hist1 = (TH1D*)file->FindObjectAny(name); sprintf(name, "h_%s_All_1", varname.c_str()); - TH1D *hist2 = (TH1D*) file->FindObjectAny(name); + TH1D* hist2 = (TH1D*)file->FindObjectAny(name); if (hist1 && hist2) { sprintf(name, "h_%s_Effy_%d", varname.c_str(), ifl); - int nbins = hist1->GetNbinsX(); - double xl = hist1->GetBinLowEdge(1); - double xh = hist1->GetBinLowEdge(nbins) + hist1->GetBinWidth(nbins); - TH1D* hist = new TH1D(name,hist1->GetTitle(),nbins,xl,xh); - for (int i=1; iGetNbinsX(); + double xl = hist1->GetBinLowEdge(1); + double xh = hist1->GetBinLowEdge(nbins) + hist1->GetBinWidth(nbins); + TH1D* hist = new TH1D(name, hist1->GetTitle(), nbins, xl, xh); + for (int i = 1; i < nbins; ++i) { double den = hist1->GetBinContent(i); - double val = (den > 0) ? (hist2->GetBinContent(i))/den : 0; - double err = (den > 0) ? (hist1->GetBinError(i))*(val/den) : 0; - hist->SetBinContent(i,val); - hist->SetBinError(i,err); + double val = (den > 0) ? (hist2->GetBinContent(i)) / den : 0; + double err = (den > 0) ? (hist1->GetBinError(i)) * (val / den) : 0; + hist->SetBinContent(i, val); + hist->SetBinError(i, err); } return hist; } else { @@ -59,52 +132,60 @@ TH1D* getEffi(TFile* file, std::string varname, unsigned int ifl) { } TCanvas* plotEffi(int type, bool approve) { - - std::string varnam[4] = {"pt","p","eta","phi"}; + std::string varnam[4] = {"pt", "p", "eta", "phi"}; std::string xvtitl[4] = {"p_{T} (GeV)", "p (GeV)", "#eta", "#phi"}; - bool irng[4] = {true, true, true, false}; - double xlowr[4] = { 0.0, 0.0, -2.2, -3.1415926}; - double xtopr[4] = {20.0, 20.0, 2.2, 3.1415926}; + bool irng[4] = {true, true, true, false}; + double xlowr[4] = {0.0, 0.0, -2.2, -3.1415926}; + double xtopr[4] = {20.0, 20.0, 2.2, 3.1415926}; TCanvas* c(0); - if (type < 0 || type > 3) type = 0; - TObjArray histArr; - for (unsigned k=0; k 3) + type = 0; + TObjArray histArr; + for (unsigned k = 0; k < nmodels; ++k) { + TFile* file = TFile::Open(filem[k].c_str()); + TH1D* hist = getEffi(file, varnam[type], k); if (hist) { hist->GetXaxis()->SetTitle(xvtitl[type].c_str()); hist->GetYaxis()->SetTitle("Efficiency"); - if (irng[type]) hist->GetXaxis()->SetRangeUser(xlowr[type],xtopr[type]); + if (irng[type]) + hist->GetXaxis()->SetRangeUser(xlowr[type], xtopr[type]); histArr.AddLast(hist); } } - if (histArr.GetEntries()>0) { - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptTitle(kFALSE); gStyle->SetPadBorderMode(0); - gStyle->SetCanvasBorderMode(0); gStyle->SetOptStat(0); + if (histArr.GetEntries() > 0) { + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(kFALSE); + gStyle->SetPadBorderMode(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetOptStat(0); char cname[50]; - sprintf (cname, "c_%sEff", varnam[type].c_str()); + sprintf(cname, "c_%sEff", varnam[type].c_str()); c = new TCanvas(cname, cname, 500, 500); - gPad->SetTopMargin(0.10); gPad->SetBottomMargin(0.10); - gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025); + gPad->SetTopMargin(0.10); + gPad->SetBottomMargin(0.10); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); - TLegend *legend = new TLegend(0.30, 0.15, 0.975, 0.30); - TPaveText *text1 = new TPaveText(0.05, 0.94, 0.35, 0.99, "brNDC"); - legend->SetBorderSize(1); legend->SetFillColor(kWhite); + TLegend* legend = new TLegend(0.30, 0.15, 0.975, 0.30); + TPaveText* text1 = new TPaveText(0.05, 0.94, 0.35, 0.99, "brNDC"); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); char texts[200]; - sprintf (texts, "CMS Preliminary"); + sprintf(texts, "CMS Preliminary"); text1->AddText(texts); - THStack *Hs = new THStack("hs2"," "); - for (int i=0; iSetLineColor(colors[i]); h->SetLineWidth(2); h->SetMarkerSize(markStyle[i]); Hs->Add(h, "hist sames"); - legend->AddEntry(h,typem[i].c_str(),"l"); + legend->AddEntry(h, typem[i].c_str(), "l"); } Hs->Draw("nostack"); c->Update(); @@ -112,148 +193,380 @@ TCanvas* plotEffi(int type, bool approve) { Hs->GetHistogram()->GetXaxis()->SetLabelSize(0.035); Hs->GetHistogram()->GetYaxis()->SetTitleOffset(1.6); Hs->GetHistogram()->GetYaxis()->SetTitle("Track Reconstruction Efficiency"); - Hs->GetHistogram()->GetYaxis()->SetRangeUser(0.0,1.2); + Hs->GetHistogram()->GetYaxis()->SetRangeUser(0.0, 1.2); if (irng[type]) - Hs->GetHistogram()->GetXaxis()->SetRangeUser(xlowr[type],xtopr[type]); + Hs->GetHistogram()->GetXaxis()->SetRangeUser(xlowr[type], xtopr[type]); c->Modified(); c->Update(); legend->Draw(""); - if (approve) text1->Draw("same"); + if (approve) + text1->Draw("same"); c->Modified(); c->Update(); } return c; } -void plotEffiAll(bool approve=false, int savePlot=-1) { - for (int var=0; var<=4; ++var) { - TCanvas* c = plotEffi(var, approve); +TCanvas* plotEffiRatio(int type, bool approve) { + std::string varnam[4] = {"pt", "p", "eta", "phi"}; + std::string xvtitl[4] = {"p_{T} (GeV)", "p (GeV)", "#eta", "#phi"}; + bool irng[4] = {true, true, true, false}; + double xlowr[4] = {0.0, 0.0, -2.2, -3.1415926}; + double xtopr[4] = {20.0, 20.0, 2.2, 3.1415926}; + + TCanvas* c(0); + if (type < 0 || type > 3) + type = 0; + TObjArray histArr; + TH1D* hist0(0); + for (unsigned k = 0; k < nmodels; ++k) { + TFile* file = TFile::Open(filem[k].c_str()); + TH1D* hist = getEffi(file, varnam[type], k); + if (hist) { + if (hist0 == nullptr) { + hist0 = hist; + } else { + int nbin = hist->GetNbinsX(); + int npt1(0), npt2(0); + double sumNum(0), sumDen(0); + for (int i = 1; i < nbin; ++i) { + double val1 = hist->GetBinContent(i); + double err1 = (val1 > 0) ? ((hist->GetBinError(i)) / val1) : 0; + double val2 = hist0->GetBinContent(i); + double err2 = (val2 > 0) ? ((hist0->GetBinError(i)) / val2) : 0; + double ratio = (val2 > 0) ? (val1 / val2) : 0; + double drat = ratio * sqrt(err1 * err1 + err2 * err2); + if ((((hist->GetBinLowEdge(i)) >= xlowr[type]) && + ((hist->GetBinLowEdge(i) + hist->GetBinWidth(i)) <= xtopr[type])) || + (!irng[type])) { + ++npt1; + if (val2 > 0) { + double temp1 = (ratio > 1.0) ? 1.0 / ratio : ratio; + double temp2 = (ratio > 1.0) ? drat / (ratio * ratio) : drat; + sumNum += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen += (1.0 / (temp2 * temp2)); + ++npt2; + } + } + hist->SetBinContent(i, ratio); + hist->SetBinError(i, drat); + } + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + std::cout << "Get Ratio of mean for " << npt2 << ":" << npt1 << ":" << nbin << " points: Mean " << sumNum + << " +- " << sumDen << " from " << filem[k] << std::endl; + hist->GetXaxis()->SetTitle(xvtitl[type].c_str()); + hist->GetYaxis()->SetTitle("Efficiency Ratio"); + hist->GetYaxis()->SetRangeUser(0.8, 1.2); + if (irng[type]) + hist->GetXaxis()->SetRangeUser(xlowr[type], xtopr[type]); + histArr.AddLast(hist); + } + } + } + + if (histArr.GetEntries() > 0) { + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(kFALSE); + gStyle->SetPadBorderMode(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetOptStat(0); + + char cname[50]; + sprintf(cname, "c_%sEffRat", varnam[type].c_str()); + c = new TCanvas(cname, cname, 500, 500); + gPad->SetTopMargin(0.10); + gPad->SetBottomMargin(0.10); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); + + TLegend* legend = new TLegend(0.30, 0.15, 0.975, 0.30); + TPaveText* text1 = new TPaveText(0.05, 0.94, 0.35, 0.99, "brNDC"); + legend->SetBorderSize(1); + legend->SetFillColor(kWhite); + char texts[200]; + sprintf(texts, "CMS Preliminary"); + text1->AddText(texts); + THStack* Hs = new THStack("hs2", " "); + for (int i = 0; i < histArr.GetEntries(); i++) { + TH1D* h = (TH1D*)histArr[i]; + h->SetLineColor(colors[i + 1]); + h->SetLineWidth(2); + h->SetMarkerSize(markStyle[i + 1]); + Hs->Add(h, "hist sames"); + legend->AddEntry(h, typem[i + 1].c_str(), "l"); + } + Hs->Draw("nostack"); + c->Update(); + Hs->GetHistogram()->GetXaxis()->SetTitle(xvtitl[type].c_str()); + Hs->GetHistogram()->GetXaxis()->SetLabelSize(0.035); + Hs->GetHistogram()->GetYaxis()->SetTitleOffset(1.6); + Hs->GetHistogram()->GetYaxis()->SetTitle("Track Reconstruction Efficiency Ratio"); + Hs->GetHistogram()->GetYaxis()->SetRangeUser(0.8, 1.2); + if (irng[type]) + Hs->GetHistogram()->GetXaxis()->SetRangeUser(xlowr[type], xtopr[type]); + c->Modified(); + c->Update(); + legend->Draw(""); + if (approve) + text1->Draw("same"); + c->Modified(); + c->Update(); + } + return c; +} + +void plotEffiAll(bool approve, bool ratio, int savePlot) { + for (int var = 0; var <= 4; ++var) { + TCanvas* c = (ratio) ? plotEffiRatio(var, approve) : plotEffi(var, approve); if (c != 0 && savePlot >= 0 && savePlot < 3) { std::string ext[3] = {"eps", "gif", "pdf"}; char name[200]; - sprintf (name, "%s.%s", c->GetName(), ext[savePlot].c_str()); + sprintf(name, "%s.%s", c->GetName(), ext[savePlot].c_str()); c->Print(name); } } } -void plotCompare(const char* infile1, const char* text1, const char* infile2, - const char* text2, int type1=-1, int type2=-1, int type3=-1, - bool logy=true, bool save=false) { - - int ndets[4] = {1, 9, 9, 15}; - int types[4] = {7, 8, 2, 3}; - std::string htype0[7] = {"PtInc", "EneInc", "EtaInc", "PhiInc", - "HitLow", "HitHigh", "HitMu"}; - int rebin0[7] = {1, 1, 1, 1, 10, 10, 10}; - std::string htype1[8] = {"Hit", "Time", "Edep", "Etot", "TimeAll", "EdepEM", - "EdepHad", "EtotG"}; - int rebin1[8] = {10, 10, 100, 100, 10, 50, 50, 50}; +void plotCompare(const char* infile1, + const char* text1, + const char* infile2, + const char* text2, + int type1, + int type2, + int type3, + bool ratio, + bool logy, + std::string postfix, + bool save, + bool normalize) { + int ndets[4] = {1, 9, 9, 15}; + int types[4] = {7, 8, 2, 3}; + std::string htype0[7] = {"PtInc", "EneInc", "EtaInc", "PhiInc", "HitLow", "HitHigh", "HitMu"}; + int rebin0[7] = {1, 1, 1, 1, 10, 10, 10}; + std::string htype1[8] = {"Hit", "Time", "Edep", "Etot", "TimeAll", "EdepEM", "EdepHad", "EtotG"}; + int rebin1[8] = {10, 10, 100, 100, 10, 50, 50, 50}; std::string htype2[2] = {"EdepCal", "EdepCalT"}; - int rebin2[2] = {50, 50}; + int rebin2[2] = {50, 50}; std::string htype3[3] = {"HitTk", "TimeTk", "EdepTk"}; - int rebin3[3] = {10, 10, 50}; + int rebin3[3] = {10, 10, 50}; + const double ymin(10.0); - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptTitle(0); gStyle->SetOptFit(0); - gStyle->SetOptStat(1110); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(0); + gStyle->SetOptFit(0); + if (ratio || normalize) + gStyle->SetOptStat(0); + else + gStyle->SetOptStat(1110); int itmin1 = (type1 >= 0) ? type1 : 0; int itmax1 = (type1 >= 0) ? type1 : 3; - TFile *file1 = new TFile(infile1); - TFile *file2 = new TFile(infile2); - std::cout << "File1: " << infile1 << ":" << file1 << " File2: " << infile2 - << ":" << file2 << std::endl; + TFile* file1 = new TFile(infile1); + TFile* file2 = new TFile(infile2); + std::cout << "File1: " << infile1 << ":" << file1 << " File2: " << infile2 << ":" << file2 << std::endl; if (file1 != 0 && file2 != 0) { - for (int it1=itmin1; it1<=itmax1; ++it1) { + for (int it1 = itmin1; it1 <= itmax1; ++it1) { int itmin2 = (type2 >= 0) ? type2 : 0; - int itmax2 = (type2 >= 0) ? type2 : ((type1 == 1) ? 3 : types[it1]-1); + int itmax2 = (type2 >= 0) ? type2 : ((type1 == 1) ? 3 : types[it1] - 1); int itmin3 = (type3 >= 0) ? type3 : 0; - int itmax3 = (type3 >= 0) ? type3 : ndets[it1]-1; - for (int it2=itmin2; it2<=itmax2; ++it2) { - int rebin(1); - if (it1 == 0) rebin = rebin0[it2]; - else if (it1 == 1) rebin = rebin1[it2]; - else if (it1 == 2) rebin = rebin2[it2]; - else if (it1 == 3) rebin = rebin3[it2]; - for (int it3=itmin3; it3<=itmax3; ++it3) { - if (type1 == 1 && (it3 == 1 || it3 == 2)) continue; - char name[20], namec[22]; - if (it1 == 0) sprintf (name, "%s", htype0[it2].c_str()); - else if (it1 == 1) sprintf (name, "%s%d", htype1[it2].c_str(), it3); - else if (it1 == 2) sprintf (name, "%s%d", htype2[it2].c_str(), it3); - else sprintf (name, "%s%d", htype3[it2].c_str(), it3); - TH1D* hist[2]; - hist[0] = (TH1D*)file1->FindObjectAny(name); - hist[1] = (TH1D*)file2->FindObjectAny(name); - if (hist[0] != 0 && hist[1] != 0) { - sprintf (namec, "c_%s", name); - TCanvas *pad = new TCanvas(namec, namec, 500, 500); - pad->SetRightMargin(0.10); - pad->SetTopMargin(0.10); - if (logy) pad->SetLogy(); - TLegend *legend = new TLegend(0.12, 0.79, 0.64, 0.89); - legend->SetFillColor(kWhite); - double ymax(0.90), dy(0.08); - for (int ih=0; ih<2; ++ih) { - hist[ih]->GetXaxis()->SetTitle(hist[ih]->GetTitle()); - hist[ih]->SetMarkerStyle(markStyle[ih]); - hist[ih]->SetMarkerColor(colors[ih]); - hist[ih]->SetLineStyle(lineStyle[ih]); - hist[ih]->SetLineColor(colors[ih]); - hist[ih]->SetLineWidth(2); - hist[ih]->GetYaxis()->SetTitleOffset(1.20); - if (rebin > 1) hist[ih]->Rebin(rebin); - if (ih == 0) { - legend->AddEntry(hist[ih],text1,"lp"); - hist[ih]->Draw(); - } else { - legend->AddEntry(hist[ih],text2,"lp"); - hist[ih]->Draw("sames"); - } - pad->Update(); - TPaveStats* st1 = (TPaveStats*)hist[ih]->GetListOfFunctions()->FindObject("stats"); - if (st1 != NULL) { - st1->SetLineColor(colors[ih]); - st1->SetTextColor(colors[ih]); - st1->SetY1NDC(ymax-dy); st1->SetY2NDC(ymax); - st1->SetX1NDC(0.65); st1->SetX2NDC(0.90); - ymax -= dy; - } - pad->Update(); - } - legend->Draw("same"); - pad->Update(); - if (save) { - sprintf (name, "%s.pdf", pad->GetName()); - pad->Print(name); - } - } - } + int itmax3 = (type3 >= 0) ? type3 : ndets[it1] - 1; + for (int it2 = itmin2; it2 <= itmax2; ++it2) { + int rebin(1); + if (it1 == 0) + rebin = rebin0[it2]; + else if (it1 == 1) + rebin = rebin1[it2]; + else if (it1 == 2) + rebin = rebin2[it2]; + else if (it1 == 3) + rebin = rebin3[it2]; + for (int it3 = itmin3; it3 <= itmax3; ++it3) { + if (type1 == 1 && (it3 == 1 || it3 == 2)) + continue; + char name[20], namec[22]; + if (it1 == 0) + sprintf(name, "%s", htype0[it2].c_str()); + else if (it1 == 1) + sprintf(name, "%s%d", htype1[it2].c_str(), it3); + else if (it1 == 2) + sprintf(name, "%s%d", htype2[it2].c_str(), it3); + else + sprintf(name, "%s%d", htype3[it2].c_str(), it3); + TH1D* hist[2]; + hist[0] = (TH1D*)file1->FindObjectAny(name); + hist[1] = (TH1D*)file2->FindObjectAny(name); + std::cout << name << " Hist " << hist[0] << ":" << hist[1] << std::endl; + if (hist[0] != 0 && hist[1] != 0) { + sprintf(namec, "c_%s%s", name, postfix.c_str()); + TCanvas* pad = new TCanvas(namec, namec, 500, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + if (logy) + pad->SetLogy(); + double xedge = (ratio) ? 0.88 : 0.64; + TLegend* legend = new TLegend(0.12, 0.79, xedge, 0.89); + legend->SetFillColor(kWhite); + if (ratio) { + int nbin = hist[0]->GetNbinsX(); + int nbinX = nbin / rebin; + sprintf(namec, "%sR", name); + TH1D* hist0 = (TH1D*)gROOT->FindObjectAny(namec); + if (hist0) + hist0->Delete(); + double xlow = hist[0]->GetBinLowEdge(1); + double xhigh = hist[0]->GetBinLowEdge(nbin) + hist[0]->GetBinWidth(nbin); + hist0 = new TH1D(namec, hist[0]->GetTitle(), nbinX, xlow, xhigh); + int npt1(0), npt2(0); + double sumNum(0), sumDen(0); + xhigh = hist0->GetBinLowEdge(1) + hist0->GetBinWidth(1); + double fact = (hist[1]->GetEntries()) / (hist[0]->GetEntries()); + bool empty(false); + for (int i = 1; i < nbinX; ++i) { + double val1(0), err1(0), val2(0), err2(0); + for (int j = 0; j < rebin; ++j) { + int i1 = (i - 1) * rebin + j + 1; + val1 += hist[0]->GetBinContent(i1); + val2 += hist[1]->GetBinContent(i1); + err1 += ((hist[0]->GetBinError(i1)) * (hist[0]->GetBinError(i1))); + err2 += ((hist[1]->GetBinError(i1)) * (hist[1]->GetBinError(i1))); + } + err1 = (val1 > 0) ? (sqrt(err1) / val1) : 0; + err2 = (val2 > 0) ? (sqrt(err2) / val2) : 0; + double ratio = (((val1 > 0) && (val2 > 0)) ? fact * (val1 / val2) : 1); + double drat = (((val1 > 0) && (val2 > 0)) ? (ratio * sqrt(err1 * err1 + err2 * err2)) : 0); + ++npt1; + if ((val1 > ymin) && (val2 > ymin)) { + double temp1 = (ratio > 1.0) ? 1.0 / ratio : ratio; + double temp2 = (ratio > 1.0) ? drat / (ratio * ratio) : drat; + sumNum += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen += (1.0 / (temp2 * temp2)); + ++npt2; + } + hist0->SetBinContent(i, ratio); + hist0->SetBinError(i, drat); + if (i <= 10) + xhigh = (hist0->GetBinLowEdge(i) + hist0->GetBinWidth(i)); + else if (val1 > 0 && val2 > 0 && (!empty)) + xhigh = (hist0->GetBinLowEdge(i) + hist0->GetBinWidth(i)); + else if ((val1 <= 0 || val2 <= 0) && i > 10) + empty = true; + } + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + std::cout << "Get Ratio of mean for " << hist[0]->GetTitle() << " with " << npt2 << ":" << npt1 << ":" + << nbinX << " points: Mean " << sumNum << " +- " << sumDen << std::endl; + if (it1 == 0) + sprintf(namec, "Ratio for %s", htype0[it2].c_str()); + else if (it1 == 1) + sprintf(namec, "Ratio for %s", htype1[it2].c_str()); + else if (it1 == 2) + sprintf(namec, "Ratio for %s", htype2[it2].c_str()); + else + sprintf(namec, "Ratio for %s", htype3[it2].c_str()); + hist0->SetMarkerStyle(markStyle[0]); + hist0->SetMarkerColor(colors[0]); + hist0->SetLineStyle(lineStyle[0]); + hist0->SetLineColor(colors[0]); + hist0->GetXaxis()->SetTitle(hist[0]->GetTitle()); + hist0->GetXaxis()->SetRangeUser(xlow, xhigh); + hist0->GetYaxis()->SetTitle(namec); + hist0->GetYaxis()->SetRangeUser(0.0, 2.0); + hist0->GetYaxis()->SetTitleOffset(1.3); + sprintf(namec, "%s vs %s", text1, text2); + legend->AddEntry(hist0, namec, "lp"); + hist0->Draw(); + pad->Update(); + TLine* line = new TLine(xlow, 1.0, xhigh, 1.0); + line->SetLineWidth(2); + line->SetLineStyle(2); + line->Draw("same"); + TPaveText* text0 = new TPaveText(0.12, 0.12, 0.65, 0.17, "brNDC"); + char texts[200]; + sprintf(texts, "Mean Deviation = %5.3f #pm %5.3f", sumNum, sumDen); + text0->SetFillColor(kWhite); + text0->AddText(texts); + text0->Draw("same"); + } else { + double ymax(0.90), dy(0.08); + for (int ih = 0; ih < 2; ++ih) { + hist[ih]->GetXaxis()->SetTitle(hist[ih]->GetTitle()); + hist[ih]->SetMarkerStyle(markStyle[ih]); + hist[ih]->SetMarkerColor(colors[ih]); + hist[ih]->SetLineStyle(lineStyle[ih]); + hist[ih]->SetLineColor(colors[ih]); + hist[ih]->SetLineWidth(2); + hist[ih]->GetYaxis()->SetTitleOffset(1.20); + if (rebin > 1) + hist[ih]->Rebin(rebin); + if (ih == 0) { + legend->AddEntry(hist[ih], text1, "lp"); + if (normalize) + hist[ih]->DrawNormalized("hist"); + else + hist[ih]->Draw(); + } else { + legend->AddEntry(hist[ih], text2, "lp"); + if (normalize) + hist[ih]->DrawNormalized("sames hist"); + else + hist[ih]->Draw("sames"); + } + pad->Update(); + TPaveStats* st1 = (TPaveStats*)hist[ih]->GetListOfFunctions()->FindObject("stats"); + if (st1 != NULL) { + st1->SetLineColor(colors[ih]); + st1->SetTextColor(colors[ih]); + st1->SetY1NDC(ymax - dy); + st1->SetY2NDC(ymax); + st1->SetX1NDC(0.65); + st1->SetX2NDC(0.90); + ymax -= dy; + } + pad->Update(); + } + } + legend->Draw("same"); + pad->Update(); + if (save) { + sprintf(name, "%s.pdf", pad->GetName()); + pad->Print(name); + } + } + } } } } } -void plotCompareAll(const std::string & cdir1="10.4.r00.g4", - const std::string & cdir2="10.4.r00.vg", - const std::string & cvers="10.4 MinBias", - const std::string & cfile= "minbias.root", - const std::string & ctype1="Native", - const std::string & ctype2="VecGeom v0.5", - bool logy=true, bool save=false) { - +void plotCompareAll(std::string cdir1, + std::string cdir2, + std::string cvers, + std::string cfil1, + std::string cfil2, + std::string ctype1, + std::string ctype2, + std::string postfix, + bool ratio, + bool logy, + bool save, + bool norm) { char infile1[200], infile2[200], text1[200], text2[200]; - sprintf (infile1, "%s/%s", cdir1.c_str(), cfile.c_str()); - sprintf (infile2, "%s/%s", cdir2.c_str(), cfile.c_str()); - sprintf (text1, "%s (%s)", cvers.c_str(), ctype1.c_str()); - sprintf (text2, "%s (%s)", cvers.c_str(), ctype2.c_str()); - plotCompare(infile1,text1,infile2,text2,1,-1,0,logy,save); - plotCompare(infile1,text1,infile2,text2,1,-1,3,logy,save); - plotCompare(infile1,text1,infile2,text2,1,-1,4,logy,save); - plotCompare(infile1,text1,infile2,text2,1,-1,5,logy,save); - plotCompare(infile1,text1,infile2,text2,1,-1,6,logy,save); - plotCompare(infile1,text1,infile2,text2,1,-1,7,logy,save); - plotCompare(infile1,text1,infile2,text2,1,-1,8,logy,save); + sprintf(infile1, "%s/%s", cdir1.c_str(), cfil1.c_str()); + sprintf(infile2, "%s/%s", cdir2.c_str(), cfil2.c_str()); + sprintf(text1, "%s (%s)", cvers.c_str(), ctype1.c_str()); + sprintf(text2, "%s (%s)", cvers.c_str(), ctype2.c_str()); + plotCompare(infile1, text1, infile2, text2, 1, -1, 0, ratio, logy, postfix, save, norm); + plotCompare(infile1, text1, infile2, text2, 1, -1, 3, ratio, logy, postfix, save, norm); + plotCompare(infile1, text1, infile2, text2, 1, -1, 4, ratio, logy, postfix, save, norm); + plotCompare(infile1, text1, infile2, text2, 1, -1, 5, ratio, logy, postfix, save, norm); + plotCompare(infile1, text1, infile2, text2, 1, -1, 6, ratio, logy, postfix, save, norm); + plotCompare(infile1, text1, infile2, text2, 1, -1, 7, ratio, logy, postfix, save, norm); + plotCompare(infile1, text1, infile2, text2, 1, -1, 8, ratio, logy, postfix, save, norm); } diff --git a/Calibration/IsolatedParticles/macros/StudyHLT.C b/Calibration/IsolatedParticles/macros/StudyHLT.C index c216f94e99126..19ae60cf5f131 100644 --- a/Calibration/IsolatedParticles/macros/StudyHLT.C +++ b/Calibration/IsolatedParticles/macros/StudyHLT.C @@ -7,14 +7,14 @@ // t.Loop() // // where -// infile string Name of the input ROOT tree file +// infile string Name of the input ROOT tree file // outfile string Name of the output ROOT histogram file // dirname string Name of the directory ("StudyHLT") // treeName string Name of the tree ("testTree") // // In addition there are useful methods: // -// void GetPUWeight(mcFile, dataFile, type, dirName) +// void GetPUWeight(mcFile, dataFile, type, dirName) // Calculates PilUp weights using distribution of number of vertex // where // mcFile string Name of the input ROOT tree MC file @@ -44,101 +44,102 @@ #include class StudyHLT { -public : - TTree *fChain; //!pointer to the analyzed TTree or TChain - Int_t fCurrent; //!current Tree number in a TChain +public: + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain // Declaration of leaf types - Int_t tr_goodRun; - Int_t tr_goodPV; - Double_t tr_eventWeight; - std::vector *tr_TrigName; - std::vector *tr_TrkPt; - std::vector *tr_TrkP; - std::vector *tr_TrkEta; - std::vector *tr_TrkPhi; - std::vector *tr_TrkID; - std::vector *tr_MaxNearP31X31; - std::vector *tr_MaxNearHcalP7x7; - std::vector *tr_FE7x7P; - std::vector *tr_FE11x11P; - std::vector *tr_FE15x15P; - std::vector *tr_SE7x7P; - std::vector *tr_SE11x11P; - std::vector *tr_SE15x15P; - std::vector *tr_H7x7; - std::vector *tr_H5x5; - std::vector *tr_H3x3; - std::vector *tr_iEta; + Int_t tr_goodRun; + Int_t tr_goodPV; + Double_t tr_eventWeight; + std::vector *tr_TrigName; + std::vector *tr_TrkPt; + std::vector *tr_TrkP; + std::vector *tr_TrkEta; + std::vector *tr_TrkPhi; + std::vector *tr_TrkID; + std::vector *tr_MaxNearP31X31; + std::vector *tr_MaxNearHcalP7x7; + std::vector *tr_FE7x7P; + std::vector *tr_FE11x11P; + std::vector *tr_FE15x15P; + std::vector *tr_SE7x7P; + std::vector *tr_SE11x11P; + std::vector *tr_SE15x15P; + std::vector *tr_H7x7; + std::vector *tr_H5x5; + std::vector *tr_H3x3; + std::vector *tr_iEta; // List of branches - TBranch *b_tr_goodRun; //! - TBranch *b_tr_goodPV; //! - TBranch *b_tr_eventWeight; //! - TBranch *b_tr_TrigName; //! - TBranch *b_tr_TrkPt; //! - TBranch *b_tr_TrkP; //! - TBranch *b_tr_TrkEta; //! - TBranch *b_tr_TrkPhi; //! - TBranch *b_tr_TrkID; //! - TBranch *b_tr_MaxNearP31X31; //! - TBranch *b_tr_MaxNearHcalP7x7; //! - TBranch *b_tr_TrkQuality; //! - TBranch *b_tr_TrkokECAL; //! - TBranch *b_tr_TrkokHCAL; //! - TBranch *b_tr_FE7x7P; //! - TBranch *b_tr_FE11x11P; //! - TBranch *b_tr_FE15x15P; //! - TBranch *b_tr_SE7x7P; //! - TBranch *b_tr_SE11x11P; //! - TBranch *b_tr_SE15x15P; //! - TBranch *b_tr_H7x7; //! - TBranch *b_tr_H5x5; //! - TBranch *b_tr_H3x3; //! - TBranch *b_tr_iEta; //! + TBranch *b_tr_goodRun; //! + TBranch *b_tr_goodPV; //! + TBranch *b_tr_eventWeight; //! + TBranch *b_tr_TrigName; //! + TBranch *b_tr_TrkPt; //! + TBranch *b_tr_TrkP; //! + TBranch *b_tr_TrkEta; //! + TBranch *b_tr_TrkPhi; //! + TBranch *b_tr_TrkID; //! + TBranch *b_tr_MaxNearP31X31; //! + TBranch *b_tr_MaxNearHcalP7x7; //! + TBranch *b_tr_TrkQuality; //! + TBranch *b_tr_TrkokECAL; //! + TBranch *b_tr_TrkokHCAL; //! + TBranch *b_tr_FE7x7P; //! + TBranch *b_tr_FE11x11P; //! + TBranch *b_tr_FE15x15P; //! + TBranch *b_tr_SE7x7P; //! + TBranch *b_tr_SE11x11P; //! + TBranch *b_tr_SE15x15P; //! + TBranch *b_tr_H7x7; //! + TBranch *b_tr_H5x5; //! + TBranch *b_tr_H3x3; //! + TBranch *b_tr_iEta; //! - StudyHLT(std::string inFile, std::string outFile, - std::string dirnam="StudyHLT", std::string treeNam="testTree"); + StudyHLT(std::string inFile, std::string outFile, std::string dirnam = "StudyHLT", std::string treeNam = "testTree"); virtual ~StudyHLT(); - virtual Int_t Cut(Long64_t entry); - virtual Int_t GetEntry(Long64_t entry); + 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); + virtual void Init(TTree *tree); + virtual void Loop(); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); - std::string outFile_; + std::string outFile_; }; -StudyHLT::StudyHLT(std::string inFile, std::string outFile, std::string dirnam, - std::string treeNam) : outFile_(outFile) { - - TFile *file = new TFile(inFile.c_str()); - TDirectory *dir = (TDirectory*)file->FindObjectAny(dirnam.c_str()); - std::cout << inFile << " file " << file << " " << dirnam << " " << dir - << std::endl; - TTree *tree = (TTree*)dir->Get(treeNam.c_str()); +StudyHLT::StudyHLT(std::string inFile, std::string outFile, std::string dirnam, std::string treeNam) + : outFile_(outFile) { + TFile *file = new TFile(inFile.c_str()); + TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnam.c_str()); + std::cout << inFile << " file " << file << " " << dirnam << " " << dir << std::endl; + TTree *tree = (TTree *)dir->Get(treeNam.c_str()); std::cout << "Tree:" << treeNam << " at " << tree << std::endl; Init(tree); } StudyHLT::~StudyHLT() { - if (!fChain) return; - delete fChain->GetCurrentFile(); + if (!fChain) + return; + delete fChain->GetCurrentFile(); } Int_t StudyHLT::GetEntry(Long64_t entry) { // Read contents of entry. - if (!fChain) return 0; + if (!fChain) + return 0; return fChain->GetEntry(entry); } Long64_t StudyHLT::LoadTree(Long64_t entry) { // Set the environment to read one entry - if (!fChain) return -5; + if (!fChain) + return -5; Long64_t centry = fChain->LoadTree(entry); - if (centry < 0) return centry; + if (centry < 0) + return centry; if (fChain->GetTreeNumber() != fCurrent) { fCurrent = fChain->GetTreeNumber(); Notify(); @@ -161,7 +162,7 @@ void StudyHLT::Init(TTree *tree) { tr_TrkP = 0; tr_TrkEta = 0; tr_TrkPhi = 0; - tr_TrkID = 0; + tr_TrkID = 0; tr_MaxNearP31X31 = 0; tr_MaxNearHcalP7x7 = 0; tr_FE7x7P = 0; @@ -175,32 +176,33 @@ void StudyHLT::Init(TTree *tree) { tr_H7x7 = 0; // Set branch addresses and branch pointers - if (!tree) return; + if (!tree) + return; fChain = tree; fCurrent = -1; fChain->SetMakeClass(1); - fChain->SetBranchAddress("tr_goodRun", &tr_goodRun, &b_tr_goodRun); - fChain->SetBranchAddress("tr_goodPV", &tr_goodPV, &b_tr_goodPV); + fChain->SetBranchAddress("tr_goodRun", &tr_goodRun, &b_tr_goodRun); + fChain->SetBranchAddress("tr_goodPV", &tr_goodPV, &b_tr_goodPV); fChain->SetBranchAddress("tr_eventWeight", &tr_eventWeight, &b_tr_eventWeight); - fChain->SetBranchAddress("tr_TrigName", &tr_TrigName, &b_tr_TrigName); - fChain->SetBranchAddress("tr_TrkPt", &tr_TrkPt, &b_tr_TrkPt); - fChain->SetBranchAddress("tr_TrkP", &tr_TrkP, &b_tr_TrkP); - fChain->SetBranchAddress("tr_TrkEta", &tr_TrkEta, &b_tr_TrkEta); - fChain->SetBranchAddress("tr_TrkPhi", &tr_TrkPhi, &b_tr_TrkPhi); - fChain->SetBranchAddress("tr_TrkID", &tr_TrkID, &b_tr_TrkID); - fChain->SetBranchAddress("tr_MaxNearP31X31", &tr_MaxNearP31X31, &b_tr_MaxNearP31X31); + fChain->SetBranchAddress("tr_TrigName", &tr_TrigName, &b_tr_TrigName); + fChain->SetBranchAddress("tr_TrkPt", &tr_TrkPt, &b_tr_TrkPt); + fChain->SetBranchAddress("tr_TrkP", &tr_TrkP, &b_tr_TrkP); + fChain->SetBranchAddress("tr_TrkEta", &tr_TrkEta, &b_tr_TrkEta); + fChain->SetBranchAddress("tr_TrkPhi", &tr_TrkPhi, &b_tr_TrkPhi); + fChain->SetBranchAddress("tr_TrkID", &tr_TrkID, &b_tr_TrkID); + fChain->SetBranchAddress("tr_MaxNearP31X31", &tr_MaxNearP31X31, &b_tr_MaxNearP31X31); fChain->SetBranchAddress("tr_MaxNearHcalP7x7", &tr_MaxNearHcalP7x7, &b_tr_MaxNearHcalP7x7); - fChain->SetBranchAddress("tr_FE7x7P", &tr_FE7x7P, &b_tr_FE7x7P); - fChain->SetBranchAddress("tr_FE11x11P", &tr_FE11x11P, &b_tr_FE11x11P); - fChain->SetBranchAddress("tr_FE15x15P", &tr_FE15x15P, &b_tr_FE15x15P); - fChain->SetBranchAddress("tr_SE7x7P", &tr_SE7x7P, &b_tr_SE7x7P); - fChain->SetBranchAddress("tr_SE11x11P", &tr_SE11x11P, &b_tr_SE11x11P); - fChain->SetBranchAddress("tr_SE15x15P", &tr_SE15x15P, &b_tr_SE15x15P); - fChain->SetBranchAddress("tr_H3x3", &tr_H3x3, &b_tr_H3x3); - fChain->SetBranchAddress("tr_H5x5", &tr_H5x5, &b_tr_H5x5); - fChain->SetBranchAddress("tr_H7x7", &tr_H7x7, &b_tr_H7x7); - fChain->SetBranchAddress("tr_iEta", &tr_iEta, &b_tr_iEta); + fChain->SetBranchAddress("tr_FE7x7P", &tr_FE7x7P, &b_tr_FE7x7P); + fChain->SetBranchAddress("tr_FE11x11P", &tr_FE11x11P, &b_tr_FE11x11P); + fChain->SetBranchAddress("tr_FE15x15P", &tr_FE15x15P, &b_tr_FE15x15P); + fChain->SetBranchAddress("tr_SE7x7P", &tr_SE7x7P, &b_tr_SE7x7P); + fChain->SetBranchAddress("tr_SE11x11P", &tr_SE11x11P, &b_tr_SE11x11P); + fChain->SetBranchAddress("tr_SE15x15P", &tr_SE15x15P, &b_tr_SE15x15P); + fChain->SetBranchAddress("tr_H3x3", &tr_H3x3, &b_tr_H3x3); + fChain->SetBranchAddress("tr_H5x5", &tr_H5x5, &b_tr_H5x5); + fChain->SetBranchAddress("tr_H7x7", &tr_H7x7, &b_tr_H7x7); + fChain->SetBranchAddress("tr_iEta", &tr_iEta, &b_tr_iEta); Notify(); } @@ -210,18 +212,19 @@ Bool_t StudyHLT::Notify() { // 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 StudyHLT::Show(Long64_t entry) { // Print contents of entry. // If entry is not specified, print current entry - if (!fChain) return; + if (!fChain) + return; fChain->Show(entry); } -Int_t StudyHLT::Cut(Long64_t ) { +Int_t StudyHLT::Cut(Long64_t) { // This function may be called from Loop. // returns 1 if entry is accepted. // returns -1 otherwise. @@ -229,196 +232,218 @@ Int_t StudyHLT::Cut(Long64_t ) { } void StudyHLT::Loop() { - //create a new root file to store the output. TFile *f = new TFile(outFile_.c_str(), "RECREATE"); - if (fChain == 0) return; + if (fChain == 0) + return; //declare the histograms in the same way as is done in studyHLT - std::string titls[6] = {"NoIso","okEcal","EcalCharIso","HcalCharIso", - "EcalNeutIso","HcalNeutIso"}; + std::string titls[6] = {"NoIso", "okEcal", "EcalCharIso", "HcalCharIso", "EcalNeutIso", "HcalNeutIso"}; char name[40], htit[400]; - TH1D *h_pt[6], *h_p[6], *h_eta[6], *h_phi[6]; - for (int i=0; i<6; ++i) { - sprintf (name, "h_pt_%s", titls[i].c_str()); - h_pt[i] = new TH1D(name,"", 400, 0, 200); - sprintf (name, "h_p_%s", titls[i].c_str()); - h_p[i] = new TH1D(name,"", 400, 0, 200); - sprintf (name, "h_eta_%s", titls[i].c_str()); - h_eta[i] = new TH1D(name,"", 60, -3.0, 3.0); - sprintf (name, "h_phi_%s", titls[i].c_str()); - h_phi[i] = new TH1D(name,"", 100, -3.15, 3.15); + TH1D *h_pt[6], *h_p[6], *h_eta[6], *h_phi[6]; + for (int i = 0; i < 6; ++i) { + sprintf(name, "h_pt_%s", titls[i].c_str()); + h_pt[i] = new TH1D(name, "", 400, 0, 200); + sprintf(name, "h_p_%s", titls[i].c_str()); + h_p[i] = new TH1D(name, "", 400, 0, 200); + sprintf(name, "h_eta_%s", titls[i].c_str()); + h_eta[i] = new TH1D(name, "", 60, -3.0, 3.0); + sprintf(name, "h_phi_%s", titls[i].c_str()); + h_phi[i] = new TH1D(name, "", 100, -3.15, 3.15); } - static const int nPBin=10, nEtaBin=4, nPVBin=4; - TH1D *h_energy[nPVBin+1][nPBin][nEtaBin][6]; - double pBin[nPBin+1] = {1.0,2.0,3.0,4.0,5.0,6.0,7.0,9.0,11.0,15.0,20.0}; - int etaBin[nEtaBin+1] = {1, 7, 13, 17, 23}; - int pvBin[nPVBin+1] = {1, 2, 3, 5, 100}; - std::string energyNames[6]={"E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})", - "E_{11x11}", "H_{5x5}", "{E_{11x11}+H_{5x5})"}; - for (int i=0; iSumw2(); - } + static const int nPBin = 10, nEtaBin = 4, nPVBin = 4; + TH1D *h_energy[nPVBin + 1][nPBin][nEtaBin][6]; + double pBin[nPBin + 1] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 9.0, 11.0, 15.0, 20.0}; + int etaBin[nEtaBin + 1] = {1, 7, 13, 17, 23}; + int pvBin[nPVBin + 1] = {1, 2, 3, 5, 100}; + std::string energyNames[6] = { + "E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})", "E_{11x11}", "H_{5x5}", "{E_{11x11}+H_{5x5})"}; + for (int i = 0; i < nPVBin; ++i) { + for (int ip = 0; ip < nPBin; ++ip) { + for (int ie = 0; ie < nEtaBin; ++ie) { + for (int j = 0; j < 6; ++j) { + sprintf(name, "h_energy_%d_%d_%d_%d", i, ip, ie, j); + if (i < nPVBin) { + sprintf(htit, + "%s/p (p=%4.1f:%4.1f; i#eta=%d:%d, PV=%d:%d)", + energyNames[j].c_str(), + pBin[ip], + pBin[ip + 1], + etaBin[ie], + (etaBin[ie + 1] - 1), + pvBin[i], + pvBin[i + 1]); + } else { + sprintf(htit, + "%s/p (p=%4.1f:%4.1f; i#eta=%d:%d, All PV)", + energyNames[j].c_str(), + pBin[ip], + pBin[ip + 1], + etaBin[ie], + (etaBin[ie + 1] - 1)); + } + h_energy[i][ip][ie][j] = new TH1D(name, htit, 500, -0.1, 4.9); + h_energy[i][ip][ie][j]->Sumw2(); + } } } } Long64_t nentries = fChain->GetEntries(); Long64_t nbytes = 0, nb = 0; - for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; - if ((tr_TrkPt->size() != tr_TrkP->size()) || - (tr_TrkPt->size() != tr_TrkEta->size()) || - (tr_TrkPt->size() != tr_TrkPhi->size())) { + if (ientry < 0) + break; + nb = fChain->GetEntry(jentry); + nbytes += nb; + if ((tr_TrkPt->size() != tr_TrkP->size()) || (tr_TrkPt->size() != tr_TrkEta->size()) || + (tr_TrkPt->size() != tr_TrkPhi->size())) { std::cout << "Error processing samples " << std::endl; break; - }//matches if + } //matches if //fill the distributions //loop over all the reco tracks - for (unsigned int itk=0; itk !=tr_TrkPt->size(); ++itk){ - - h_pt[0]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight); - h_p[0]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight); - h_eta[0]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight); - h_phi[0]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight); - if ( (*tr_TrkPt)[itk] > 1.0 && abs((*tr_TrkEta)[itk]) < 2.5) { - h_pt[1]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight); - h_p[1]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight); - h_eta[1]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight); - h_phi[1]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight); - //condition of charged Isolation - if ( (*tr_MaxNearP31X31)[itk] < 0) { - h_pt[2]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight); - h_p[2]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight); - h_eta[2]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight); - h_phi[2]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight); - //condition for HCal Charged Isolation - if ( (*tr_MaxNearHcalP7x7)[itk] < 0) { - h_pt[3]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight); - h_p[3]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight); - h_eta[3]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight); - h_phi[3]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight); - //condition of Neutral Isolation - if ( (*tr_SE11x11P)[itk] && (*tr_SE15x15P)[itk] && - ( (*tr_FE15x15P)[itk] - (*tr_FE11x11P)[itk]) < 2.0 ) { - h_pt[4]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight); - h_p[4]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight); - h_eta[4]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight); - h_phi[4]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight); - if ( ((*tr_H7x7)[itk] - (*tr_H5x5)[itk] ) < 2.0) { - h_pt[5]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight); - h_p[5]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight); - h_eta[5]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight); - h_phi[5]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight); - int ip(-1), ie(-1), nPV(-1); - for (int i=0; i= pBin[i]) && - ((*tr_TrkP)[itk] < pBin[i+1])) { ip = i; break; } - } - for (int i=0; i= etaBin[i]) && - ((*tr_iEta)[itk] < etaBin[i+1])) { ie = i; break; } - } - for (int i=0; i= pvBin[i]) && - ((tr_goodPV) < pvBin[i+1])) { nPV = i; break; } - } - double den = 1.0/((*tr_TrkP)[itk]); - if ((ip >= 0) && (ie >= 0) && ((*tr_FE7x7P)[itk] > 0.02) - && ((*tr_H3x3)[itk] > 0.1) && (nPV >= 0)) { - h_energy[nPV][ip][ie][0]->Fill(den*(*tr_FE7x7P)[itk], tr_eventWeight); - h_energy[nPV][ip][ie][1]->Fill(den*(*tr_H3x3)[itk], tr_eventWeight); - h_energy[nPV][ip][ie][2]->Fill(den*((*tr_FE7x7P)[itk]+(*tr_H3x3)[itk]), tr_eventWeight); - h_energy[nPV][ip][ie][3]->Fill(den*(*tr_FE11x11P)[itk], tr_eventWeight); - h_energy[nPV][ip][ie][4]->Fill(den*(*tr_H5x5)[itk], tr_eventWeight); - h_energy[nPV][ip][ie][5]->Fill(den*((*tr_FE11x11P)[itk]+(*tr_H5x5)[itk]), tr_eventWeight); - h_energy[nPVBin][ip][ie][0]->Fill(den*(*tr_FE7x7P)[itk], tr_eventWeight); - h_energy[nPVBin][ip][ie][1]->Fill(den*(*tr_H3x3)[itk], tr_eventWeight); - h_energy[nPVBin][ip][ie][2]->Fill(den*((*tr_FE7x7P)[itk]+(*tr_H3x3)[itk]), tr_eventWeight); - h_energy[nPVBin][ip][ie][3]->Fill(den*(*tr_FE11x11P)[itk], tr_eventWeight); - h_energy[nPVBin][ip][ie][4]->Fill(den*(*tr_H5x5)[itk], tr_eventWeight); - h_energy[nPVBin][ip][ie][5]->Fill(den*((*tr_FE11x11P)[itk]+(*tr_H5x5)[itk]), tr_eventWeight); - } - }//HCal Neutral Iso - }//neutral isolation - }//HCal Charged Iso - }//charged Iso + for (unsigned int itk = 0; itk != tr_TrkPt->size(); ++itk) { + h_pt[0]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight); + h_p[0]->Fill(((*tr_TrkP)[itk]), tr_eventWeight); + h_eta[0]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight); + h_phi[0]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight); + if ((*tr_TrkPt)[itk] > 1.0 && abs((*tr_TrkEta)[itk]) < 2.5) { + h_pt[1]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight); + h_p[1]->Fill(((*tr_TrkP)[itk]), tr_eventWeight); + h_eta[1]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight); + h_phi[1]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight); + //condition of charged Isolation + if ((*tr_MaxNearP31X31)[itk] < 0) { + h_pt[2]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight); + h_p[2]->Fill(((*tr_TrkP)[itk]), tr_eventWeight); + h_eta[2]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight); + h_phi[2]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight); + //condition for HCal Charged Isolation + if ((*tr_MaxNearHcalP7x7)[itk] < 0) { + h_pt[3]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight); + h_p[3]->Fill(((*tr_TrkP)[itk]), tr_eventWeight); + h_eta[3]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight); + h_phi[3]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight); + //condition of Neutral Isolation + if ((*tr_SE11x11P)[itk] && (*tr_SE15x15P)[itk] && ((*tr_FE15x15P)[itk] - (*tr_FE11x11P)[itk]) < 2.0) { + h_pt[4]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight); + h_p[4]->Fill(((*tr_TrkP)[itk]), tr_eventWeight); + h_eta[4]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight); + h_phi[4]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight); + if (((*tr_H7x7)[itk] - (*tr_H5x5)[itk]) < 2.0) { + h_pt[5]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight); + h_p[5]->Fill(((*tr_TrkP)[itk]), tr_eventWeight); + h_eta[5]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight); + h_phi[5]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight); + int ip(-1), ie(-1), nPV(-1); + for (int i = 0; i < nPBin; ++i) { + if (((*tr_TrkP)[itk] >= pBin[i]) && ((*tr_TrkP)[itk] < pBin[i + 1])) { + ip = i; + break; + } + } + for (int i = 0; i < nEtaBin; ++i) { + if (((*tr_iEta)[itk] >= etaBin[i]) && ((*tr_iEta)[itk] < etaBin[i + 1])) { + ie = i; + break; + } + } + for (int i = 0; i < nPVBin; ++i) { + if (((tr_goodPV) >= pvBin[i]) && ((tr_goodPV) < pvBin[i + 1])) { + nPV = i; + break; + } + } + double den = 1.0 / ((*tr_TrkP)[itk]); + if ((ip >= 0) && (ie >= 0) && ((*tr_FE7x7P)[itk] > 0.02) && ((*tr_H3x3)[itk] > 0.1) && (nPV >= 0)) { + h_energy[nPV][ip][ie][0]->Fill(den * (*tr_FE7x7P)[itk], tr_eventWeight); + h_energy[nPV][ip][ie][1]->Fill(den * (*tr_H3x3)[itk], tr_eventWeight); + h_energy[nPV][ip][ie][2]->Fill(den * ((*tr_FE7x7P)[itk] + (*tr_H3x3)[itk]), tr_eventWeight); + h_energy[nPV][ip][ie][3]->Fill(den * (*tr_FE11x11P)[itk], tr_eventWeight); + h_energy[nPV][ip][ie][4]->Fill(den * (*tr_H5x5)[itk], tr_eventWeight); + h_energy[nPV][ip][ie][5]->Fill(den * ((*tr_FE11x11P)[itk] + (*tr_H5x5)[itk]), tr_eventWeight); + h_energy[nPVBin][ip][ie][0]->Fill(den * (*tr_FE7x7P)[itk], tr_eventWeight); + h_energy[nPVBin][ip][ie][1]->Fill(den * (*tr_H3x3)[itk], tr_eventWeight); + h_energy[nPVBin][ip][ie][2]->Fill(den * ((*tr_FE7x7P)[itk] + (*tr_H3x3)[itk]), tr_eventWeight); + h_energy[nPVBin][ip][ie][3]->Fill(den * (*tr_FE11x11P)[itk], tr_eventWeight); + h_energy[nPVBin][ip][ie][4]->Fill(den * (*tr_H5x5)[itk], tr_eventWeight); + h_energy[nPVBin][ip][ie][5]->Fill(den * ((*tr_FE11x11P)[itk] + (*tr_H5x5)[itk]), tr_eventWeight); + } + } //HCal Neutral Iso + } //neutral isolation + } //HCal Charged Iso + } //charged Iso } - }//end for loop + } //end for loop } f->Write(); f->Close(); } -void GetPUWeight(std::string mcFile, std::string dataFile, int type=0, - std::string dirName="StudyHLT") { - - std::string hName = (type == 0) ? "h_numberPV" : "h_goodPV"; - TFile *file1 = new TFile(mcFile.c_str()); - TDirectory *dir1 = (TDirectory*)file1->FindObjectAny(dirName.c_str()); - TH1D *histM = (TH1D*)dir1->Get(hName.c_str()); - TFile *file2 = new TFile(dataFile.c_str()); - TDirectory *dir2 = (TDirectory*)file2->FindObjectAny(dirName.c_str()); - TH1D *histD = (TH1D*)dir2->Get(hName.c_str()); - double scale = histM->Integral()/histD->Integral(); -//std::cout << "Scale " << scale << std::endl; +void GetPUWeight(std::string mcFile, std::string dataFile, int type = 0, std::string dirName = "StudyHLT") { + std::string hName = (type == 0) ? "h_numberPV" : "h_goodPV"; + TFile *file1 = new TFile(mcFile.c_str()); + TDirectory *dir1 = (TDirectory *)file1->FindObjectAny(dirName.c_str()); + TH1D *histM = (TH1D *)dir1->Get(hName.c_str()); + TFile *file2 = new TFile(dataFile.c_str()); + TDirectory *dir2 = (TDirectory *)file2->FindObjectAny(dirName.c_str()); + TH1D *histD = (TH1D *)dir2->Get(hName.c_str()); + double scale = histM->Integral() / histD->Integral(); std::vector weight; - for (int k=1; k<=histM->GetNbinsX(); ++k) { + for (int k = 1; k <= histM->GetNbinsX(); ++k) { double num = histD->GetBinContent(k); double den = histM->GetBinContent(k); - double rat = (den > 0) ? (scale*num/den) : 0; + double rat = (den > 0) ? (scale * num / den) : 0; weight.push_back(rat); -// std::cout << "Bin[" << k << "] " << num << ":" << den << ":" << rat << std::endl; } char buff[100]; - for (int k=0; kGetNbinsX(); k+=10) { - sprintf(buff,"%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,", - weight[k+0],weight[k+1],weight[k+2],weight[k+3],weight[k+4], - weight[k+5],weight[k+6],weight[k+7],weight[k+8],weight[k+9]); + for (int k = 0; k < histM->GetNbinsX(); k += 10) { + sprintf(buff, + "%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,", + weight[k + 0], + weight[k + 1], + weight[k + 2], + weight[k + 3], + weight[k + 4], + weight[k + 5], + weight[k + 6], + weight[k + 7], + weight[k + 8], + weight[k + 9]); std::cout << buff << std::endl; } } -TCanvas* PlotHist(std::string fileName, int type) { - std::string names[5] = {"h_goodPV", "h_numberPV", "h_maxNearP_Ecal", - "h_maxNearP_Hcal","h_p_HcalNeutIso"}; - std::string xtitl[5] = {"# Good Primary Vertex", "# Primary Vertex", - "Highest Track Momentum in ECAL Isolation Zone (GeV)", - "Highest Track Momentum in HCAL Isolation Zone (GeV)", - "Track Momentum (GeV)"}; - double xmin[5] = {0,0,-2,-2,0}; - double xmax[5] = {10,10,10,10,60}; +TCanvas *PlotHist(std::string fileName, int type) { + std::string names[5] = {"h_goodPV", "h_numberPV", "h_maxNearP_Ecal", "h_maxNearP_Hcal", "h_p_HcalNeutIso"}; + std::string xtitl[5] = {"# Good Primary Vertex", + "# Primary Vertex", + "Highest Track Momentum in ECAL Isolation Zone (GeV)", + "Highest Track Momentum in HCAL Isolation Zone (GeV)", + "Track Momentum (GeV)"}; + double xmin[5] = {0, 0, -2, -2, 0}; + double xmax[5] = {10, 10, 10, 10, 60}; - TCanvas* pad(0); - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); + TCanvas *pad(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); gStyle->SetOptTitle(0); - gStyle->SetOptStat(1110); gStyle->SetOptFit(0); - TFile* file = new TFile(fileName.c_str()); + gStyle->SetOptStat(1110); + gStyle->SetOptFit(0); + TFile *file = new TFile(fileName.c_str()); if (file) { - TH1D *hist = (TH1D*)file->FindObjectAny(names[type].c_str()); + TH1D *hist = (TH1D *)file->FindObjectAny(names[type].c_str()); if (hist) { char name[100]; - sprintf(name,"%s",names[type].c_str()); - pad = new TCanvas(name,name,700,500); - pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); + sprintf(name, "%s", names[type].c_str()); + pad = new TCanvas(name, name, 700, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); pad->SetLogy(); hist->SetLineColor(1); hist->SetMarkerColor(1); @@ -427,19 +452,23 @@ TCanvas* PlotHist(std::string fileName, int type) { hist->GetYaxis()->SetTitle("Events"); hist->GetYaxis()->SetLabelOffset(0.005); hist->GetYaxis()->SetTitleOffset(1.20); - hist->GetXaxis()->SetRangeUser(xmin[type],xmax[type]); + hist->GetXaxis()->SetRangeUser(xmin[type], xmax[type]); hist->Draw(); - pad->Modified(); pad->Update(); - TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats"); + pad->Modified(); + pad->Update(); + TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats"); std::cout << "Pad " << pad << " st " << st1 << std::endl; if (st1 != NULL) { - st1->SetFillColor(kWhite); - st1->SetLineColor(1); - st1->SetTextColor(1); - st1->SetY1NDC(0.78); st1->SetY2NDC(0.90); - st1->SetX1NDC(0.60); st1->SetX2NDC(0.90); + st1->SetFillColor(kWhite); + st1->SetLineColor(1); + st1->SetTextColor(1); + st1->SetY1NDC(0.78); + st1->SetY2NDC(0.90); + st1->SetX1NDC(0.60); + st1->SetX2NDC(0.90); } - pad->Modified(); pad->Update(); + pad->Modified(); + pad->Update(); } } return pad; diff --git a/Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.cc b/Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.cc index eb1d0e340e3ad..7ac884a2a6fdc 100644 --- a/Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.cc +++ b/Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.cc @@ -2764,7 +2764,7 @@ void IsolatedTracksNxN::printTrack(const reco::Track *pTrack) { std::ostringstream st1; st1 << "default "; for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::TRACK_HITS); i++) { - p.printHitPattern(reco::HitPattern::TRACK_HITS, i, std::cout); + p.printHitPattern(reco::HitPattern::TRACK_HITS, i, st1); } edm::LogVerbatim("IsoTrack") << st1.str(); std::ostringstream st2;