CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/DQMOffline/EGamma/plugins/PhotonAnalyzer.cc

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <iomanip>
00003 //
00004 
00005 #include "DQMOffline/EGamma/plugins/PhotonAnalyzer.h"
00006 
00007 
00019 using namespace std;
00020 
00021 
00022 PhotonAnalyzer::PhotonAnalyzer( const edm::ParameterSet& pset )
00023 {
00024 
00025     fName_                  = pset.getUntrackedParameter<string>("Name");
00026     verbosity_              = pset.getUntrackedParameter<int>("Verbosity");
00027 
00028     prescaleFactor_         = pset.getUntrackedParameter<int>("prescaleFactor",1);
00029 
00030     photonProducer_         = pset.getParameter<string>("phoProducer");
00031     photonCollection_       = pset.getParameter<string>("photonCollection");
00032 
00033     barrelRecHitProducer_   = pset.getParameter<string>("barrelRecHitProducer");
00034     barrelRecHitCollection_ = pset.getParameter<string>("barrelRecHitCollection");
00035 
00036     endcapRecHitProducer_   = pset.getParameter<string>("endcapRecHitProducer");
00037     endcapRecHitCollection_ = pset.getParameter<string>("endcapRecHitCollection");
00038 
00039     triggerEvent_           = pset.getParameter<edm::InputTag>("triggerEvent");
00040 
00041     minPhoEtCut_            = pset.getParameter<double>("minPhoEtCut");
00042     invMassEtCut_           = pset.getParameter<double>("invMassEtCut");
00043 
00044     cutStep_                = pset.getParameter<double>("cutStep");
00045     numberOfSteps_          = pset.getParameter<int>("numberOfSteps");
00046 
00047     useBinning_             = pset.getParameter<bool>("useBinning");
00048     useTriggerFiltering_    = pset.getParameter<bool>("useTriggerFiltering");
00049 
00050     minimalSetOfHistos_ = pset.getParameter<bool>("minimalSetOfHistos");
00051     excludeBkgHistos_       = pset.getParameter<bool>("excludeBkgHistos");
00052 
00053     standAlone_             = pset.getParameter<bool>("standAlone");
00054     outputFileName_         = pset.getParameter<string>("OutputFileName");
00055 
00056     isolationStrength_      = pset.getParameter<int>("isolationStrength");
00057 
00058     isHeavyIon_             = pset.getUntrackedParameter<bool>("isHeavyIon",false);
00059 
00060     parameters_ = pset;
00061 
00062     histo_index_photons_ = 0;
00063     histo_index_conversions_ = 0;
00064     histo_index_efficiency_ = 0;
00065     histo_index_invMass_ = 0;
00066 }
00067 
00068 
00069 
00070 PhotonAnalyzer::~PhotonAnalyzer() {}
00071 
00072 
00073 void PhotonAnalyzer::beginJob()
00074 {
00075 
00076   nEvt_=0;
00077   nEntry_=0;
00078 
00079   dbe_ = 0;
00080   dbe_ = edm::Service<DQMStore>().operator->();
00081 
00082 
00083 
00084   double eMin = parameters_.getParameter<double>("eMin");
00085   double eMax = parameters_.getParameter<double>("eMax");
00086   int    eBin = parameters_.getParameter<int>("eBin");
00087 
00088   double etMin = parameters_.getParameter<double>("etMin");
00089   double etMax = parameters_.getParameter<double>("etMax");
00090   int    etBin = parameters_.getParameter<int>("etBin");
00091 
00092   double sumMin = parameters_.getParameter<double>("sumMin");
00093   double sumMax = parameters_.getParameter<double>("sumMax");
00094   int    sumBin = parameters_.getParameter<int>("sumBin");
00095 
00096   double etaMin = parameters_.getParameter<double>("etaMin");
00097   double etaMax = parameters_.getParameter<double>("etaMax");
00098   int    etaBin = parameters_.getParameter<int>("etaBin");
00099 
00100   double phiMin = parameters_.getParameter<double>("phiMin");
00101   double phiMax = parameters_.getParameter<double>("phiMax");
00102   int    phiBin = parameters_.getParameter<int>("phiBin");
00103 
00104   double r9Min = parameters_.getParameter<double>("r9Min");
00105   double r9Max = parameters_.getParameter<double>("r9Max");
00106   int    r9Bin = parameters_.getParameter<int>("r9Bin");
00107 
00108   double hOverEMin = parameters_.getParameter<double>("hOverEMin");
00109   double hOverEMax = parameters_.getParameter<double>("hOverEMax");
00110   int    hOverEBin = parameters_.getParameter<int>("hOverEBin");
00111 
00112   double xMin = parameters_.getParameter<double>("xMin");
00113   double xMax = parameters_.getParameter<double>("xMax");
00114   int    xBin = parameters_.getParameter<int>("xBin");
00115 
00116   double yMin = parameters_.getParameter<double>("yMin");
00117   double yMax = parameters_.getParameter<double>("yMax");
00118   int    yBin = parameters_.getParameter<int>("yBin");
00119 
00120   double numberMin = parameters_.getParameter<double>("numberMin");
00121   double numberMax = parameters_.getParameter<double>("numberMax");
00122   int    numberBin = parameters_.getParameter<int>("numberBin");
00123 
00124   double zMin = parameters_.getParameter<double>("zMin");
00125   double zMax = parameters_.getParameter<double>("zMax");
00126   int    zBin = parameters_.getParameter<int>("zBin");
00127 
00128   double rMin = parameters_.getParameter<double>("rMin");
00129   double rMax = parameters_.getParameter<double>("rMax");
00130   int    rBin = parameters_.getParameter<int>("rBin");
00131 
00132   double dPhiTracksMin = parameters_.getParameter<double>("dPhiTracksMin");
00133   double dPhiTracksMax = parameters_.getParameter<double>("dPhiTracksMax");
00134   int    dPhiTracksBin = parameters_.getParameter<int>("dPhiTracksBin");
00135 
00136   double dEtaTracksMin = parameters_.getParameter<double>("dEtaTracksMin");
00137   double dEtaTracksMax = parameters_.getParameter<double>("dEtaTracksMax");
00138   int    dEtaTracksBin = parameters_.getParameter<int>("dEtaTracksBin");
00139 
00140   double sigmaIetaMin = parameters_.getParameter<double>("sigmaIetaMin");
00141   double sigmaIetaMax = parameters_.getParameter<double>("sigmaIetaMax");
00142   int    sigmaIetaBin = parameters_.getParameter<int>("sigmaIetaBin");
00143 
00144   double eOverPMin = parameters_.getParameter<double>("eOverPMin");
00145   double eOverPMax = parameters_.getParameter<double>("eOverPMax");
00146   int    eOverPBin = parameters_.getParameter<int>("eOverPBin");
00147 
00148   double chi2Min = parameters_.getParameter<double>("chi2Min");
00149   double chi2Max = parameters_.getParameter<double>("chi2Max");
00150   int    chi2Bin = parameters_.getParameter<int>("chi2Bin");
00151 
00152 
00153   int reducedEtBin = etBin/4;
00154   int reducedEtaBin = etaBin/4;
00155   int reducedSumBin = sumBin/4;
00156   int reducedR9Bin = r9Bin/4;
00157 
00158 
00159   parts_.push_back("AllEcal");
00160   parts_.push_back("Barrel");
00161   parts_.push_back("Endcaps");
00162 
00163   types_.push_back("All");
00164   types_.push_back("GoodCandidate");
00165   if (!excludeBkgHistos_) types_.push_back("Background");
00166 
00167 
00168 
00170 
00171   if (dbe_) {
00172 
00173     dbe_->setCurrentFolder("Egamma/PhotonAnalyzer");
00174 
00175     //int values stored in MEs to keep track of how many histograms are in each folder
00176     totalNumberOfHistos_efficiencyFolder =  dbe_->bookInt("numberOfHistogramsInEfficiencyFolder");
00177     totalNumberOfHistos_photonsFolder =     dbe_->bookInt("numberOfHistogramsInPhotonsFolder");
00178     totalNumberOfHistos_conversionsFolder = dbe_->bookInt("numberOfHistogramsInConversionsFolder");
00179     totalNumberOfHistos_invMassFolder =     dbe_->bookInt("numberOfHistogramsInInvMassFolder");
00180 
00181 
00182     //Efficiency histograms
00183 
00184     dbe_->setCurrentFolder("Egamma/PhotonAnalyzer/Efficiencies");
00185 
00186     //don't number these histograms with the "bookHisto" method, since they'll be erased in the offline client
00187     h_phoEta_Loose_ = dbe_->book1D("phoEtaLoose","Loose Photon #eta",etaBin,etaMin,etaMax);
00188     h_phoEta_Tight_ = dbe_->book1D("phoEtaTight","Tight Photon #eta",etaBin,etaMin,etaMax);
00189     h_phoEt_Loose_  = dbe_->book1D("phoEtLoose", "Loose Photon E_{T}",etBin,etMin,etMax);
00190     h_phoEt_Tight_  = dbe_->book1D("phoEtTight", "Tight Photon E_{T}",etBin,etMin,etMax);
00191 
00192 
00193     h_phoEta_preHLT_  = dbe_->book1D("phoEtaPreHLT", "Photon #eta: before HLT",etaBin,etaMin,etaMax);
00194     h_phoEta_postHLT_ = dbe_->book1D("phoEtaPostHLT","Photon #eta: after HLT",etaBin,etaMin,etaMax);
00195     h_phoEt_preHLT_   = dbe_->book1D("phoEtPreHLT",  "Photon E_{T}: before HLT",etBin,etMin,etMax);
00196     h_phoEt_postHLT_  = dbe_->book1D("phoEtPostHLT", "Photon E_{T}: after HLT",etBin,etMin,etMax);
00197 
00198     h_convEta_Loose_ = dbe_->book1D("convEtaLoose","Converted Loose Photon #eta",etaBin,etaMin,etaMax);
00199     h_convEta_Tight_ = dbe_->book1D("convEtaTight","Converted Tight Photon #eta",etaBin,etaMin,etaMax);
00200     h_convEt_Loose_  = dbe_->book1D("convEtLoose", "Converted Loose Photon E_{T}",etBin,etMin,etMax);
00201     h_convEt_Tight_  = dbe_->book1D("convEtTight", "Converted Tight Photon E_{T}",etBin,etMin,etMax);
00202 
00203     h_phoEta_Vertex_ = dbe_->book1D("phoEtaVertex","Converted Photons before valid vertex cut: #eta",etaBin,etaMin,etaMax);
00204 
00205 
00206     vector<MonitorElement*> temp1DVectorEta;
00207     vector<MonitorElement*> temp1DVectorPhi;
00208     vector<vector<MonitorElement*> > temp2DVectorPhi;
00209 
00210 
00211     for(int cut = 0; cut != numberOfSteps_; ++cut){ //looping over Et cut values
00212       for(uint type=0;type!=types_.size();++type){  //looping over isolation type
00213         currentFolder_.str("");
00214         currentFolder_ << "Egamma/PhotonAnalyzer/" << types_[type] << "Photons/Et above " << (cut+1)*cutStep_ << " GeV/Conversions";
00215         dbe_->setCurrentFolder(currentFolder_.str());
00216 
00217         temp1DVectorEta.push_back(dbe_->book1D("phoConvEtaForEfficiency","Converted Photon #eta;#eta",etaBin,etaMin,etaMax));
00218         for(uint part=0;part!=parts_.size();++part){
00219           temp1DVectorPhi.push_back(dbe_->book1D("phoConvPhiForEfficiency"+parts_[part],"Converted Photon #phi;#phi",phiBin,phiMin,phiMax));
00220         }
00221         temp2DVectorPhi.push_back(temp1DVectorPhi);
00222         temp1DVectorPhi.clear();
00223       }
00224       h_phoConvEtaForEfficiency_.push_back(temp1DVectorEta);
00225       temp1DVectorEta.clear();
00226       h_phoConvPhiForEfficiency_.push_back(temp2DVectorPhi);
00227       temp2DVectorPhi.clear();
00228     }
00229 
00230 
00231 
00232 
00233     //Invariant mass plots
00234 
00235     dbe_->setCurrentFolder("Egamma/PhotonAnalyzer/InvMass");
00236 
00237     h_invMassAllPhotons_    = bookHisto("invMassAllIsolatedPhotons","Two photon invariant mass: All isolated photons;M (GeV)",etBin,etMin,etMax);
00238     h_invMassPhotonsEBarrel_    = bookHisto("invMassIsoPhotonsEBarrel", "Two photon invariant mass: isolated photons in barrel; M (GeV)",etBin,etMin,etMax);
00239     h_invMassPhotonsEEndcap_    = bookHisto("invMassIsoPhotonsEEndcap", "Two photon invariant mass: isolated photons in endcap; M (GeV)",etBin,etMin,etMax);
00240     
00241     h_invMassZeroWithTracks_= bookHisto("invMassZeroWithTracks",    "Two photon invariant mass: Neither has tracks;M (GeV)",  etBin,etMin,etMax);
00242     h_invMassOneWithTracks_ = bookHisto("invMassOneWithTracks",     "Two photon invariant mass: Only one has tracks;M (GeV)", etBin,etMin,etMax);
00243     h_invMassTwoWithTracks_ = bookHisto("invMassTwoWithTracks",     "Two photon invariant mass: Both have tracks;M (GeV)",    etBin,etMin,etMax);
00244     
00245 
00247 
00248     //ENERGY VARIABLES
00249 
00250     book3DHistoVector(h_phoE_, "1D","phoE","Energy;E (GeV)",eBin,eMin,eMax);
00251     book3DHistoVector(h_phoEt_, "1D","phoEt","E_{T};E_{T} (GeV)", etBin,etMin,etMax);
00252 
00253     //NUMBER OF PHOTONS
00254 
00255     book3DHistoVector(h_nPho_, "1D","nPho","Number of Photons per Event;# #gamma",numberBin,numberMin,numberMax);
00256 
00257     //GEOMETRICAL VARIABLES
00258 
00259     //photon eta/phi
00260     book2DHistoVector(h_phoEta_, "1D","phoEta","#eta;#eta",etaBin,etaMin,etaMax) ;
00261     book3DHistoVector(h_phoPhi_, "1D","phoPhi","#phi;#phi",phiBin,phiMin,phiMax) ;
00262 
00263     //supercluster eta/phi
00264     book2DHistoVector(h_scEta_, "1D","scEta","SuperCluster #eta;#eta",etaBin,etaMin,etaMax) ;
00265     book3DHistoVector(h_scPhi_, "1D","scPhi","SuperCluster #phi;#phi",phiBin,phiMin,phiMax) ;
00266 
00267     //SHOWER SHAPE VARIABLES
00268 
00269     //r9
00270     book3DHistoVector(h_r9_, "1D","r9","R9;R9",r9Bin,r9Min, r9Max);
00271     book2DHistoVector(h_r9VsEt_, "2D","r9VsEt2D","R9 vs E_{T};E_{T} (GeV);R9",reducedEtBin,etMin,etMax,reducedR9Bin,r9Min,r9Max);
00272     book2DHistoVector(p_r9VsEt_, "Profile","r9VsEt","Avg R9 vs E_{T};E_{T} (GeV);R9",etBin,etMin,etMax,r9Bin,r9Min,r9Max);
00273     book2DHistoVector(h_r9VsEta_, "2D","r9VsEta2D","R9 vs #eta;#eta;R9",reducedEtaBin,etaMin,etaMax,reducedR9Bin,r9Min,r9Max);
00274     book2DHistoVector(p_r9VsEta_, "Profile","r9VsEta","Avg R9 vs #eta;#eta;R9",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
00275 
00276     //sigma ieta ieta
00277     book3DHistoVector(h_phoSigmaIetaIeta_,   "1D","phoSigmaIetaIeta","#sigma_{i#etai#eta};#sigma_{i#etai#eta}",sigmaIetaBin,sigmaIetaMin,sigmaIetaMax);
00278     book2DHistoVector(h_sigmaIetaIetaVsEta_, "2D","sigmaIetaIetaVsEta2D","#sigma_{i#etai#eta} vs #eta;#eta;#sigma_{i#etai#eta}",reducedEtaBin,etaMin,etaMax,sigmaIetaBin,sigmaIetaMin,sigmaIetaMax);
00279     book2DHistoVector(p_sigmaIetaIetaVsEta_, "Profile","sigmaIetaIetaVsEta","Avg #sigma_{i#etai#eta} vs #eta;#eta;#sigma_{i#etai#eta}",etaBin,etaMin,etaMax,sigmaIetaBin,sigmaIetaMin,sigmaIetaMax);
00280 
00281     //e1x5
00282     book2DHistoVector(h_e1x5VsEt_,  "2D","e1x5VsEt2D","E1x5 vs E_{T};E_{T} (GeV);E1X5 (GeV)",reducedEtBin,etMin,etMax,reducedEtBin,etMin,etMax);
00283     book2DHistoVector(p_e1x5VsEt_,  "Profile","e1x5VsEt","Avg E1x5 vs E_{T};E_{T} (GeV);E1X5 (GeV)",etBin,etMin,etMax,etBin,etMin,etMax);
00284     book2DHistoVector(h_e1x5VsEta_, "2D","e1x5VsEta2D","E1x5 vs #eta;#eta;E1X5 (GeV)",reducedEtaBin,etaMin,etaMax,reducedEtBin,etMin,etMax);
00285     book2DHistoVector(p_e1x5VsEta_, "Profile","e1x5VsEta","Avg E1x5 vs #eta;#eta;E1X5 (GeV)",etaBin,etaMin,etaMax,etBin,etMin,etMax);
00286 
00287     //e2x5
00288     book2DHistoVector(h_e2x5VsEt_,  "2D","e2x5VsEt2D","E2x5 vs E_{T};E_{T} (GeV);E2X5 (GeV)",reducedEtBin,etMin,etMax,reducedEtBin,etMin,etMax);
00289     book2DHistoVector(p_e2x5VsEt_,  "Profile","e2x5VsEt","Avg E2x5 vs E_{T};E_{T} (GeV);E2X5 (GeV)",etBin,etMin,etMax,etBin,etMin,etMax);
00290     book2DHistoVector(h_e2x5VsEta_, "2D","e2x5VsEta2D","E2x5 vs #eta;#eta;E2X5 (GeV)",reducedEtaBin,etaMin,etaMax,reducedEtBin,etMin,etMax);
00291     book2DHistoVector(p_e2x5VsEta_, "Profile","e2x5VsEta","Avg E2x5 vs #eta;#eta;E2X5 (GeV)",etaBin,etaMin,etaMax,etBin,etMin,etMax);
00292 
00293     //r1x5
00294     book2DHistoVector(h_r1x5VsEt_,  "2D","r1x5VsEt2D","R1x5 vs E_{T};E_{T} (GeV);R1X5",reducedEtBin,etMin,etMax,reducedR9Bin,r9Min,r9Max);
00295     book2DHistoVector(p_r1x5VsEt_,  "Profile","r1x5VsEt","Avg R1x5 vs E_{T};E_{T} (GeV);R1X5",etBin,etMin,etMax,r9Bin,r9Min,r9Max);
00296     book2DHistoVector(h_r1x5VsEta_, "2D","r1x5VsEta2D","R1x5 vs #eta;#eta;R1X5",reducedEtaBin,etaMin,etaMax,reducedR9Bin,r9Min,r9Max);
00297     book2DHistoVector(p_r1x5VsEta_, "Profile","r1x5VsEta","Avg R1x5 vs #eta;#eta;R1X5",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
00298 
00299     //r2x5
00300     book2DHistoVector(    h_r2x5VsEt_  ,"2D","r2x5VsEt2D","R2x5 vs E_{T};E_{T} (GeV);R2X5",reducedEtBin,etMin,etMax,reducedR9Bin,r9Min,r9Max);
00301     book2DHistoVector(    p_r2x5VsEt_  ,"Profile","r2x5VsEt","Avg R2x5 vs E_{T};E_{T} (GeV);R2X5",etBin,etMin,etMax,r9Bin,r9Min,r9Max);
00302     book2DHistoVector(    h_r2x5VsEta_ ,"2D","r2x5VsEta2D","R2x5 vs #eta;#eta;R2X5",reducedEtaBin,etaMin,etaMax,reducedR9Bin,r9Min,r9Max);
00303     book2DHistoVector(    p_r2x5VsEta_ ,"Profile","r2x5VsEta","Avg R2x5 vs #eta;#eta;R2X5",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
00304 
00305     //maxEXtalOver3x3
00306     book2DHistoVector(    h_maxEXtalOver3x3VsEt_  ,"2D","maxEXtalOver3x3VsEt2D","(Max Xtal E)/E3x3 vs E_{T};E_{T} (GeV);(Max Xtal E)/E3x3",reducedEtBin,etMin,etMax,r9Bin,r9Min,r9Max);
00307     book2DHistoVector(    p_maxEXtalOver3x3VsEt_  ,"Profile","maxEXtalOver3x3VsEt","Avg (Max Xtal E)/E3x3 vs E_{T};E_{T} (GeV);(Max Xtal E)/E3x3",etBin,etMin,etMax,r9Bin,r9Min,r9Max);
00308     book2DHistoVector(    h_maxEXtalOver3x3VsEta_ ,"2D","maxEXtalOver3x3VsEta2D","(Max Xtal E)/E3x3 vs #eta;#eta;(Max Xtal E)/E3x3",reducedEtaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
00309     book2DHistoVector(    p_maxEXtalOver3x3VsEta_ ,"Profile","maxEXtalOver3x3VsEta","Avg (Max Xtal E)/E3x3 vs #eta;#eta;(Max Xtal E)/E3x3",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
00310 
00311 
00312     //TRACK ISOLATION VARIABLES
00313 
00314     //nTrackIsolSolid
00315     book2DHistoVector(    h_nTrackIsolSolid_       ,"1D","nIsoTracksSolid","Number Of Tracks in the Solid Iso Cone;# tracks",numberBin,numberMin,numberMax);
00316     book2DHistoVector(    h_nTrackIsolSolidVsEt_   ,"2D","nIsoTracksSolidVsEt2D","Number Of Tracks in the Solid Iso Cone vs E_{T};E_{T};# tracks",reducedEtBin,etMin, etMax,numberBin,numberMin,numberMax);
00317     book2DHistoVector(    p_nTrackIsolSolidVsEt_   ,"Profile","nIsoTracksSolidVsEt","Avg Number Of Tracks in the Solid Iso Cone vs E_{T};E_{T};# tracks",etBin,etMin,etMax,numberBin,numberMin,numberMax);
00318     book2DHistoVector(    h_nTrackIsolSolidVsEta_  ,"2D","nIsoTracksSolidVsEta2D","Number Of Tracks in the Solid Iso Cone vs #eta;#eta;# tracks",reducedEtaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
00319     book2DHistoVector(    p_nTrackIsolSolidVsEta_  ,"Profile","nIsoTracksSolidVsEta","Avg Number Of Tracks in the Solid Iso Cone vs #eta;#eta;# tracks",etaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
00320 
00321     //nTrackIsolHollow
00322     book2DHistoVector(    h_nTrackIsolHollow_      ,"1D","nIsoTracksHollow","Number Of Tracks in the Hollow Iso Cone;# tracks",numberBin,numberMin,numberMax);
00323     book2DHistoVector(    h_nTrackIsolHollowVsEt_  ,"2D","nIsoTracksHollowVsEt2D","Number Of Tracks in the Hollow Iso Cone vs E_{T};E_{T};# tracks",reducedEtBin,etMin, etMax,numberBin,numberMin,numberMax);
00324     book2DHistoVector(    p_nTrackIsolHollowVsEt_  ,"Profile","nIsoTracksHollowVsEt","Avg Number Of Tracks in the Hollow Iso Cone vs E_{T};E_{T};# tracks",etBin,etMin,etMax,numberBin,numberMin,numberMax);
00325     book2DHistoVector(    h_nTrackIsolHollowVsEta_ ,"2D","nIsoTracksHollowVsEta2D","Number Of Tracks in the Hollow Iso Cone vs #eta;#eta;# tracks",reducedEtaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
00326     book2DHistoVector(    p_nTrackIsolHollowVsEta_ ,"Profile","nIsoTracksHollowVsEta","Avg Number Of Tracks in the Hollow Iso Cone vs #eta;#eta;# tracks",etaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
00327 
00328     //trackPtSumSolid
00329     book2DHistoVector(    h_trackPtSumSolid_       ,"1D","isoPtSumSolid","Track P_{T} Sum in the Solid Iso Cone;P_{T} (GeV)",sumBin,sumMin,sumMax);
00330     book2DHistoVector(    h_trackPtSumSolidVsEt_   ,"2D","isoPtSumSolidVsEt2D","Track P_{T} Sum in the Solid Iso Cone;E_{T} (GeV);P_{T} (GeV)",reducedEtBin,etMin, etMax,reducedSumBin,sumMin,sumMax);
00331     book2DHistoVector(    p_trackPtSumSolidVsEt_   ,"Profile","isoPtSumSolidVsEt","Avg Track P_{T} Sum in the Solid Iso Cone vs E_{T};E_{T} (GeV);P_{T} (GeV)",etBin,etMin,etMax,sumBin,sumMin,sumMax);
00332     book2DHistoVector(    h_trackPtSumSolidVsEta_  ,"2D","isoPtSumSolidVsEta2D","Track P_{T} Sum in the Solid Iso Cone;#eta;P_{T} (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
00333     book2DHistoVector(    p_trackPtSumSolidVsEta_  ,"Profile","isoPtSumSolidVsEta","Avg Track P_{T} Sum in the Solid Iso Cone vs #eta;#eta;P_{T} (GeV)",etaBin,etaMin, etaMax,sumBin,sumMin,sumMax);
00334 
00335     //trackPtSumHollow
00336     book2DHistoVector(    h_trackPtSumHollow_      ,"1D","isoPtSumHollow","Track P_{T} Sum in the Hollow Iso Cone;P_{T} (GeV)",sumBin,sumMin,sumMax);
00337     book2DHistoVector(    h_trackPtSumHollowVsEt_  ,"2D","isoPtSumHollowVsEt2D","Track P_{T} Sum in the Hollow Iso Cone;E_{T} (GeV);P_{T} (GeV)",reducedEtBin,etMin, etMax,reducedSumBin,sumMin,sumMax);
00338     book2DHistoVector(    p_trackPtSumHollowVsEt_  ,"Profile","isoPtSumHollowVsEt","Avg Track P_{T} Sum in the Hollow Iso Cone vs E_{T};E_{T} (GeV);P_{T} (GeV)",etBin,etMin,etMax,sumBin,sumMin,sumMax);
00339     book2DHistoVector(    h_trackPtSumHollowVsEta_ ,"2D","isoPtSumHollowVsEta2D","Track P_{T} Sum in the Hollow Iso Cone;#eta;P_{T} (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
00340     book2DHistoVector(    p_trackPtSumHollowVsEta_ ,"Profile","isoPtSumHollowVsEta","Avg Track P_{T} Sum in the Hollow Iso Cone vs #eta;#eta;P_{T} (GeV)",etaBin,etaMin, etaMax,sumBin,sumMin,sumMax);
00341 
00342 
00343     //CALORIMETER ISOLATION VARIABLES
00344 
00345     //ecal sum
00346     book2DHistoVector(    h_ecalSum_      ,"1D","ecalSum","Ecal Sum in the Iso Cone;E (GeV)",sumBin,sumMin,sumMax);
00347     book2DHistoVector(    h_ecalSumEBarrel_,"1D","ecalSumEBarrel","Ecal Sum in the IsoCone for Barrel;E (GeV)",sumBin,sumMin,sumMax);
00348     book2DHistoVector(    h_ecalSumEEndcap_,"1D","ecalSumEEndcap","Ecal Sum in the IsoCone for Endcap;E (GeV)",sumBin,sumMin,sumMax);
00349     book2DHistoVector(    h_ecalSumVsEt_  ,"2D","ecalSumVsEt2D","Ecal Sum in the Iso Cone;E_{T} (GeV);E (GeV)",reducedEtBin,etMin, etMax,reducedSumBin,sumMin,sumMax);
00350     book3DHistoVector(    p_ecalSumVsEt_  ,"Profile","ecalSumVsEt","Avg Ecal Sum in the Iso Cone vs E_{T};E_{T} (GeV);E (GeV)",etBin,etMin, etMax,sumBin,sumMin,sumMax);
00351     book2DHistoVector(    h_ecalSumVsEta_ ,"2D","ecalSumVsEta2D","Ecal Sum in the Iso Cone;#eta;E (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
00352     book2DHistoVector(    p_ecalSumVsEta_ ,"Profile","ecalSumVsEta","Avg Ecal Sum in the Iso Cone vs #eta;#eta;E (GeV)",etaBin,etaMin, etaMax,sumBin,sumMin,sumMax);
00353 
00354     //hcal sum
00355     book2DHistoVector(    h_hcalSum_      ,"1D","hcalSum","Hcal Sum in the Iso Cone;E (GeV)",sumBin,sumMin,sumMax);
00356     book2DHistoVector(    h_hcalSumEBarrel_,"1D","hcalSumEBarrel","Hcal Sum in the IsoCone for Barrel;E (GeV)",sumBin,sumMin,sumMax);
00357     book2DHistoVector(    h_hcalSumEEndcap_,"1D","hcalSumEEndcap","Hcal Sum in the IsoCone for Endcap;E (GeV)",sumBin,sumMin,sumMax);
00358     book2DHistoVector(    h_hcalSumVsEt_  ,"2D","hcalSumVsEt2D","Hcal Sum in the Iso Cone;E_{T} (GeV);E (GeV)",reducedEtBin,etMin, etMax,reducedSumBin,sumMin,sumMax);
00359     book3DHistoVector(    p_hcalSumVsEt_  ,"Profile","hcalSumVsEt","Avg Hcal Sum in the Iso Cone vs E_{T};E_{T} (GeV);E (GeV)",etBin,etMin, etMax,sumBin,sumMin,sumMax);
00360     book2DHistoVector(    h_hcalSumVsEta_ ,"2D","hcalSumVsEta2D","Hcal Sum in the Iso Cone;#eta;E (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
00361     book2DHistoVector(    p_hcalSumVsEta_ ,"Profile","hcalSumVsEta","Avg Hcal Sum in the Iso Cone vs #eta;#eta;E (GeV)",etaBin,etaMin, etaMax,sumBin,sumMin,sumMax);
00362 
00363     //h over e
00364     book3DHistoVector(    h_hOverE_       ,"1D","hOverE","H/E;H/E",hOverEBin,hOverEMin,hOverEMax);
00365     book2DHistoVector(    p_hOverEVsEt_   ,"Profile","hOverEVsEt","Avg H/E vs Et;E_{T} (GeV);H/E",etBin,etMin,etMax,hOverEBin,hOverEMin,hOverEMax);
00366     book2DHistoVector(    p_hOverEVsEta_  ,"Profile","hOverEVsEta","Avg H/E vs #eta;#eta;H/E",etaBin,etaMin,etaMax,hOverEBin,hOverEMin,hOverEMax);
00367     book3DHistoVector(    h_h1OverE_      ,"1D","h1OverE","H/E for Depth 1;H/E",hOverEBin,hOverEMin,hOverEMax);
00368     book3DHistoVector(    h_h2OverE_      ,"1D","h2OverE","H/E for Depth 2;H/E",hOverEBin,hOverEMin,hOverEMax);
00369 
00370 
00371     //OTHER VARIABLES
00372 
00373     //bad channel histograms
00374     book2DHistoVector(    h_phoEt_BadChannels_  ,"1D","phoEtBadChannels", "Fraction Containing Bad Channels: E_{T};E_{T} (GeV)",etBin,etMin,etMax);
00375     book2DHistoVector(    h_phoEta_BadChannels_ ,"1D","phoEtaBadChannels","Fraction Containing Bad Channels: #eta;#eta",etaBin,etaMin,etaMax);
00376     book2DHistoVector(    h_phoPhi_BadChannels_ ,"1D","phoPhiBadChannels","Fraction Containing Bad Channels: #phi;#phi",phiBin,phiMin,phiMax);
00377 
00378 
00380 
00381     dbe_->setCurrentFolder("Egamma/PhotonAnalyzer/AllPhotons/Et Above 0 GeV/Conversions");
00382 
00383     //ENERGY VARIABLES
00384 
00385     book3DHistoVector(    h_phoConvE_  ,"1D","phoConvE","E;E (GeV)",eBin,eMin,eMax);
00386     book3DHistoVector(    h_phoConvEt_ ,"1D","phoConvEt","E_{T};E_{T} (GeV)",etBin,etMin,etMax);
00387 
00388     //GEOMETRICAL VARIABLES
00389 
00390     book2DHistoVector(    h_phoConvEta_ ,"1D","phoConvEta","#eta;#eta",etaBin,etaMin,etaMax);
00391     book3DHistoVector(    h_phoConvPhi_ ,"1D","phoConvPhi","#phi;#phi",phiBin,phiMin,phiMax);
00392 
00393     //NUMBER OF PHOTONS
00394 
00395     book3DHistoVector(    h_nConv_ ,"1D","nConv","Number Of Conversions per Event ;# conversions",numberBin,numberMin,numberMax);
00396 
00397     //SHOWER SHAPE VARIABLES
00398 
00399     book3DHistoVector(    h_phoConvR9_ ,"1D","phoConvR9","R9;R9",r9Bin,r9Min,r9Max);
00400 
00401     //TRACK RELATED VARIABLES
00402 
00403     book3DHistoVector(    h_eOverPTracks_ ,"1D","eOverPTracks","E/P;E/P",eOverPBin,eOverPMin,eOverPMax);
00404     book3DHistoVector(    h_pOverETracks_ ,"1D","pOverETracks","P/E;P/E",eOverPBin,eOverPMin,eOverPMax);
00405 
00406     book3DHistoVector(    h_dPhiTracksAtVtx_  ,"1D","dPhiTracksAtVtx", "#Delta#phi of Tracks at Vertex;#Delta#phi",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
00407     book3DHistoVector(    h_dPhiTracksAtEcal_ ,"1D","dPhiTracksAtEcal", "Abs(#Delta#phi) of Tracks at Ecal;#Delta#phi",dPhiTracksBin,0.,dPhiTracksMax);
00408     book3DHistoVector(    h_dEtaTracksAtEcal_ ,"1D","dEtaTracksAtEcal", "#Delta#eta of Tracks at Ecal;#Delta#eta",dEtaTracksBin,dEtaTracksMin,dEtaTracksMax);
00409 
00410     book3DHistoVector(    h_dCotTracks_      ,"1D","dCotTracks","#Deltacot(#theta) of Tracks;#Deltacot(#theta)",dEtaTracksBin,dEtaTracksMin,dEtaTracksMax);
00411     book2DHistoVector(    p_dCotTracksVsEta_ ,"Profile","dCotTracksVsEta","Avg #Deltacot(#theta) of Tracks vs #eta;#eta;#Deltacot(#theta)",etaBin,etaMin,etaMax,dEtaTracksBin,dEtaTracksMin,dEtaTracksMax);
00412 
00413     book2DHistoVector(    p_nHitsVsEta_ ,"Profile","nHitsVsEta","Avg Number of Hits per Track vs #eta;#eta;# hits",etaBin,etaMin,etaMax,etaBin,0,16);
00414 
00415     book2DHistoVector(    h_tkChi2_      ,"1D","tkChi2","#chi^{2} of Track Fitting;#chi^{2}",chi2Bin,chi2Min,chi2Max);
00416     book2DHistoVector(    p_tkChi2VsEta_ ,"Profile","tkChi2VsEta","Avg #chi^{2} of Track Fitting vs #eta;#eta;#chi^{2}",etaBin,etaMin,etaMax,chi2Bin,chi2Min,chi2Max);
00417 
00418     //VERTEX RELATED VARIABLES
00419 
00420     book2DHistoVector(    h_convVtxRvsZ_ ,"2D","convVtxRvsZ","Vertex Position;Z (cm);R (cm)",500,zMin,zMax,rBin,rMin,rMax);
00421     book2DHistoVector(    h_convVtxZEndcap_    ,"1D","convVtxZEndcap",   "Vertex Position: #eta > 1.5;Z (cm)",zBin,zMin,zMax);
00422     book2DHistoVector(    h_convVtxZ_    ,"1D","convVtxZ",   "Vertex Position;Z (cm)",zBin,zMin,zMax);
00423     book2DHistoVector(    h_convVtxR_    ,"1D","convVtxR",   "Vertex Position: #eta < 1;R (cm)",rBin,rMin,rMax);
00424     book2DHistoVector(    h_convVtxYvsX_ ,"2D","convVtxYvsX","Vertex Position: #eta < 1;X (cm);Y (cm)",xBin,xMin,xMax,yBin,yMin,yMax);
00425 
00426 
00427 
00428     book2DHistoVector(    h_vertexChi2Prob_ ,"1D","vertexChi2Prob","#chi^{2} Probability of Vertex Fitting;#chi^{2}",100,0.,1.0);
00429 
00430 
00431   }//end if(dbe_)
00432 
00433 
00434 }//end BeginJob
00435 
00436 
00437 
00438 void PhotonAnalyzer::analyze( const edm::Event& e, const edm::EventSetup& esup )
00439 {
00440   using namespace edm;
00441 
00442   if (nEvt_% prescaleFactor_ ) return;
00443   nEvt_++;
00444   LogInfo("PhotonAnalyzer") << "PhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
00445 
00446   // Get the trigger results
00447   bool validTriggerEvent=true;
00448   edm::Handle<trigger::TriggerEvent> triggerEventHandle;
00449   trigger::TriggerEvent triggerEvent;
00450   e.getByLabel(triggerEvent_,triggerEventHandle);
00451   if(!triggerEventHandle.isValid()) {
00452     edm::LogInfo("PhotonAnalyzer") << "Error! Can't get the product "<< triggerEvent_.label() << endl;
00453     validTriggerEvent=false;
00454   }
00455   if(validTriggerEvent) triggerEvent = *(triggerEventHandle.product());
00456 
00457   // Get the reconstructed photons
00458   bool validPhotons=true;
00459   Handle<reco::PhotonCollection> photonHandle;
00460   reco::PhotonCollection photonCollection;
00461   e.getByLabel(photonProducer_, photonCollection_ , photonHandle);
00462   if ( !photonHandle.isValid()) {
00463     edm::LogInfo("PhotonAnalyzer") << "Error! Can't get the product "<< photonCollection_ << endl;
00464     validPhotons=false;
00465   }
00466   if(validPhotons) photonCollection = *(photonHandle.product());
00467 
00468   // Get the PhotonId objects
00469   bool validloosePhotonID=true;
00470   Handle<edm::ValueMap<bool> > loosePhotonFlag;
00471   edm::ValueMap<bool> loosePhotonID;
00472   e.getByLabel("PhotonIDProd", "PhotonCutBasedIDLoose", loosePhotonFlag);
00473   if ( !loosePhotonFlag.isValid()) {
00474     edm::LogInfo("PhotonAnalyzer") << "Error! Can't get the product "<< "PhotonCutBasedIDLoose" << endl;
00475     validloosePhotonID=false;
00476   }
00477   if (validloosePhotonID) loosePhotonID = *(loosePhotonFlag.product());
00478 
00479   bool validtightPhotonID=true;
00480   Handle<edm::ValueMap<bool> > tightPhotonFlag;
00481   edm::ValueMap<bool> tightPhotonID;
00482   e.getByLabel("PhotonIDProd", "PhotonCutBasedIDTight", tightPhotonFlag);
00483   if ( !tightPhotonFlag.isValid()) {
00484     edm::LogInfo("PhotonAnalyzer") << "Error! Can't get the product "<< "PhotonCutBasedIDTight" << endl;
00485     validtightPhotonID=false;
00486   }
00487   if (validtightPhotonID) tightPhotonID = *(tightPhotonFlag.product());
00488 
00489 
00490 
00491   // Create array to hold #photons/event information
00492   int nPho[100][3][3];
00493 
00494   for (int cut=0; cut!=100; ++cut){
00495     for (unsigned int type=0; type!=types_.size(); ++type){
00496       for (unsigned int part=0; part!=parts_.size(); ++part){
00497         nPho[cut][type][part] = 0;
00498       }
00499     }
00500   }
00501   // Create array to hold #conversions/event information
00502   int nConv[100][3][3];
00503 
00504   for (int cut=0; cut!=100; ++cut){
00505     for (unsigned int type=0; type!=types_.size(); ++type){
00506       for (unsigned int part=0; part!=parts_.size(); ++part){
00507         nConv[cut][type][part] = 0;
00508       }
00509     }
00510   }
00511 
00512 
00513 
00514   //Prepare list of photon-related HLT filter names
00515 
00516   vector<int> Keys;
00517 
00518   for(uint filterIndex=0;filterIndex<triggerEvent.sizeFilters();++filterIndex){  //loop over all trigger filters in event (i.e. filters passed)
00519 
00520     string label = triggerEvent.filterTag(filterIndex).label();
00521 
00522     if(label.find( "Photon" ) != string::npos ) {  //get photon-related filters
00523 
00524       for(uint filterKeyIndex=0;filterKeyIndex<triggerEvent.filterKeys(filterIndex).size();++filterKeyIndex){  //loop over keys to objects passing this filter
00525         Keys.push_back(triggerEvent.filterKeys(filterIndex)[filterKeyIndex]);  //add keys to a vector for later reference
00526       }
00527 
00528     }
00529 
00530   }
00531 
00532   // sort Keys vector in ascending order
00533   // and erases duplicate entries from the vector
00534   sort(Keys.begin(),Keys.end());
00535   for ( uint i=0 ; i<Keys.size() ; )
00536    {
00537     if (i!=(Keys.size()-1))
00538      {
00539       if (Keys[i]==Keys[i+1]) Keys.erase(Keys.begin()+i+1) ;
00540       else ++i ;
00541      }
00542     else ++i ;
00543    }
00544 
00545   //We now have a vector of unique keys to TriggerObjects passing a photon-related filter
00546 
00547   int photonCounter = 0;
00548 
00550 
00551   for( reco::PhotonCollection::const_iterator  iPho = photonCollection.begin(); iPho != photonCollection.end(); iPho++) {
00552 
00553 
00554     //for HLT efficiency plots
00555 
00556     h_phoEta_preHLT_->Fill((*iPho).eta());
00557     h_phoEt_preHLT_->Fill( (*iPho).et());
00558 
00559 
00560     double deltaR=1000.;
00561     double deltaRMin=1000.;
00562     double deltaRMax=0.05;//sets deltaR threshold for matching photons to trigger objects
00563 
00564 
00565     for (vector<int>::const_iterator objectKey=Keys.begin();objectKey!=Keys.end();objectKey++){  //loop over keys to objects that fired photon triggers
00566 
00567       deltaR = reco::deltaR(triggerEvent.getObjects()[(*objectKey)].eta(),triggerEvent.getObjects()[(*objectKey)].phi(),(*iPho).superCluster()->eta(),(*iPho).superCluster()->phi());
00568       if(deltaR < deltaRMin) deltaRMin = deltaR;
00569 
00570     }
00571 
00572     if(deltaRMin > deltaRMax) {  //photon fails delta R cut
00573       if(useTriggerFiltering_) continue;  //throw away photons that haven't passed any photon filters
00574     }
00575 
00576     if(deltaRMin <= deltaRMax) { //photon passes delta R cut
00577       h_phoEta_postHLT_->Fill((*iPho).eta() );
00578       h_phoEt_postHLT_->Fill( (*iPho).et() );
00579     }
00580 
00581     if ((*iPho).et()  < minPhoEtCut_) continue;
00582 
00583     nEntry_++;
00584 
00585     edm::Ref<reco::PhotonCollection> photonref(photonHandle, photonCounter);
00586     photonCounter++;
00587 
00588     bool isLoosePhoton(false), isTightPhoton(false);
00589     if ( !isHeavyIon_ ) {
00590       isLoosePhoton = (loosePhotonID)[photonref];
00591       isTightPhoton = (tightPhotonID)[photonref];
00592     }
00593 
00594 
00595     //find out which part of the Ecal contains the photon
00596 
00597     bool  phoIsInBarrel=false;
00598     bool  phoIsInEndcap=false;
00599     float etaPho = (*iPho).superCluster()->eta();
00600     if ( fabs(etaPho) <  1.479 )
00601       phoIsInBarrel=true;
00602     else {
00603       phoIsInEndcap=true;
00604     }
00605 
00606     int part = 0;
00607     if ( phoIsInBarrel ) part = 1;
00608     if ( phoIsInEndcap ) part = 2;
00609 
00611     bool isIsolated=false;
00612     if ( isolationStrength_ == 0)  isIsolated = isLoosePhoton;
00613     if ( isolationStrength_ == 1)  isIsolated = isTightPhoton;
00614 
00615     int type=0;
00616     if ( isIsolated ) type=1;
00617     if ( !excludeBkgHistos_ && !isIsolated ) type=2;
00618 
00619 
00620     //get rechit collection containing this photon
00621     bool validEcalRecHits=true;
00622     edm::Handle<EcalRecHitCollection>   ecalRecHitHandle;
00623     EcalRecHitCollection ecalRecHitCollection;
00624     if ( phoIsInBarrel ) {
00625       // Get handle to barrel rec hits
00626       e.getByLabel(barrelRecHitProducer_, barrelRecHitCollection_, ecalRecHitHandle);
00627       if (!ecalRecHitHandle.isValid()) {
00628         edm::LogError("PhotonAnalyzer") << "Error! Can't get the product "<<barrelRecHitProducer_;
00629         validEcalRecHits=false;
00630       }
00631     }
00632     else if ( phoIsInEndcap ) {
00633       // Get handle to endcap rec hits
00634       e.getByLabel(endcapRecHitProducer_, endcapRecHitCollection_, ecalRecHitHandle);
00635       if (!ecalRecHitHandle.isValid()) {
00636         edm::LogError("PhotonAnalyzer") << "Error! Can't get the product "<<endcapRecHitProducer_;
00637         validEcalRecHits=false;
00638       }
00639     }
00640     if (validEcalRecHits) ecalRecHitCollection = *(ecalRecHitHandle.product());
00641 
00642 
00643     //if ((*iPho).isEBEEGap()) continue;  //cut out gap photons
00644 
00645 
00646     //filling histograms to make isolation efficiencies
00647     if(isLoosePhoton){
00648       h_phoEta_Loose_->Fill((*iPho).eta());
00649       h_phoEt_Loose_->Fill( (*iPho).et() );
00650     }
00651     if(isTightPhoton){
00652       h_phoEta_Tight_->Fill((*iPho).eta());
00653       h_phoEt_Tight_->Fill( (*iPho).et() );
00654     }
00655 
00656 
00657 
00658     for (int cut = 0; cut !=numberOfSteps_; ++cut) {  //loop over different transverse energy cuts
00659       double Et =  (*iPho).et();
00660       bool passesCuts = false;
00661 
00662       //sorting the photon into the right Et-dependant folder
00663       if ( useBinning_ && Et > (cut+1)*cutStep_ && ( (Et < (cut+2)*cutStep_)  | (cut == numberOfSteps_-1) ) ){
00664         passesCuts = true;
00665       }
00666       else if ( !useBinning_ && Et > (cut+1)*cutStep_ ){
00667         passesCuts = true;
00668       }
00669 
00670       if (passesCuts){
00671 
00672         //filling isolation variable histograms
00673 
00674         //tracker isolation variables
00675         fill2DHistoVector(h_nTrackIsolSolid_, (*iPho).nTrkSolidConeDR04(), cut,type);
00676         fill2DHistoVector(h_nTrackIsolHollow_,(*iPho).nTrkHollowConeDR04(),cut,type);
00677 
00678         fill2DHistoVector(h_nTrackIsolSolidVsEta_, (*iPho).eta(),(*iPho).nTrkSolidConeDR04(), cut,type);
00679         fill2DHistoVector(p_nTrackIsolSolidVsEta_, (*iPho).eta(),(*iPho).nTrkSolidConeDR04(), cut,type);
00680         fill2DHistoVector(h_nTrackIsolHollowVsEta_,(*iPho).eta(),(*iPho).nTrkHollowConeDR04(),cut,type);
00681         fill2DHistoVector(p_nTrackIsolHollowVsEta_,(*iPho).eta(),(*iPho).nTrkHollowConeDR04(),cut,type);
00682 
00683         fill2DHistoVector(h_nTrackIsolSolidVsEt_,  (*iPho).et(), (*iPho).nTrkSolidConeDR04(), cut,type);
00684         fill2DHistoVector(p_nTrackIsolSolidVsEt_,  (*iPho).et(), (*iPho).nTrkSolidConeDR04(), cut,type);
00685         fill2DHistoVector(h_nTrackIsolHollowVsEt_, (*iPho).et(), (*iPho).nTrkHollowConeDR04(),cut,type);
00686         fill2DHistoVector(p_nTrackIsolHollowVsEt_, (*iPho).et(), (*iPho).nTrkHollowConeDR04(),cut,type);
00687 
00689         fill2DHistoVector(h_trackPtSumSolid_, (*iPho).trkSumPtSolidConeDR04(),cut,type);
00690         fill2DHistoVector(h_trackPtSumHollow_,(*iPho).trkSumPtSolidConeDR04(),cut,type);
00691 
00692         fill2DHistoVector(h_trackPtSumSolidVsEta_, (*iPho).eta(),(*iPho).trkSumPtSolidConeDR04(), cut,type);
00693         fill2DHistoVector(p_trackPtSumSolidVsEta_, (*iPho).eta(),(*iPho).trkSumPtSolidConeDR04(), cut,type);
00694         fill2DHistoVector(h_trackPtSumHollowVsEta_,(*iPho).eta(),(*iPho).trkSumPtHollowConeDR04(),cut,type);
00695         fill2DHistoVector(p_trackPtSumHollowVsEta_,(*iPho).eta(),(*iPho).trkSumPtHollowConeDR04(),cut,type);
00696 
00697         fill2DHistoVector(h_trackPtSumSolidVsEt_,  (*iPho).et(), (*iPho).trkSumPtSolidConeDR04(), cut,type);
00698         fill2DHistoVector(p_trackPtSumSolidVsEt_,  (*iPho).et(), (*iPho).trkSumPtSolidConeDR04(), cut,type);
00699         fill2DHistoVector(h_trackPtSumHollowVsEt_, (*iPho).et(), (*iPho).trkSumPtHollowConeDR04(),cut,type);
00700         fill2DHistoVector(p_trackPtSumHollowVsEt_, (*iPho).et(), (*iPho).trkSumPtHollowConeDR04(),cut,type);
00701         //calorimeter isolation variables
00702 
00703         fill2DHistoVector(h_ecalSum_,(*iPho).ecalRecHitSumEtConeDR04(),cut,type);
00704     if(iPho->isEB()){fill2DHistoVector(h_ecalSumEBarrel_,(*iPho).ecalRecHitSumEtConeDR04(),cut,type);}
00705     if(iPho->isEE()){fill2DHistoVector(h_ecalSumEEndcap_,(*iPho).ecalRecHitSumEtConeDR04(),cut,type);}
00706         fill2DHistoVector(h_ecalSumVsEta_,(*iPho).eta(),(*iPho).ecalRecHitSumEtConeDR04(),cut,type);
00707         fill2DHistoVector(p_ecalSumVsEta_,(*iPho).eta(),(*iPho).ecalRecHitSumEtConeDR04(),cut,type);
00708         fill2DHistoVector(h_ecalSumVsEt_, (*iPho).et(), (*iPho).ecalRecHitSumEtConeDR04(),cut,type);
00709         fill3DHistoVector(p_ecalSumVsEt_, (*iPho).et(), (*iPho).ecalRecHitSumEtConeDR04(),cut,type,part);
00710 
00712 
00713         fill2DHistoVector(h_hcalSum_,(*iPho).hcalTowerSumEtConeDR04(),cut,type);
00714     if(iPho->isEB()){fill2DHistoVector(h_hcalSumEBarrel_,(*iPho).hcalTowerSumEtConeDR04(),cut,type);}
00715     if(iPho->isEE()){fill2DHistoVector(h_hcalSumEEndcap_,(*iPho).hcalTowerSumEtConeDR04(),cut,type);}
00716         fill2DHistoVector(h_hcalSumVsEta_,(*iPho).eta(),(*iPho).hcalTowerSumEtConeDR04(),cut,type);
00717         fill2DHistoVector(p_hcalSumVsEta_,(*iPho).eta(),(*iPho).hcalTowerSumEtConeDR04(),cut,type);
00718         fill2DHistoVector(h_hcalSumVsEt_, (*iPho).et(), (*iPho).hcalTowerSumEtConeDR04(),cut,type);
00719         fill3DHistoVector(p_hcalSumVsEt_, (*iPho).et(), (*iPho).hcalTowerSumEtConeDR04(),cut,type,part);
00720 
00721         fill3DHistoVector(h_hOverE_,(*iPho).hadronicOverEm(),cut,type,part);
00722         fill2DHistoVector(p_hOverEVsEta_,(*iPho).eta(),(*iPho).hadronicOverEm(),cut,type);
00723         fill2DHistoVector(p_hOverEVsEt_, (*iPho).et(), (*iPho).hadronicOverEm(),cut,type);
00724 
00725         fill3DHistoVector(h_h1OverE_,(*iPho).hadronicDepth1OverEm(),cut,type,part);
00726         fill3DHistoVector(h_h2OverE_,(*iPho).hadronicDepth2OverEm(),cut,type,part);
00727 
00728         //filling photon histograms
00729 
00730         nPho[cut][0][0]++;
00731         nPho[cut][0][part]++;
00732         nPho[cut][type][0]++;
00733         nPho[cut][type][part]++;
00734 
00735         //energy variables
00736 
00737         fill3DHistoVector(h_phoE_, (*iPho).energy(),cut,type,part);
00738         fill3DHistoVector(h_phoEt_,(*iPho).et(),    cut,type,part);
00739 
00740         //geometrical variables
00741 
00742         fill2DHistoVector(h_phoEta_,(*iPho).eta(),cut,type);
00743         fill2DHistoVector(h_scEta_, (*iPho).superCluster()->eta(),cut,type);
00744 
00745         fill3DHistoVector(h_phoPhi_,(*iPho).phi(),cut,type,part);
00746         fill3DHistoVector(h_scPhi_, (*iPho).superCluster()->phi(),cut,type,part);
00747 
00748         //shower shape variables
00749 
00750         fill3DHistoVector(h_r9_,(*iPho).r9(),cut,type,part);
00751         fill2DHistoVector(h_r9VsEta_,(*iPho).eta(),(*iPho).r9(),cut,type);
00752         fill2DHistoVector(p_r9VsEta_,(*iPho).eta(),(*iPho).r9(),cut,type);
00753         fill2DHistoVector(h_r9VsEt_, (*iPho).et(), (*iPho).r9(),cut,type);
00754         fill2DHistoVector(p_r9VsEt_, (*iPho).et(), (*iPho).r9(),cut,type);
00755 
00756         fill2DHistoVector(h_e1x5VsEta_,(*iPho).eta(),(*iPho).e1x5(),cut,type);
00757         fill2DHistoVector(p_e1x5VsEta_,(*iPho).eta(),(*iPho).e1x5(),cut,type);
00758         fill2DHistoVector(h_e1x5VsEt_, (*iPho).et(), (*iPho).e1x5(),cut,type);
00759         fill2DHistoVector(p_e1x5VsEt_, (*iPho).et(), (*iPho).e1x5(),cut,type);
00760 
00761         fill2DHistoVector(h_e2x5VsEta_,(*iPho).eta(),(*iPho).e2x5(),cut,type);
00762         fill2DHistoVector(p_e2x5VsEta_,(*iPho).eta(),(*iPho).e2x5(),cut,type);
00763         fill2DHistoVector(h_e2x5VsEt_, (*iPho).et(), (*iPho).e2x5(),cut,type);
00764         fill2DHistoVector(p_e2x5VsEt_, (*iPho).et(), (*iPho).e2x5(),cut,type);
00765 
00766         fill2DHistoVector(h_maxEXtalOver3x3VsEta_,(*iPho).eta(),(*iPho).maxEnergyXtal()/(*iPho).e3x3(),cut,type);
00767         fill2DHistoVector(p_maxEXtalOver3x3VsEta_,(*iPho).eta(),(*iPho).maxEnergyXtal()/(*iPho).e3x3(),cut,type);
00768         fill2DHistoVector(h_maxEXtalOver3x3VsEt_, (*iPho).et(), (*iPho).maxEnergyXtal()/(*iPho).e3x3(),cut,type);
00769         fill2DHistoVector(p_maxEXtalOver3x3VsEt_, (*iPho).et(), (*iPho).maxEnergyXtal()/(*iPho).e3x3(),cut,type);
00770 
00771 
00772         fill2DHistoVector(h_r1x5VsEta_,(*iPho).eta(),(*iPho).r1x5(),cut,type);
00773         fill2DHistoVector(p_r1x5VsEta_,(*iPho).eta(),(*iPho).r1x5(),cut,type);
00774         fill2DHistoVector(h_r1x5VsEt_, (*iPho).et(), (*iPho).r1x5(),cut,type);
00775         fill2DHistoVector(p_r1x5VsEt_, (*iPho).et(), (*iPho).r1x5(),cut,type);
00776 
00777         fill2DHistoVector(h_r2x5VsEta_,(*iPho).eta(),(*iPho).r2x5(),cut,type);
00778         fill2DHistoVector(p_r2x5VsEta_,(*iPho).eta(),(*iPho).r2x5(),cut,type);
00779         fill2DHistoVector(h_r2x5VsEt_, (*iPho).et(), (*iPho).r2x5(),cut,type);
00780         fill2DHistoVector(p_r2x5VsEt_, (*iPho).et(), (*iPho).r2x5(),cut,type);
00781 
00782         fill3DHistoVector(h_phoSigmaIetaIeta_,(*iPho).sigmaIetaIeta(),cut,type,part);
00783         fill2DHistoVector(h_sigmaIetaIetaVsEta_,(*iPho).eta(),(*iPho).sigmaIetaIeta(),cut,type);
00784         fill2DHistoVector(p_sigmaIetaIetaVsEta_,(*iPho).eta(),(*iPho).sigmaIetaIeta(),cut,type);
00785 
00786 
00787 
00788         //filling histograms for photons containing a bad ECAL channel
00789 
00790         bool atLeastOneDeadChannel=false;
00791         for(reco::CaloCluster_iterator bcIt = (*iPho).superCluster()->clustersBegin();bcIt != (*iPho).superCluster()->clustersEnd(); ++bcIt) { //loop over basic clusters in SC
00792           for(vector< pair<DetId, float> >::const_iterator rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) { //loop over rec hits in basic cluster
00793 
00794             for(EcalRecHitCollection::const_iterator it = ecalRecHitCollection.begin(); it !=  ecalRecHitCollection.end(); ++it) { //loop over all rec hits to find the right ones
00795               if  (rhIt->first ==  (*it).id() ) { //found the matching rechit
00796                 if (  (*it).recoFlag() == 9 ) { //has a bad channel
00797                   atLeastOneDeadChannel=true;
00798                   break;
00799                 }
00800               }
00801             }
00802           }
00803         }
00804         if ( atLeastOneDeadChannel ) {
00805           fill2DHistoVector(h_phoPhi_BadChannels_,(*iPho).phi(),cut,type);
00806           fill2DHistoVector(h_phoEta_BadChannels_,(*iPho).eta(),cut,type);
00807           fill2DHistoVector(h_phoEt_BadChannels_, (*iPho).et(), cut,type);
00808         }
00809 
00810 
00811         // filling conversion-related histograms
00812         if((*iPho).hasConversionTracks()){
00813           nConv[cut][0][0]++;
00814           nConv[cut][0][part]++;
00815           nConv[cut][type][0]++;
00816           nConv[cut][type][part]++;
00817         }
00818 
00819         //loop over conversions (don't forget, we're still inside the photon loop,
00820         // i.e. these are all the conversions for this ONE photon, not for all the photons in the event)
00821         reco::ConversionRefVector conversions = (*iPho).conversions();
00822         for (unsigned int iConv=0; iConv<conversions.size(); iConv++) {
00823 
00824           reco::ConversionRef aConv=conversions[iConv];
00825 
00826           if ( aConv->nTracks() <2 ) continue;
00827 
00828           //fill histogram for denominator of vertex reconstruction efficiency plot
00829           if(cut==0) h_phoEta_Vertex_->Fill(aConv->refittedPairMomentum().eta());
00830 
00831           if ( !(aConv->conversionVertex().isValid()) ) continue;
00832 
00833           float chi2Prob = ChiSquaredProbability( aConv->conversionVertex().chi2(), aConv->conversionVertex().ndof() );
00834 
00835           if(chi2Prob<0.0005) continue;
00836 
00837           fill2DHistoVector(h_vertexChi2Prob_,chi2Prob,cut,type);
00838 
00839 
00840 
00841           fill3DHistoVector(h_phoConvE_, (*iPho).energy(),cut,type,part);
00842           fill3DHistoVector(h_phoConvEt_,(*iPho).et(),cut,type,part);
00843           fill3DHistoVector(h_phoConvR9_,(*iPho).r9(),cut,type,part);
00844 
00845           if (cut==0 && isLoosePhoton){
00846             h_convEta_Loose_->Fill((*iPho).eta());
00847             h_convEt_Loose_->Fill( (*iPho).et() );
00848           }
00849           if (cut==0 && isTightPhoton){
00850             h_convEta_Tight_->Fill((*iPho).eta());
00851             h_convEt_Tight_->Fill( (*iPho).et() );
00852           }
00853 
00854           fill2DHistoVector(h_phoConvEta_,aConv->refittedPairMomentum().eta(),cut,type);
00855           fill3DHistoVector(h_phoConvPhi_,aConv->refittedPairMomentum().phi(),cut,type,part);
00856 
00857         
00858           //we use the photon position because we'll be dividing it by a photon histogram (not a conversion histogram)
00859           fill2DHistoVector(h_phoConvEtaForEfficiency_,(*iPho).eta(),cut,type);
00860           fill3DHistoVector(h_phoConvPhiForEfficiency_,(*iPho).phi(),cut,type,part);
00861 
00862 
00863           //vertex histograms
00864           double convR= sqrt(aConv->conversionVertex().position().perp2());
00865           double scalar = aConv->conversionVertex().position().x()*aConv->refittedPairMomentum().x() + aConv->conversionVertex().position().y()*aConv->refittedPairMomentum().y();
00866           if ( scalar < 0 ) convR= -convR;
00867 
00868           fill2DHistoVector(h_convVtxRvsZ_,aConv->conversionVertex().position().z(), convR,cut,type);//trying to "see" R-Z view of tracker
00869           fill2DHistoVector(h_convVtxZ_,aConv->conversionVertex().position().z(), cut,type);
00870 
00871           if(fabs(aConv->caloCluster()[0]->eta()) > 1.5){//trying to "see" tracker endcaps
00872             fill2DHistoVector(h_convVtxZEndcap_,aConv->conversionVertex().position().z(), cut,type);
00873           }
00874           else if(fabs(aConv->caloCluster()[0]->eta()) < 1){//trying to "see" tracker barrel
00875             fill2DHistoVector(h_convVtxR_,convR,cut,type);
00876             fill2DHistoVector(h_convVtxYvsX_,aConv->conversionVertex().position().x(),aConv->conversionVertex().position().y(),cut,type);
00877           }
00878 
00879           const std::vector<edm::RefToBase<reco::Track> > tracks = aConv->tracks();
00880 
00881 
00882           for (unsigned int i=0; i<tracks.size(); i++) {
00883             fill2DHistoVector(h_tkChi2_,tracks[i]->normalizedChi2(),cut,type);
00884             fill2DHistoVector(p_tkChi2VsEta_,aConv->caloCluster()[0]->eta(),tracks[i]->normalizedChi2(),cut,type);
00885             fill2DHistoVector(p_dCotTracksVsEta_,aConv->caloCluster()[0]->eta(),aConv->pairCotThetaSeparation(),cut,type);
00886             fill2DHistoVector(p_nHitsVsEta_,aConv->caloCluster()[0]->eta(),float(tracks[i]->numberOfValidHits()),cut,type);
00887           }
00888 
00889           //calculating delta eta and delta phi of the two tracks
00890 
00891           float  DPhiTracksAtVtx = -99;
00892           float  dPhiTracksAtEcal= -99;
00893           float  dEtaTracksAtEcal= -99;
00894 
00895           float phiTk1= aConv->tracksPin()[0].phi();
00896           float phiTk2= aConv->tracksPin()[1].phi();
00897           DPhiTracksAtVtx = phiTk1-phiTk2;
00898           DPhiTracksAtVtx = phiNormalization( DPhiTracksAtVtx );
00899 
00900           if (aConv->bcMatchingWithTracks()[0].isNonnull() && aConv->bcMatchingWithTracks()[1].isNonnull() ) {
00901             float recoPhi1 = aConv->ecalImpactPosition()[0].phi();
00902             float recoPhi2 = aConv->ecalImpactPosition()[1].phi();
00903             float recoEta1 = aConv->ecalImpactPosition()[0].eta();
00904             float recoEta2 = aConv->ecalImpactPosition()[1].eta();
00905 
00906             recoPhi1 = phiNormalization(recoPhi1);
00907             recoPhi2 = phiNormalization(recoPhi2);
00908 
00909             dPhiTracksAtEcal = recoPhi1 -recoPhi2;
00910             dPhiTracksAtEcal = phiNormalization( dPhiTracksAtEcal );
00911             dEtaTracksAtEcal = recoEta1 -recoEta2;
00912 
00913           }
00914 
00915         
00916           fill3DHistoVector(h_dPhiTracksAtVtx_,DPhiTracksAtVtx,cut,type,part);
00917           fill3DHistoVector(h_dPhiTracksAtEcal_,fabs(dPhiTracksAtEcal),cut,type,part);
00918           fill3DHistoVector(h_dEtaTracksAtEcal_, dEtaTracksAtEcal,cut,type,part);
00919           fill3DHistoVector(h_eOverPTracks_,aConv->EoverPrefittedTracks(),cut,type,part);
00920           fill3DHistoVector(h_pOverETracks_,1./aConv->EoverPrefittedTracks(),cut,type,part);
00921           fill3DHistoVector(h_dCotTracks_,aConv->pairCotThetaSeparation(),cut,type,part);
00922 
00923 
00924         }//end loop over conversions
00925 
00926       }//end loop over photons passing cuts
00927     }//end loop over transverse energy cuts
00928 
00929 
00930 
00931 
00932 
00933     //make invariant mass plots
00934 
00935     if (isIsolated && iPho->et()>=invMassEtCut_){
00936 
00937       for (reco::PhotonCollection::const_iterator iPho2=iPho+1; iPho2!=photonCollection.end(); iPho2++){
00938 
00939         edm::Ref<reco::PhotonCollection> photonref2(photonHandle, photonCounter); //note: it's correct to use photonCounter and not photonCounter+1
00940                                                                                   //since it has already been incremented earlier
00941 
00942         bool  isTightPhoton2(false), isLoosePhoton2(false);
00943 
00944         if ( !isHeavyIon_ ) {
00945           isTightPhoton2 = (tightPhotonID)[photonref2];
00946           isLoosePhoton2 = (loosePhotonID)[photonref2];
00947         }
00948 
00949         bool isIsolated2=false;
00950         if ( isolationStrength_ == 0)  isIsolated2 = isLoosePhoton2;
00951         if ( isolationStrength_ == 1)  isIsolated2 = isTightPhoton2;
00952 
00953         reco::ConversionRefVector conversions = (*iPho).conversions();
00954         reco::ConversionRefVector conversions2 = (*iPho2).conversions();
00955 
00956         if(isIsolated2 && iPho2->et()>=invMassEtCut_){
00957 
00958           math::XYZTLorentzVector p12 = iPho->p4()+iPho2->p4();
00959           float gamgamMass2 = p12.Dot(p12);
00960 
00961 
00962           h_invMassAllPhotons_ -> Fill(sqrt( gamgamMass2 ));
00963       if(iPho->isEB() && iPho2->isEB()){h_invMassPhotonsEBarrel_ -> Fill(sqrt( gamgamMass2 ));}
00964       if(iPho->isEE() || iPho2->isEE()){h_invMassPhotonsEEndcap_ -> Fill(sqrt( gamgamMass2 ));}      
00965 
00966           if(conversions.size()!=0 && conversions[0]->nTracks() >= 2){
00967             if(conversions2.size()!=0 && conversions2[0]->nTracks() >= 2) h_invMassTwoWithTracks_ -> Fill(sqrt( gamgamMass2 ));
00968             else h_invMassOneWithTracks_ -> Fill(sqrt( gamgamMass2 ));
00969           }
00970           else if(conversions2.size()!=0 && conversions2[0]->nTracks() >= 2) h_invMassOneWithTracks_ -> Fill(sqrt( gamgamMass2 ));
00971           else h_invMassZeroWithTracks_ -> Fill(sqrt( gamgamMass2 ));
00972         }
00973 
00974       }
00975 
00976     }
00977 
00978 
00979 
00980   }
00981 
00982 
00983   //filling number of photons/conversions per event histograms
00984   for (int cut = 0; cut !=numberOfSteps_; ++cut) {
00985     for(uint type=0;type!=types_.size();++type){
00986       for(uint part=0;part!=parts_.size();++part){
00987         h_nPho_[cut][type][part]-> Fill (float(nPho[cut][type][part]));
00988         h_nConv_[cut][type][part]-> Fill (float(nConv[cut][type][part]));
00989       }
00990     }
00991   }
00992 
00993 }//End of Analyze method
00994 
00995 void PhotonAnalyzer::endRun(const edm::Run& run, const edm::EventSetup& setup)
00996 {
00997   if(!standAlone_){
00998 
00999     dbe_->setCurrentFolder("Egamma/PhotonAnalyzer");
01000 
01001     //keep track of how many histos are in each folder
01002     totalNumberOfHistos_efficiencyFolder->Fill(histo_index_efficiency_);
01003     totalNumberOfHistos_invMassFolder->Fill(histo_index_invMass_);
01004     totalNumberOfHistos_photonsFolder->Fill(histo_index_photons_);
01005     totalNumberOfHistos_conversionsFolder->Fill(histo_index_conversions_);
01006 
01007   }
01008 
01009 }
01010 
01011 
01012 void PhotonAnalyzer::endJob()
01013 {
01014   //dbe_->showDirStructure();
01015   if(standAlone_){
01016     dbe_->setCurrentFolder("Egamma/PhotonAnalyzer");
01017 
01018     //keep track of how many histos are in each folder
01019     totalNumberOfHistos_efficiencyFolder->Fill(histo_index_efficiency_);
01020     totalNumberOfHistos_invMassFolder->Fill(histo_index_invMass_);
01021     totalNumberOfHistos_photonsFolder->Fill(histo_index_photons_);
01022     totalNumberOfHistos_conversionsFolder->Fill(histo_index_conversions_);
01023 
01024 
01025     dbe_->save(outputFileName_);
01026   }
01027 
01028 
01029 }
01030 
01032 
01033 
01034 
01035 float PhotonAnalyzer::phiNormalization(float & phi)
01036 {
01037  const float PI    = 3.1415927;
01038  const float TWOPI = 2.0*PI;
01039 
01040  if(phi >  PI) {phi = phi - TWOPI;}
01041  if(phi < -PI) {phi = phi + TWOPI;}
01042 
01043  return phi;
01044 }
01045 
01046 
01047 void  PhotonAnalyzer::fill2DHistoVector(vector<vector<MonitorElement*> >& histoVector,double x, double y, int cut, int type){
01048 
01049   histoVector[cut][0]->Fill(x,y);
01050   if(histoVector[cut].size()>1)   histoVector[cut][type]->Fill(x,y); //don't try to fill 2D histos that are only in the "AllPhotons" folder
01051 
01052 }
01053 
01054 void  PhotonAnalyzer::fill2DHistoVector(vector<vector<MonitorElement*> >& histoVector, double x, int cut, int type){
01055 
01056   histoVector[cut][0]->Fill(x);
01057   histoVector[cut][type]->Fill(x);
01058 
01059 }
01060 
01061 void  PhotonAnalyzer::fill3DHistoVector(vector<vector<vector<MonitorElement*> > >& histoVector,double x, int cut, int type, int part){
01062 
01063   histoVector[cut][0][0]->Fill(x);
01064   histoVector[cut][0][part]->Fill(x);
01065   histoVector[cut][type][0]->Fill(x);
01066   histoVector[cut][type][part]->Fill(x);
01067 
01068 }
01069 
01070 void  PhotonAnalyzer::fill3DHistoVector(vector<vector<vector<MonitorElement*> > >& histoVector,double x, double y, int cut, int type, int part){
01071 
01072   histoVector[cut][0][0]->Fill(x,y);
01073   histoVector[cut][0][part]->Fill(x,y);
01074   histoVector[cut][type][0]->Fill(x,y);
01075   histoVector[cut][type][part]->Fill(x,y);
01076 
01077 }
01078 
01079 
01080 MonitorElement* PhotonAnalyzer::bookHisto(string histoName, string title, int bin, double min, double max)
01081 {
01082 
01083   int histo_index = 0;
01084   stringstream histo_number_stream;
01085 
01086   //determining which folder we're in
01087   if(dbe_->pwd().find( "InvMass" ) != string::npos){
01088     histo_index_invMass_++;
01089     histo_index = histo_index_invMass_;
01090   }
01091   if(dbe_->pwd().find( "Efficiencies" ) != string::npos){
01092     histo_index_efficiency_++;
01093     histo_index = histo_index_efficiency_;
01094   }
01095 
01096   histo_number_stream << "h_";
01097   if(histo_index<10)   histo_number_stream << "0";
01098   histo_number_stream << histo_index;
01099 
01100   return dbe_->book1D(histo_number_stream.str()+"_"+histoName,title,bin,min,max);
01101 
01102 }
01103 
01104 
01105 void PhotonAnalyzer::book2DHistoVector(vector<vector<MonitorElement*> > &temp2DVector,
01106                                        string histoType, string histoName, string title,
01107                                                                              int xbin, double xmin,double xmax,
01108                                                                              int ybin, double ymin, double ymax)
01109 {
01110   int histo_index = 0;
01111 
01112   vector<MonitorElement*> temp1DVector;
01113 //   vector<vector<MonitorElement*> > temp2DVector;
01114 
01115   //determining which folder we're in
01116   bool conversionPlot = false;
01117   if(dbe_->pwd().find( "Conversions" ) != string::npos) conversionPlot = true;
01118   bool TwoDPlot = false;
01119   if(histoName.find( "2D" ) != string::npos) TwoDPlot = true;
01120 
01121   if(conversionPlot){
01122     histo_index_conversions_++;
01123     histo_index = histo_index_conversions_;
01124   }
01125   else{
01126     histo_index_photons_++;
01127     histo_index = histo_index_photons_;
01128   }
01129 
01130   stringstream histo_number_stream;
01131   histo_number_stream << "h_";
01132   if(histo_index<10)   histo_number_stream << "0";
01133   histo_number_stream << histo_index << "_";
01134 
01135 
01136 
01137   for(int cut = 0; cut != numberOfSteps_; ++cut){ //looping over Et cut values
01138 
01139     for(uint type=0;type!=types_.size();++type){  //looping over isolation type
01140 
01141       currentFolder_.str("");
01142       currentFolder_ << "Egamma/PhotonAnalyzer/" << types_[type] << "Photons/Et above " << (cut+1)*cutStep_ << " GeV";
01143       if(conversionPlot) currentFolder_ << "/Conversions";
01144 
01145       dbe_->setCurrentFolder(currentFolder_.str());
01146 
01147       string kind;
01148       if(conversionPlot) kind = " Conversions: ";
01149       else kind = " Photons: ";
01150 
01151       if(histoType=="1D")           temp1DVector.push_back(dbe_->book1D(      histo_number_stream.str()+histoName,types_[type]+kind+title,xbin,xmin,xmax));
01152       else if(histoType=="2D"){
01153         if((TwoDPlot && type==0) || !TwoDPlot){//only book the 2D plots in the "AllPhotons" folder
01154                                     temp1DVector.push_back(dbe_->book2D(      histo_number_stream.str()+histoName,types_[type]+kind+title,xbin,xmin,xmax,ybin,ymin,ymax));
01155         }
01156       }
01157       else if(histoType=="Profile") temp1DVector.push_back(dbe_->bookProfile( histo_number_stream.str()+histoName,types_[type]+kind+title,xbin,xmin,xmax,ybin,ymin,ymax,""));
01158       else cout << "bad histoType\n";
01159     }
01160 
01161     temp2DVector.push_back(temp1DVector);
01162     temp1DVector.clear();
01163   }
01164 
01165 //   return temp2DVector;
01166 
01167 }
01168 
01169 
01170 void PhotonAnalyzer::book3DHistoVector(vector<vector<vector<MonitorElement*> > > &temp3DVector,
01171                                        string histoType, string histoName, string title,
01172                                                                              int xbin, double xmin,double xmax,
01173                                                                              int ybin, double ymin, double ymax)
01174 {
01175 
01176 
01177   int histo_index = 0;
01178 
01179   vector<MonitorElement*> temp1DVector;
01180   vector<vector<MonitorElement*> > temp2DVector;
01181 //   vector<vector<vector<MonitorElement*> > > temp3DVector;
01182 
01183 
01184   //determining which folder we're in
01185   bool conversionPlot = false;
01186   if(dbe_->pwd().find( "Conversions" ) != string::npos) conversionPlot = true;
01187 
01188 
01189   if(conversionPlot){
01190     histo_index_conversions_++;
01191     histo_index = histo_index_conversions_;
01192   }
01193   else{
01194     histo_index_photons_++;
01195     histo_index = histo_index_photons_;
01196   }
01197 
01198 
01199 
01200   stringstream histo_number_stream;
01201   histo_number_stream << "h_";
01202   if(histo_index<10)   histo_number_stream << "0";
01203   histo_number_stream << histo_index << "_";
01204 
01205   for(int cut = 0; cut != numberOfSteps_; ++cut){     //looping over Et cut values
01206 
01207     for(uint type=0;type!=types_.size();++type){      //looping over isolation type
01208 
01209       for(uint part=0;part!=parts_.size();++part){    //looping over different parts of the ecal
01210 
01211         currentFolder_.str("");
01212         currentFolder_ << "Egamma/PhotonAnalyzer/" << types_[type] << "Photons/Et above " << (cut+1)*cutStep_ << " GeV";
01213         if(conversionPlot) currentFolder_ << "/Conversions";
01214 
01215         dbe_->setCurrentFolder(currentFolder_.str());
01216 
01217         string kind;
01218         if(conversionPlot) kind = " Conversions: ";
01219         else kind = " Photons: ";
01220 
01221         if(histoType=="1D")           temp1DVector.push_back(dbe_->book1D(      histo_number_stream.str()+histoName+parts_[part],types_[type]+kind+parts_[part]+": "+title,xbin,xmin,xmax));
01222         else if(histoType=="2D")      temp1DVector.push_back(dbe_->book2D(      histo_number_stream.str()+histoName+parts_[part],types_[type]+kind+parts_[part]+": "+title,xbin,xmin,xmax,ybin,ymin,ymax));
01223         else if(histoType=="Profile") temp1DVector.push_back(dbe_->bookProfile( histo_number_stream.str()+histoName+parts_[part],types_[type]+kind+parts_[part]+": "+title,xbin,xmin,xmax,ybin,ymin,ymax,""));
01224         else cout << "bad histoType\n";
01225 
01226 
01227       }
01228 
01229       temp2DVector.push_back(temp1DVector);
01230       temp1DVector.clear();
01231     }
01232 
01233     temp3DVector.push_back(temp2DVector);
01234     temp2DVector.clear();
01235   }
01236 
01237   //  return temp3DVector;
01238 }