CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PhotonAnalyzer.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 //
4 
5 #include "DQMOffline/EGamma/plugins/PhotonAnalyzer.h"
6 
7 
18 using namespace std;
19 
20 
21 PhotonAnalyzer::PhotonAnalyzer( const edm::ParameterSet& pset )
22 {
23 
24  fName_ = pset.getParameter<string>("analyzerName");
25  verbosity_ = pset.getUntrackedParameter<int>("Verbosity");
26 
27  prescaleFactor_ = pset.getUntrackedParameter<int>("prescaleFactor",1);
28 
29  photon_token_ = consumes<vector<reco::Photon> >(pset.getParameter<edm::InputTag>("phoProducer"));
30 
31  barrelRecHit_token_ = consumes<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >(pset.getParameter<edm::InputTag>("barrelRecHitProducer"));
32 
33  PhotonIDLoose_token_ = consumes<edm::ValueMap<bool> >(pset.getParameter<edm::InputTag>("photonIDLoose"));
34  PhotonIDTight_token_ = consumes<edm::ValueMap<bool> >(pset.getParameter<edm::InputTag>("photonIDTight"));
35 
36  endcapRecHit_token_ = consumes<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >(pset.getParameter<edm::InputTag>("endcapRecHitProducer"));
37 
38  triggerEvent_token_ = consumes<trigger::TriggerEvent>(pset.getParameter<edm::InputTag>("triggerEvent"));
39 
40  offline_pvToken_ = consumes<reco::VertexCollection>(pset.getUntrackedParameter<edm::InputTag>("offlinePV", edm::InputTag("offlinePrimaryVertices")));
41 
42 
43  minPhoEtCut_ = pset.getParameter<double>("minPhoEtCut");
44  photonMaxEta_ = pset.getParameter<double>("maxPhoEta");
45  invMassEtCut_ = pset.getParameter<double>("invMassEtCut");
46  cutStep_ = pset.getParameter<double>("cutStep");
47  numberOfSteps_ = pset.getParameter<int>("numberOfSteps");
48 
49  useBinning_ = pset.getParameter<bool>("useBinning");
50  useTriggerFiltering_ = pset.getParameter<bool>("useTriggerFiltering");
51 
52  minimalSetOfHistos_ = pset.getParameter<bool>("minimalSetOfHistos");
53  excludeBkgHistos_ = pset.getParameter<bool>("excludeBkgHistos");
54 
55  standAlone_ = pset.getParameter<bool>("standAlone");
56  outputFileName_ = pset.getParameter<string>("OutputFileName");
57 
58  isolationStrength_ = pset.getParameter<int>("isolationStrength");
59 
60  isHeavyIon_ = pset.getUntrackedParameter<bool>("isHeavyIon",false);
61 
62  parameters_ = pset;
63 
64  histo_index_photons_ = 0;
65  histo_index_conversions_ = 0;
66  histo_index_efficiency_ = 0;
67  histo_index_invMass_ = 0;
68 }
69 
70 
71 
72 PhotonAnalyzer::~PhotonAnalyzer() {}
73 
74 
75 void PhotonAnalyzer::beginJob()
76 {
77 
78  nEvt_=0;
79  nEntry_=0;
80 
81  dbe_ = 0;
82  dbe_ = edm::Service<DQMStore>().operator->();
83 
84 
85 
86  double eMin = parameters_.getParameter<double>("eMin");
87  double eMax = parameters_.getParameter<double>("eMax");
88  int eBin = parameters_.getParameter<int>("eBin");
89 
90  double etMin = parameters_.getParameter<double>("etMin");
91  double etMax = parameters_.getParameter<double>("etMax");
92  int etBin = parameters_.getParameter<int>("etBin");
93 
94  double sumMin = parameters_.getParameter<double>("sumMin");
95  double sumMax = parameters_.getParameter<double>("sumMax");
96  int sumBin = parameters_.getParameter<int>("sumBin");
97 
98  double etaMin = parameters_.getParameter<double>("etaMin");
99  double etaMax = parameters_.getParameter<double>("etaMax");
100  int etaBin = parameters_.getParameter<int>("etaBin");
101 
102  double phiMin = parameters_.getParameter<double>("phiMin");
103  double phiMax = parameters_.getParameter<double>("phiMax");
104  int phiBin = parameters_.getParameter<int>("phiBin");
105 
106  double r9Min = parameters_.getParameter<double>("r9Min");
107  double r9Max = parameters_.getParameter<double>("r9Max");
108  int r9Bin = parameters_.getParameter<int>("r9Bin");
109 
110  double hOverEMin = parameters_.getParameter<double>("hOverEMin");
111  double hOverEMax = parameters_.getParameter<double>("hOverEMax");
112  int hOverEBin = parameters_.getParameter<int>("hOverEBin");
113 
114  double xMin = parameters_.getParameter<double>("xMin");
115  double xMax = parameters_.getParameter<double>("xMax");
116  int xBin = parameters_.getParameter<int>("xBin");
117 
118  double yMin = parameters_.getParameter<double>("yMin");
119  double yMax = parameters_.getParameter<double>("yMax");
120  int yBin = parameters_.getParameter<int>("yBin");
121 
122  double numberMin = parameters_.getParameter<double>("numberMin");
123  double numberMax = parameters_.getParameter<double>("numberMax");
124  int numberBin = parameters_.getParameter<int>("numberBin");
125 
126  double zMin = parameters_.getParameter<double>("zMin");
127  double zMax = parameters_.getParameter<double>("zMax");
128  int zBin = parameters_.getParameter<int>("zBin");
129 
130  double rMin = parameters_.getParameter<double>("rMin");
131  double rMax = parameters_.getParameter<double>("rMax");
132  int rBin = parameters_.getParameter<int>("rBin");
133 
134  double dPhiTracksMin = parameters_.getParameter<double>("dPhiTracksMin");
135  double dPhiTracksMax = parameters_.getParameter<double>("dPhiTracksMax");
136  int dPhiTracksBin = parameters_.getParameter<int>("dPhiTracksBin");
137 
138  double dEtaTracksMin = parameters_.getParameter<double>("dEtaTracksMin");
139  double dEtaTracksMax = parameters_.getParameter<double>("dEtaTracksMax");
140  int dEtaTracksBin = parameters_.getParameter<int>("dEtaTracksBin");
141 
142  double sigmaIetaMin = parameters_.getParameter<double>("sigmaIetaMin");
143  double sigmaIetaMax = parameters_.getParameter<double>("sigmaIetaMax");
144  int sigmaIetaBin = parameters_.getParameter<int>("sigmaIetaBin");
145 
146  double eOverPMin = parameters_.getParameter<double>("eOverPMin");
147  double eOverPMax = parameters_.getParameter<double>("eOverPMax");
148  int eOverPBin = parameters_.getParameter<int>("eOverPBin");
149 
150  double chi2Min = parameters_.getParameter<double>("chi2Min");
151  double chi2Max = parameters_.getParameter<double>("chi2Max");
152  int chi2Bin = parameters_.getParameter<int>("chi2Bin");
153 
154 
155  int reducedEtBin = etBin/4;
156  int reducedEtaBin = etaBin/4;
157  int reducedSumBin = sumBin/4;
158  int reducedR9Bin = r9Bin/4;
159 
160 
161  parts_.push_back("AllEcal");
162  parts_.push_back("Barrel");
163  parts_.push_back("Endcaps");
164 
165  types_.push_back("All");
166  types_.push_back("GoodCandidate");
167  if (!excludeBkgHistos_) types_.push_back("Background");
168 
169 
170 
172 
173  if (dbe_) {
174 
175  dbe_->setCurrentFolder("Egamma/"+fName_+"/");
176 
177  //int values stored in MEs to keep track of how many histograms are in each folder
178  totalNumberOfHistos_efficiencyFolder = dbe_->bookInt("numberOfHistogramsInEfficiencyFolder");
179  totalNumberOfHistos_photonsFolder = dbe_->bookInt("numberOfHistogramsInPhotonsFolder");
180  totalNumberOfHistos_conversionsFolder = dbe_->bookInt("numberOfHistogramsInConversionsFolder");
181  totalNumberOfHistos_invMassFolder = dbe_->bookInt("numberOfHistogramsInInvMassFolder");
182 
183 
184  //Efficiency histograms
185 
186  dbe_->setCurrentFolder("Egamma/"+fName_+"/Efficiencies");
187 
188  //don't number these histograms with the "bookHisto" method, since they'll be erased in the offline client
189  h_phoEta_Loose_ = dbe_->book1D("phoEtaLoose","Loose Photon #eta",etaBin,etaMin,etaMax);
190  h_phoEta_Tight_ = dbe_->book1D("phoEtaTight","Tight Photon #eta",etaBin,etaMin,etaMax);
191  h_phoEt_Loose_ = dbe_->book1D("phoEtLoose", "Loose Photon E_{T}",etBin,etMin,etMax);
192  h_phoEt_Tight_ = dbe_->book1D("phoEtTight", "Tight Photon E_{T}",etBin,etMin,etMax);
193 
194 
195  h_phoEta_preHLT_ = dbe_->book1D("phoEtaPreHLT", "Photon #eta: before HLT",etaBin,etaMin,etaMax);
196  h_phoEta_postHLT_ = dbe_->book1D("phoEtaPostHLT","Photon #eta: after HLT",etaBin,etaMin,etaMax);
197  h_phoEt_preHLT_ = dbe_->book1D("phoEtPreHLT", "Photon E_{T}: before HLT",etBin,etMin,etMax);
198  h_phoEt_postHLT_ = dbe_->book1D("phoEtPostHLT", "Photon E_{T}: after HLT",etBin,etMin,etMax);
199 
200  h_convEta_Loose_ = dbe_->book1D("convEtaLoose","Converted Loose Photon #eta",etaBin,etaMin,etaMax);
201  h_convEta_Tight_ = dbe_->book1D("convEtaTight","Converted Tight Photon #eta",etaBin,etaMin,etaMax);
202  h_convEt_Loose_ = dbe_->book1D("convEtLoose", "Converted Loose Photon E_{T}",etBin,etMin,etMax);
203  h_convEt_Tight_ = dbe_->book1D("convEtTight", "Converted Tight Photon E_{T}",etBin,etMin,etMax);
204 
205  h_phoEta_Vertex_ = dbe_->book1D("phoEtaVertex","Converted Photons before valid vertex cut: #eta",etaBin,etaMin,etaMax);
206 
207 
208  vector<MonitorElement*> temp1DVectorEta;
209  vector<MonitorElement*> temp1DVectorPhi;
210  vector<vector<MonitorElement*> > temp2DVectorPhi;
211 
212 
213  for(int cut = 0; cut != numberOfSteps_; ++cut){ //looping over Et cut values
214  for(uint type=0;type!=types_.size();++type){ //looping over isolation type
215  currentFolder_.str("");
216  currentFolder_ << "Egamma/"+fName_+"/" << types_[type] << "Photons/Et above " << (cut+1)*cutStep_ << " GeV/Conversions";
217  dbe_->setCurrentFolder(currentFolder_.str());
218 
219  temp1DVectorEta.push_back(dbe_->book1D("phoConvEtaForEfficiency","Converted Photon #eta;#eta",etaBin,etaMin,etaMax));
220  for(uint part=0;part!=parts_.size();++part){
221  temp1DVectorPhi.push_back(dbe_->book1D("phoConvPhiForEfficiency"+parts_[part],"Converted Photon #phi;#phi",phiBin,phiMin,phiMax));
222  }
223  temp2DVectorPhi.push_back(temp1DVectorPhi);
224  temp1DVectorPhi.clear();
225  }
226  h_phoConvEtaForEfficiency_.push_back(temp1DVectorEta);
227  temp1DVectorEta.clear();
228  h_phoConvPhiForEfficiency_.push_back(temp2DVectorPhi);
229  temp2DVectorPhi.clear();
230  }
231 
232 
233 
234 
235  //Invariant mass plots
236 
237  dbe_->setCurrentFolder("Egamma/"+fName_+"/InvMass");
238 
239  h_invMassAllPhotons_ = bookHisto("invMassAllIsolatedPhotons","Two photon invariant mass: All isolated photons;M (GeV)",etBin,etMin,etMax);
240  h_invMassPhotonsEBarrel_ = bookHisto("invMassIsoPhotonsEBarrel", "Two photon invariant mass: isolated photons in barrel; M (GeV)",etBin,etMin,etMax);
241  h_invMassPhotonsEEndcap_ = bookHisto("invMassIsoPhotonsEEndcap", "Two photon invariant mass: isolated photons in endcap; M (GeV)",etBin,etMin,etMax);
242 
243  h_invMassZeroWithTracks_= bookHisto("invMassZeroWithTracks", "Two photon invariant mass: Neither has tracks;M (GeV)", etBin,etMin,etMax);
244  h_invMassOneWithTracks_ = bookHisto("invMassOneWithTracks", "Two photon invariant mass: Only one has tracks;M (GeV)", etBin,etMin,etMax);
245  h_invMassTwoWithTracks_ = bookHisto("invMassTwoWithTracks", "Two photon invariant mass: Both have tracks;M (GeV)", etBin,etMin,etMax);
246 
247 
248  h_nRecoVtx_ = bookHisto("nOfflineVtx","# of Offline Vertices",80, -0.5, 79.5);
249 
251 
252  //ENERGY VARIABLES
253 
254  book3DHistoVector(h_phoE_, "1D","phoE","Energy;E (GeV)",eBin,eMin,eMax);
255  book3DHistoVector(h_phoSigmaEoverE_, "1D","phoSigmaEoverE","#sigma_{E}/E; #sigma_{E}/E", 100,0.,0.08);
256  book3DHistoVector(p_phoSigmaEoverEvsNVtx_, "Profile","phoSigmaEoverEvsNVtx","#sigma_{E}/E vs NVtx; N_{vtx}; #sigma_{E}/E",80, -0.5, 79.5, 100,0., 0.08);
257  book3DHistoVector(h_phoEt_, "1D","phoEt","E_{T};E_{T} (GeV)", etBin,etMin,etMax);
258 
259 
260  //NUMBER OF PHOTONS
261 
262  book3DHistoVector(h_nPho_, "1D","nPho","Number of Photons per Event;# #gamma",numberBin,numberMin,numberMax);
263 
264  //GEOMETRICAL VARIABLES
265 
266  //photon eta/phi
267  book2DHistoVector(h_phoEta_, "1D","phoEta","#eta;#eta",etaBin,etaMin,etaMax) ;
268  book3DHistoVector(h_phoPhi_, "1D","phoPhi","#phi;#phi",phiBin,phiMin,phiMax) ;
269 
270  //supercluster eta/phi
271  book2DHistoVector(h_scEta_, "1D","scEta","SuperCluster #eta;#eta",etaBin,etaMin,etaMax) ;
272  book3DHistoVector(h_scPhi_, "1D","scPhi","SuperCluster #phi;#phi",phiBin,phiMin,phiMax) ;
273 
274  //SHOWER SHAPE VARIABLES
275 
276  //r9
277  book3DHistoVector(h_r9_, "1D","r9","R9;R9",r9Bin,r9Min, r9Max);
278  if (standAlone_) book2DHistoVector(h_r9VsEt_, "2D","r9VsEt2D","R9 vs E_{T};E_{T} (GeV);R9",reducedEtBin,etMin,etMax,reducedR9Bin,r9Min,r9Max);
279  book2DHistoVector(p_r9VsEt_, "Profile","r9VsEt","Avg R9 vs E_{T};E_{T} (GeV);R9",etBin,etMin,etMax,r9Bin,r9Min,r9Max);
280  if (standAlone_) book2DHistoVector(h_r9VsEta_, "2D","r9VsEta2D","R9 vs #eta;#eta;R9",reducedEtaBin,etaMin,etaMax,reducedR9Bin,r9Min,r9Max);
281  book2DHistoVector(p_r9VsEta_, "Profile","r9VsEta","Avg R9 vs #eta;#eta;R9",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
282 
283  //sigma ieta ieta
284  book3DHistoVector(h_phoSigmaIetaIeta_, "1D","phoSigmaIetaIeta","#sigma_{i#etai#eta};#sigma_{i#etai#eta}",sigmaIetaBin,sigmaIetaMin,sigmaIetaMax);
285  if (standAlone_) book2DHistoVector(h_sigmaIetaIetaVsEta_, "2D","sigmaIetaIetaVsEta2D","#sigma_{i#etai#eta} vs #eta;#eta;#sigma_{i#etai#eta}",reducedEtaBin,etaMin,etaMax,sigmaIetaBin,sigmaIetaMin,sigmaIetaMax);
286  book2DHistoVector(p_sigmaIetaIetaVsEta_, "Profile","sigmaIetaIetaVsEta","Avg #sigma_{i#etai#eta} vs #eta;#eta;#sigma_{i#etai#eta}",etaBin,etaMin,etaMax,sigmaIetaBin,sigmaIetaMin,sigmaIetaMax);
287 
288  //e1x5
289  if (standAlone_) book2DHistoVector(h_e1x5VsEt_, "2D","e1x5VsEt2D","E1x5 vs E_{T};E_{T} (GeV);E1X5 (GeV)",reducedEtBin,etMin,etMax,reducedEtBin,etMin,etMax);
290  book2DHistoVector(p_e1x5VsEt_, "Profile","e1x5VsEt","Avg E1x5 vs E_{T};E_{T} (GeV);E1X5 (GeV)",etBin,etMin,etMax,etBin,etMin,etMax);
291  if (standAlone_) book2DHistoVector(h_e1x5VsEta_, "2D","e1x5VsEta2D","E1x5 vs #eta;#eta;E1X5 (GeV)",reducedEtaBin,etaMin,etaMax,reducedEtBin,etMin,etMax);
292  book2DHistoVector(p_e1x5VsEta_, "Profile","e1x5VsEta","Avg E1x5 vs #eta;#eta;E1X5 (GeV)",etaBin,etaMin,etaMax,etBin,etMin,etMax);
293 
294  //e2x5
295  if (standAlone_) book2DHistoVector(h_e2x5VsEt_, "2D","e2x5VsEt2D","E2x5 vs E_{T};E_{T} (GeV);E2X5 (GeV)",reducedEtBin,etMin,etMax,reducedEtBin,etMin,etMax);
296  book2DHistoVector(p_e2x5VsEt_, "Profile","e2x5VsEt","Avg E2x5 vs E_{T};E_{T} (GeV);E2X5 (GeV)",etBin,etMin,etMax,etBin,etMin,etMax);
297  if (standAlone_) book2DHistoVector(h_e2x5VsEta_, "2D","e2x5VsEta2D","E2x5 vs #eta;#eta;E2X5 (GeV)",reducedEtaBin,etaMin,etaMax,reducedEtBin,etMin,etMax);
298  book2DHistoVector(p_e2x5VsEta_, "Profile","e2x5VsEta","Avg E2x5 vs #eta;#eta;E2X5 (GeV)",etaBin,etaMin,etaMax,etBin,etMin,etMax);
299 
300  //r1x5
301  if (standAlone_) book2DHistoVector(h_r1x5VsEt_, "2D","r1x5VsEt2D","R1x5 vs E_{T};E_{T} (GeV);R1X5",reducedEtBin,etMin,etMax,reducedR9Bin,r9Min,r9Max);
302  book2DHistoVector(p_r1x5VsEt_, "Profile","r1x5VsEt","Avg R1x5 vs E_{T};E_{T} (GeV);R1X5",etBin,etMin,etMax,r9Bin,r9Min,r9Max);
303  if (standAlone_) book2DHistoVector(h_r1x5VsEta_, "2D","r1x5VsEta2D","R1x5 vs #eta;#eta;R1X5",reducedEtaBin,etaMin,etaMax,reducedR9Bin,r9Min,r9Max);
304  book2DHistoVector(p_r1x5VsEta_, "Profile","r1x5VsEta","Avg R1x5 vs #eta;#eta;R1X5",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
305 
306  //r2x5
307  if (standAlone_) book2DHistoVector( h_r2x5VsEt_ ,"2D","r2x5VsEt2D","R2x5 vs E_{T};E_{T} (GeV);R2X5",reducedEtBin,etMin,etMax,reducedR9Bin,r9Min,r9Max);
308  book2DHistoVector( p_r2x5VsEt_ ,"Profile","r2x5VsEt","Avg R2x5 vs E_{T};E_{T} (GeV);R2X5",etBin,etMin,etMax,r9Bin,r9Min,r9Max);
309  if (standAlone_) book2DHistoVector( h_r2x5VsEta_ ,"2D","r2x5VsEta2D","R2x5 vs #eta;#eta;R2X5",reducedEtaBin,etaMin,etaMax,reducedR9Bin,r9Min,r9Max);
310  book2DHistoVector( p_r2x5VsEta_ ,"Profile","r2x5VsEta","Avg R2x5 vs #eta;#eta;R2X5",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
311 
312  //maxEXtalOver3x3
313  if (standAlone_) 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);
314  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);
315  if (standAlone_) book2DHistoVector( h_maxEXtalOver3x3VsEta_ ,"2D","maxEXtalOver3x3VsEta2D","(Max Xtal E)/E3x3 vs #eta;#eta;(Max Xtal E)/E3x3",reducedEtaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
316  book2DHistoVector( p_maxEXtalOver3x3VsEta_ ,"Profile","maxEXtalOver3x3VsEta","Avg (Max Xtal E)/E3x3 vs #eta;#eta;(Max Xtal E)/E3x3",etaBin,etaMin,etaMax,r9Bin,r9Min,r9Max);
317 
318  //TRACK ISOLATION VARIABLES
319  //nTrackIsolSolid
320  book2DHistoVector( h_nTrackIsolSolid_ ,"1D","nIsoTracksSolid","Number Of Tracks in the Solid Iso Cone;# tracks",numberBin,numberMin,numberMax);
321  if (standAlone_) 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);
322  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);
323  if (standAlone_) book2DHistoVector( h_nTrackIsolSolidVsEta_ ,"2D","nIsoTracksSolidVsEta2D","Number Of Tracks in the Solid Iso Cone vs #eta;#eta;# tracks",reducedEtaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
324  book2DHistoVector( p_nTrackIsolSolidVsEta_ ,"Profile","nIsoTracksSolidVsEta","Avg Number Of Tracks in the Solid Iso Cone vs #eta;#eta;# tracks",etaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
325 
326  //nTrackIsolHollow
327  book2DHistoVector( h_nTrackIsolHollow_ ,"1D","nIsoTracksHollow","Number Of Tracks in the Hollow Iso Cone;# tracks",numberBin,numberMin,numberMax);
328  if (standAlone_) 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);
329  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);
330  if (standAlone_) book2DHistoVector( h_nTrackIsolHollowVsEta_ ,"2D","nIsoTracksHollowVsEta2D","Number Of Tracks in the Hollow Iso Cone vs #eta;#eta;# tracks",reducedEtaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
331  book2DHistoVector( p_nTrackIsolHollowVsEta_ ,"Profile","nIsoTracksHollowVsEta","Avg Number Of Tracks in the Hollow Iso Cone vs #eta;#eta;# tracks",etaBin,etaMin, etaMax,numberBin,numberMin,numberMax);
332 
333  //trackPtSumSolid
334  book2DHistoVector( h_trackPtSumSolid_ ,"1D","isoPtSumSolid","Track P_{T} Sum in the Solid Iso Cone;P_{T} (GeV)",sumBin,sumMin,sumMax);
335  if (standAlone_) 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);
336  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);
337  if (standAlone_) book2DHistoVector( h_trackPtSumSolidVsEta_ ,"2D","isoPtSumSolidVsEta2D","Track P_{T} Sum in the Solid Iso Cone;#eta;P_{T} (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
338  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);
339 
340  //trackPtSumHollow
341  book2DHistoVector( h_trackPtSumHollow_ ,"1D","isoPtSumHollow","Track P_{T} Sum in the Hollow Iso Cone;P_{T} (GeV)",sumBin,sumMin,sumMax);
342  if (standAlone_) 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);
343  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);
344  if (standAlone_) book2DHistoVector( h_trackPtSumHollowVsEta_ ,"2D","isoPtSumHollowVsEta2D","Track P_{T} Sum in the Hollow Iso Cone;#eta;P_{T} (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
345  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);
346 
347 
348  //CALORIMETER ISOLATION VARIABLES
349 
350  //ecal sum
351  book2DHistoVector( h_ecalSum_ ,"1D","ecalSum","Ecal Sum in the Iso Cone;E (GeV)",sumBin,sumMin,sumMax);
352  book2DHistoVector( h_ecalSumEBarrel_,"1D","ecalSumEBarrel","Ecal Sum in the IsoCone for Barrel;E (GeV)",sumBin,sumMin,sumMax);
353  book2DHistoVector( h_ecalSumEEndcap_,"1D","ecalSumEEndcap","Ecal Sum in the IsoCone for Endcap;E (GeV)",sumBin,sumMin,sumMax);
354  if (standAlone_) book2DHistoVector( h_ecalSumVsEt_ ,"2D","ecalSumVsEt2D","Ecal Sum in the Iso Cone;E_{T} (GeV);E (GeV)",reducedEtBin,etMin, etMax,reducedSumBin,sumMin,sumMax);
355  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);
356  if (standAlone_) book2DHistoVector( h_ecalSumVsEta_ ,"2D","ecalSumVsEta2D","Ecal Sum in the Iso Cone;#eta;E (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
357  book2DHistoVector( p_ecalSumVsEta_ ,"Profile","ecalSumVsEta","Avg Ecal Sum in the Iso Cone vs #eta;#eta;E (GeV)",etaBin,etaMin, etaMax,sumBin,sumMin,sumMax);
358 
359  //hcal sum
360  book2DHistoVector( h_hcalSum_ ,"1D","hcalSum","Hcal Sum in the Iso Cone;E (GeV)",sumBin,sumMin,sumMax);
361  book2DHistoVector( h_hcalSumEBarrel_,"1D","hcalSumEBarrel","Hcal Sum in the IsoCone for Barrel;E (GeV)",sumBin,sumMin,sumMax);
362  book2DHistoVector( h_hcalSumEEndcap_,"1D","hcalSumEEndcap","Hcal Sum in the IsoCone for Endcap;E (GeV)",sumBin,sumMin,sumMax);
363  if (standAlone_) book2DHistoVector( h_hcalSumVsEt_ ,"2D","hcalSumVsEt2D","Hcal Sum in the Iso Cone;E_{T} (GeV);E (GeV)",reducedEtBin,etMin, etMax,reducedSumBin,sumMin,sumMax);
364  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);
365  if (standAlone_) book2DHistoVector( h_hcalSumVsEta_ ,"2D","hcalSumVsEta2D","Hcal Sum in the Iso Cone;#eta;E (GeV)",reducedEtaBin,etaMin, etaMax,reducedSumBin,sumMin,sumMax);
366  book2DHistoVector( p_hcalSumVsEta_ ,"Profile","hcalSumVsEta","Avg Hcal Sum in the Iso Cone vs #eta;#eta;E (GeV)",etaBin,etaMin, etaMax,sumBin,sumMin,sumMax);
367 
368  //h over e
369  book3DHistoVector( h_hOverE_ ,"1D","hOverE","H/E;H/E",hOverEBin,hOverEMin,hOverEMax);
370  book2DHistoVector( p_hOverEVsEt_ ,"Profile","hOverEVsEt","Avg H/E vs Et;E_{T} (GeV);H/E",etBin,etMin,etMax,hOverEBin,hOverEMin,hOverEMax);
371  book2DHistoVector( p_hOverEVsEta_ ,"Profile","hOverEVsEta","Avg H/E vs #eta;#eta;H/E",etaBin,etaMin,etaMax,hOverEBin,hOverEMin,hOverEMax);
372  book3DHistoVector( h_h1OverE_ ,"1D","h1OverE","H/E for Depth 1;H/E",hOverEBin,hOverEMin,hOverEMax);
373  book3DHistoVector( h_h2OverE_ ,"1D","h2OverE","H/E for Depth 2;H/E",hOverEBin,hOverEMin,hOverEMax);
374 
375  // pf isolation
376  book2DHistoVector( h_phoIsoBarrel_,"1D","phoIsoBarrel","PF photon iso Barrel;E (GeV)",reducedEtBin,etMin,25.);
377  book2DHistoVector( h_phoIsoEndcap_,"1D","phoIsoEndcap","PF photon iso Endcap;E (GeV)",reducedEtBin,etMin,25.);
378  book2DHistoVector( h_chHadIsoBarrel_,"1D","chHadIsoBarrel","PF charged Had iso Barrel;E (GeV)",reducedEtBin,etMin,25.);
379  book2DHistoVector( h_chHadIsoEndcap_,"1D","chHadIsoEndcap","PF charged Had iso Endcap;E (GeV)",reducedEtBin,etMin,25.);
380  book2DHistoVector( h_nHadIsoBarrel_,"1D","neutralHadIsoBarrel","PF neutral Had iso Barrel;E (GeV)",reducedEtBin,etMin,25.);
381  book2DHistoVector( h_nHadIsoEndcap_,"1D","neutralHadIsoEndcap","PF neutral Had iso Endcap;E (GeV)",reducedEtBin,etMin,25.);
382 
383 
384 
385  //OTHER VARIABLES
386  //bad channel histograms
387  book2DHistoVector( h_phoEt_BadChannels_ ,"1D","phoEtBadChannels", "Fraction Containing Bad Channels: E_{T};E_{T} (GeV)",etBin,etMin,etMax);
388  book2DHistoVector( h_phoEta_BadChannels_ ,"1D","phoEtaBadChannels","Fraction Containing Bad Channels: #eta;#eta",etaBin,etaMin,etaMax);
389  book2DHistoVector( h_phoPhi_BadChannels_ ,"1D","phoPhiBadChannels","Fraction Containing Bad Channels: #phi;#phi",phiBin,phiMin,phiMax);
390 
391 
393 
394  dbe_->setCurrentFolder("Egamma/"+fName_+"/AllPhotons/Et Above 0 GeV/Conversions");
395 
396  //ENERGY VARIABLES
397 
398  book3DHistoVector( h_phoConvE_ ,"1D","phoConvE","E;E (GeV)",eBin,eMin,eMax);
399  book3DHistoVector( h_phoConvEt_ ,"1D","phoConvEt","E_{T};E_{T} (GeV)",etBin,etMin,etMax);
400 
401  //GEOMETRICAL VARIABLES
402 
403  book2DHistoVector( h_phoConvEta_ ,"1D","phoConvEta","#eta;#eta",etaBin,etaMin,etaMax);
404  book3DHistoVector( h_phoConvPhi_ ,"1D","phoConvPhi","#phi;#phi",phiBin,phiMin,phiMax);
405 
406  //NUMBER OF PHOTONS
407 
408  book3DHistoVector( h_nConv_ ,"1D","nConv","Number Of Conversions per Event ;# conversions",numberBin,numberMin,numberMax);
409 
410  //SHOWER SHAPE VARIABLES
411 
412  book3DHistoVector( h_phoConvR9_ ,"1D","phoConvR9","R9;R9",r9Bin,r9Min,r9Max);
413 
414  //TRACK RELATED VARIABLES
415 
416  book3DHistoVector( h_eOverPTracks_ ,"1D","eOverPTracks","E/P;E/P",eOverPBin,eOverPMin,eOverPMax);
417  book3DHistoVector( h_pOverETracks_ ,"1D","pOverETracks","P/E;P/E",eOverPBin,eOverPMin,eOverPMax);
418 
419  book3DHistoVector( h_dPhiTracksAtVtx_ ,"1D","dPhiTracksAtVtx", "#Delta#phi of Tracks at Vertex;#Delta#phi",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
420  book3DHistoVector( h_dPhiTracksAtEcal_ ,"1D","dPhiTracksAtEcal", "Abs(#Delta#phi) of Tracks at Ecal;#Delta#phi",dPhiTracksBin,0.,dPhiTracksMax);
421  book3DHistoVector( h_dEtaTracksAtEcal_ ,"1D","dEtaTracksAtEcal", "#Delta#eta of Tracks at Ecal;#Delta#eta",dEtaTracksBin,dEtaTracksMin,dEtaTracksMax);
422 
423  book3DHistoVector( h_dCotTracks_ ,"1D","dCotTracks","#Deltacot(#theta) of Tracks;#Deltacot(#theta)",dEtaTracksBin,dEtaTracksMin,dEtaTracksMax);
424  book2DHistoVector( p_dCotTracksVsEta_ ,"Profile","dCotTracksVsEta","Avg #Deltacot(#theta) of Tracks vs #eta;#eta;#Deltacot(#theta)",etaBin,etaMin,etaMax,dEtaTracksBin,dEtaTracksMin,dEtaTracksMax);
425 
426  book2DHistoVector( p_nHitsVsEta_ ,"Profile","nHitsVsEta","Avg Number of Hits per Track vs #eta;#eta;# hits",etaBin,etaMin,etaMax,etaBin,0,16);
427 
428  book2DHistoVector( h_tkChi2_ ,"1D","tkChi2","#chi^{2} of Track Fitting;#chi^{2}",chi2Bin,chi2Min,chi2Max);
429  book2DHistoVector( p_tkChi2VsEta_ ,"Profile","tkChi2VsEta","Avg #chi^{2} of Track Fitting vs #eta;#eta;#chi^{2}",etaBin,etaMin,etaMax,chi2Bin,chi2Min,chi2Max);
430 
431  //VERTEX RELATED VARIABLES
432 
433  book2DHistoVector( h_convVtxRvsZ_ ,"2D","convVtxRvsZ","Vertex Position;Z (cm);R (cm)",500,zMin,zMax,rBin,rMin,rMax);
434  book2DHistoVector( h_convVtxZEndcap_ ,"1D","convVtxZEndcap", "Vertex Position: #eta > 1.5;Z (cm)",zBin,zMin,zMax);
435  book2DHistoVector( h_convVtxZ_ ,"1D","convVtxZ", "Vertex Position;Z (cm)",zBin,zMin,zMax);
436  book2DHistoVector( h_convVtxR_ ,"1D","convVtxR", "Vertex Position: #eta < 1;R (cm)",rBin,rMin,rMax);
437  book2DHistoVector( h_convVtxYvsX_ ,"2D","convVtxYvsX","Vertex Position: #eta < 1;X (cm);Y (cm)",xBin,xMin,xMax,yBin,yMin,yMax);
438 
439 
440 
441  book2DHistoVector( h_vertexChi2Prob_ ,"1D","vertexChi2Prob","#chi^{2} Probability of Vertex Fitting;#chi^{2}",100,0.,1.0);
442 
443 
444  }//end if(dbe_)
445 
446 
447 }//end BeginJob
448 
449 
450 
451 void PhotonAnalyzer::analyze( const edm::Event& e, const edm::EventSetup& esup )
452 {
453  using namespace edm;
454 
455  if (nEvt_% prescaleFactor_ ) return;
456  nEvt_++;
457  LogInfo(fName_) << "PhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
458 
459  // Get the trigger results
460  bool validTriggerEvent=true;
461  edm::Handle<trigger::TriggerEvent> triggerEventHandle;
462  trigger::TriggerEvent triggerEvent;
463  e.getByToken(triggerEvent_token_,triggerEventHandle);
464  if(!triggerEventHandle.isValid()) {
465  edm::LogInfo(fName_) << "Error! Can't get the product: triggerEvent_" << endl;
466  validTriggerEvent=false;
467  }
468  if(validTriggerEvent) triggerEvent = *(triggerEventHandle.product());
469 
470  // Get the reconstructed photons
471  // bool validPhotons=true;
472  Handle<reco::PhotonCollection> photonHandle;
473  reco::PhotonCollection photonCollection;
474  e.getByToken(photon_token_ , photonHandle);
475  if ( !photonHandle.isValid()) {
476  edm::LogInfo(fName_) << "Error! Can't get the product: photon_token_" << endl;
477  // validPhotons=false;
478  }
479  // if(validPhotons) photonCollection = *(photonHandle.product());
480 
481  // Get the PhotonId objects
482  bool validloosePhotonID=true;
483  Handle<edm::ValueMap<bool> > loosePhotonFlag;
484  edm::ValueMap<bool> loosePhotonID;
485  e.getByToken(PhotonIDLoose_token_, loosePhotonFlag);
486  if ( !loosePhotonFlag.isValid()) {
487  edm::LogInfo(fName_) << "Error! Can't get the product: PhotonIDLoose_token_" << endl;
488  validloosePhotonID=false;
489  }
490  if (validloosePhotonID) loosePhotonID = *(loosePhotonFlag.product());
491 
492  bool validtightPhotonID=true;
493  Handle<edm::ValueMap<bool> > tightPhotonFlag;
494  edm::ValueMap<bool> tightPhotonID;
495  e.getByToken(PhotonIDTight_token_, tightPhotonFlag);
496  if ( !tightPhotonFlag.isValid()) {
497  edm::LogInfo(fName_) << "Error! Can't get the product: PhotonIDTight_token_" << endl;
498  validtightPhotonID=false;
499  }
500  if (validtightPhotonID) tightPhotonID = *(tightPhotonFlag.product());
501 
502 
503  edm::Handle<reco::VertexCollection> vtxH;
504  if ( !isHeavyIon_) {
505  e.getByToken(offline_pvToken_, vtxH);
506  h_nRecoVtx_ ->Fill (float(vtxH->size()));
507  }
508 
509  // Create array to hold #photons/event information
510  int nPho[100][3][3];
511 
512  for (int cut=0; cut!=100; ++cut){
513  for (unsigned int type=0; type!=types_.size(); ++type){
514  for (unsigned int part=0; part!=parts_.size(); ++part){
515  nPho[cut][type][part] = 0;
516  }
517  }
518  }
519  // Create array to hold #conversions/event information
520  int nConv[100][3][3];
521 
522  for (int cut=0; cut!=100; ++cut){
523  for (unsigned int type=0; type!=types_.size(); ++type){
524  for (unsigned int part=0; part!=parts_.size(); ++part){
525  nConv[cut][type][part] = 0;
526  }
527  }
528  }
529 
530 
531 
532  //Prepare list of photon-related HLT filter names
533 
534  vector<int> Keys;
535 
536  for(uint filterIndex=0;filterIndex<triggerEvent.sizeFilters();++filterIndex){ //loop over all trigger filters in event (i.e. filters passed)
537 
538  string label = triggerEvent.filterTag(filterIndex).label();
539 
540  if(label.find( "Photon" ) != string::npos ) { //get photon-related filters
541 
542  for(uint filterKeyIndex=0;filterKeyIndex<triggerEvent.filterKeys(filterIndex).size();++filterKeyIndex){ //loop over keys to objects passing this filter
543  Keys.push_back(triggerEvent.filterKeys(filterIndex)[filterKeyIndex]); //add keys to a vector for later reference
544  }
545 
546  }
547 
548  }
549 
550  // sort Keys vector in ascending order
551  // and erases duplicate entries from the vector
552  sort(Keys.begin(),Keys.end());
553  for ( uint i=0 ; i<Keys.size() ; )
554  {
555  if (i!=(Keys.size()-1))
556  {
557  if (Keys[i]==Keys[i+1]) Keys.erase(Keys.begin()+i+1) ;
558  else ++i ;
559  }
560  else ++i ;
561  }
562 
563  //We now have a vector of unique keys to TriggerObjects passing a photon-related filter
564 
565  // old int photonCounter = 0;
566 
568  for(unsigned int iPho=0; iPho < photonHandle->size(); iPho++) {
569  reco::PhotonRef aPho(reco::PhotonRef(photonHandle, iPho));
570  // for( reco::PhotonCollection::const_iterator iPho = photonCollection.begin(); iPho != photonCollection.end(); iPho++) {
571 
572 
573  //for HLT efficiency plots
574 
575  h_phoEta_preHLT_->Fill(aPho->eta());
576  h_phoEt_preHLT_->Fill( aPho->et());
577 
578 
579  double deltaR=1000.;
580  double deltaRMin=1000.;
581  double deltaRMax=0.05;//sets deltaR threshold for matching photons to trigger objects
582 
583 
584  for (vector<int>::const_iterator objectKey=Keys.begin();objectKey!=Keys.end();objectKey++){ //loop over keys to objects that fired photon triggers
585 
586  deltaR = reco::deltaR(triggerEvent.getObjects()[(*objectKey)].eta(),triggerEvent.getObjects()[(*objectKey)].phi(),aPho->superCluster()->eta(),aPho->superCluster()->phi());
587  if(deltaR < deltaRMin) deltaRMin = deltaR;
588 
589  }
590 
591  if(deltaRMin > deltaRMax) { //photon fails delta R cut
592  if(useTriggerFiltering_) continue; //throw away photons that haven't passed any photon filters
593  }
594 
595  if(deltaRMin <= deltaRMax) { //photon passes delta R cut
596  h_phoEta_postHLT_->Fill(aPho->eta() );
597  h_phoEt_postHLT_->Fill( aPho->et() );
598  }
599 
600  // if (aPho->et() < minPhoEtCut_) continue;
601  bool isLoosePhoton(false), isTightPhoton(false);
602  if ( photonSelection (aPho) ) isLoosePhoton=true ;
603 
604 
605  nEntry_++;
606 
607  // old edm::Ref<reco::PhotonCollection> photonref(photonHandle, photonCounter);
608  // old photonCounter++;
609 
610 
611  // old if ( !isHeavyIon_ ) {
612  // isLoosePhoton = (loosePhotonID)[photonref];
613  // isTightPhoton = (tightPhotonID)[photonref];
614  // }
615 
616 
617  //find out which part of the Ecal contains the photon
618  bool phoIsInBarrel=false;
619  bool phoIsInEndcap=false;
620  float etaPho = aPho->superCluster()->eta();
621  if ( fabs(etaPho) < 1.479 )
622  phoIsInBarrel=true;
623  else {
624  phoIsInEndcap=true;
625  }
626 
627  int part = 0;
628  if ( phoIsInBarrel ) part = 1;
629  if ( phoIsInEndcap ) part = 2;
630 
632  bool isIsolated=false;
633  if ( isolationStrength_ == 0) isIsolated = isLoosePhoton;
634  if ( isolationStrength_ == 1) isIsolated = isTightPhoton;
635 
636  int type=0;
637  if ( isIsolated ) type=1;
638  if ( !excludeBkgHistos_ && !isIsolated ) type=2;
639 
640 
641  //get rechit collection containing this photon
642  bool validEcalRecHits=true;
643  edm::Handle<EcalRecHitCollection> ecalRecHitHandle;
644  EcalRecHitCollection ecalRecHitCollection;
645  if ( phoIsInBarrel ) {
646  // Get handle to barrel rec hits
647  e.getByToken(barrelRecHit_token_, ecalRecHitHandle);
648  if (!ecalRecHitHandle.isValid()) {
649  edm::LogError(fName_) << "Error! Can't get the product: barrelRecHit_token_";
650  validEcalRecHits=false;
651  }
652  }
653  else if ( phoIsInEndcap ) {
654  // Get handle to endcap rec hits
655  e.getByToken(endcapRecHit_token_, ecalRecHitHandle);
656  if (!ecalRecHitHandle.isValid()) {
657  edm::LogError(fName_) << "Error! Can't get the product: endcapRecHit_token";
658  validEcalRecHits=false;
659  }
660  }
661  if (validEcalRecHits) ecalRecHitCollection = *(ecalRecHitHandle.product());
662 
663 
664  //if (aPho->isEBEEGap()) continue; //cut out gap photons
665 
666 
667  //filling histograms to make isolation efficiencies
668  if(isLoosePhoton){
669  h_phoEta_Loose_->Fill(aPho->eta());
670  h_phoEt_Loose_->Fill( aPho->et() );
671  }
672  if(isTightPhoton){
673  h_phoEta_Tight_->Fill(aPho->eta());
674  h_phoEt_Tight_->Fill( aPho->et() );
675  }
676 
677 
678 
679  for (int cut = 0; cut !=numberOfSteps_; ++cut) { //loop over different transverse energy cuts
680  double Et = aPho->et();
681  bool passesCuts = false;
682 
683  //sorting the photon into the right Et-dependant folder
684  if ( useBinning_ && Et > (cut+1)*cutStep_ && ( (Et < (cut+2)*cutStep_) | (cut == numberOfSteps_-1) ) ){
685  passesCuts = true;
686  }
687  else if ( !useBinning_ && Et > (cut+1)*cutStep_ ){
688  passesCuts = true;
689  }
690 
691  if (passesCuts){
692 
693  //filling isolation variable histograms
694 
695  //tracker isolation variables
696  fill2DHistoVector(h_nTrackIsolSolid_, aPho->nTrkSolidConeDR04(), cut,type);
697  fill2DHistoVector(h_nTrackIsolHollow_,aPho->nTrkHollowConeDR04(),cut,type);
698 
699  if (standAlone_) fill2DHistoVector(h_nTrackIsolSolidVsEta_, aPho->eta(),aPho->nTrkSolidConeDR04(), cut,type);
700  fill2DHistoVector(p_nTrackIsolSolidVsEta_, aPho->eta(),aPho->nTrkSolidConeDR04(), cut,type);
701  if (standAlone_) fill2DHistoVector(h_nTrackIsolHollowVsEta_,aPho->eta(),aPho->nTrkHollowConeDR04(),cut,type);
702  fill2DHistoVector(p_nTrackIsolHollowVsEta_,aPho->eta(),aPho->nTrkHollowConeDR04(),cut,type);
703 
704  if (standAlone_) fill2DHistoVector(h_nTrackIsolSolidVsEt_, aPho->et(), aPho->nTrkSolidConeDR04(), cut,type);
705  fill2DHistoVector(p_nTrackIsolSolidVsEt_, aPho->et(), aPho->nTrkSolidConeDR04(), cut,type);
706  if (standAlone_) fill2DHistoVector(h_nTrackIsolHollowVsEt_, aPho->et(), aPho->nTrkHollowConeDR04(),cut,type);
707  fill2DHistoVector(p_nTrackIsolHollowVsEt_, aPho->et(), aPho->nTrkHollowConeDR04(),cut,type);
708 
710  fill2DHistoVector(h_trackPtSumSolid_, aPho->trkSumPtSolidConeDR04(),cut,type);
711  fill2DHistoVector(h_trackPtSumHollow_,aPho->trkSumPtSolidConeDR04(),cut,type);
712 
713  if (standAlone_) fill2DHistoVector(h_trackPtSumSolidVsEta_, aPho->eta(),aPho->trkSumPtSolidConeDR04(), cut,type);
714  fill2DHistoVector(p_trackPtSumSolidVsEta_, aPho->eta(),aPho->trkSumPtSolidConeDR04(), cut,type);
715  if (standAlone_) fill2DHistoVector(h_trackPtSumHollowVsEta_,aPho->eta(),aPho->trkSumPtHollowConeDR04(),cut,type);
716  fill2DHistoVector(p_trackPtSumHollowVsEta_,aPho->eta(),aPho->trkSumPtHollowConeDR04(),cut,type);
717 
718  if (standAlone_) fill2DHistoVector(h_trackPtSumSolidVsEt_, aPho->et(), aPho->trkSumPtSolidConeDR04(), cut,type);
719  fill2DHistoVector(p_trackPtSumSolidVsEt_, aPho->et(), aPho->trkSumPtSolidConeDR04(), cut,type);
720  if (standAlone_) fill2DHistoVector(h_trackPtSumHollowVsEt_, aPho->et(), aPho->trkSumPtHollowConeDR04(),cut,type);
721  fill2DHistoVector(p_trackPtSumHollowVsEt_, aPho->et(), aPho->trkSumPtHollowConeDR04(),cut,type);
722  //calorimeter isolation variables
723 
724  fill2DHistoVector(h_ecalSum_,aPho->ecalRecHitSumEtConeDR04(),cut,type);
725  if(aPho->isEB()){fill2DHistoVector(h_ecalSumEBarrel_,aPho->ecalRecHitSumEtConeDR04(),cut,type);}
726  if(aPho->isEE()){fill2DHistoVector(h_ecalSumEEndcap_,aPho->ecalRecHitSumEtConeDR04(),cut,type);}
727  if (standAlone_) fill2DHistoVector(h_ecalSumVsEta_,aPho->eta(),aPho->ecalRecHitSumEtConeDR04(),cut,type);
728  fill2DHistoVector(p_ecalSumVsEta_,aPho->eta(),aPho->ecalRecHitSumEtConeDR04(),cut,type);
729  if (standAlone_) fill2DHistoVector(h_ecalSumVsEt_, aPho->et(), aPho->ecalRecHitSumEtConeDR04(),cut,type);
730  fill3DHistoVector(p_ecalSumVsEt_, aPho->et(), aPho->ecalRecHitSumEtConeDR04(),cut,type,part);
731 
733 
734  fill2DHistoVector(h_hcalSum_,aPho->hcalTowerSumEtConeDR04(),cut,type);
735  if(aPho->isEB()){fill2DHistoVector(h_hcalSumEBarrel_,aPho->hcalTowerSumEtConeDR04(),cut,type);}
736  if(aPho->isEE()){fill2DHistoVector(h_hcalSumEEndcap_,aPho->hcalTowerSumEtConeDR04(),cut,type);}
737  if (standAlone_) fill2DHistoVector(h_hcalSumVsEta_,aPho->eta(),aPho->hcalTowerSumEtConeDR04(),cut,type);
738  fill2DHistoVector(p_hcalSumVsEta_,aPho->eta(),aPho->hcalTowerSumEtConeDR04(),cut,type);
739  if (standAlone_) fill2DHistoVector(h_hcalSumVsEt_, aPho->et(), aPho->hcalTowerSumEtConeDR04(),cut,type);
740  fill3DHistoVector(p_hcalSumVsEt_, aPho->et(), aPho->hcalTowerSumEtConeDR04(),cut,type,part);
741 
742  fill3DHistoVector(h_hOverE_,aPho->hadronicOverEm(),cut,type,part);
743  fill2DHistoVector(p_hOverEVsEta_,aPho->eta(),aPho->hadronicOverEm(),cut,type);
744  fill2DHistoVector(p_hOverEVsEt_, aPho->et(), aPho->hadronicOverEm(),cut,type);
745 
746  fill3DHistoVector(h_h1OverE_,aPho->hadronicDepth1OverEm(),cut,type,part);
747  fill3DHistoVector(h_h2OverE_,aPho->hadronicDepth2OverEm(),cut,type,part);
748 
749 
750  // filling pf isolation variables
751  if(aPho->isEB()) {
752  fill2DHistoVector( h_phoIsoBarrel_, aPho->photonIso(),cut,type);
753  fill2DHistoVector( h_chHadIsoBarrel_, aPho->chargedHadronIso(),cut,type);
754  fill2DHistoVector( h_nHadIsoBarrel_, aPho->neutralHadronIso(),cut,type);
755  }
756  if(aPho->isEE()) {
757  fill2DHistoVector( h_phoIsoEndcap_, aPho->photonIso(),cut,type);
758  fill2DHistoVector( h_chHadIsoEndcap_, aPho->chargedHadronIso(),cut,type);
759  fill2DHistoVector( h_nHadIsoEndcap_, aPho->neutralHadronIso(),cut,type);
760  }
761 
762 
763  //filling photon histograms
764 
765  nPho[cut][0][0]++;
766  nPho[cut][0][part]++;
767  nPho[cut][type][0]++;
768  nPho[cut][type][part]++;
769 
770  //energy variables
771 
772  fill3DHistoVector(h_phoE_, aPho->energy(),cut,type,part);
773  fill3DHistoVector(h_phoSigmaEoverE_, aPho->getCorrectedEnergyError(aPho->getCandidateP4type())/aPho->energy(),cut,type,part);
774 
775  if ( !isHeavyIon_) fill3DHistoVector(p_phoSigmaEoverEvsNVtx_, float(vtxH->size()), aPho->getCorrectedEnergyError(aPho->getCandidateP4type())/aPho->energy(),cut,type,part);
776 
777  fill3DHistoVector(h_phoEt_,aPho->et(), cut,type,part);
778 
779  //geometrical variables
780 
781  fill2DHistoVector(h_phoEta_,aPho->eta(),cut,type);
782  fill2DHistoVector(h_scEta_, aPho->superCluster()->eta(),cut,type);
783 
784  fill3DHistoVector(h_phoPhi_,aPho->phi(),cut,type,part);
785  fill3DHistoVector(h_scPhi_, aPho->superCluster()->phi(),cut,type,part);
786 
787  //shower shape variables
788 
789  fill3DHistoVector(h_r9_,aPho->r9(),cut,type,part);
790  if (standAlone_) fill2DHistoVector(h_r9VsEta_,aPho->eta(),aPho->r9(),cut,type);
791  fill2DHistoVector(p_r9VsEta_,aPho->eta(),aPho->r9(),cut,type);
792  if (standAlone_) fill2DHistoVector(h_r9VsEt_, aPho->et(), aPho->r9(),cut,type);
793  fill2DHistoVector(p_r9VsEt_, aPho->et(), aPho->r9(),cut,type);
794 
795  if (standAlone_) fill2DHistoVector(h_e1x5VsEta_,aPho->eta(),aPho->e1x5(),cut,type);
796  fill2DHistoVector(p_e1x5VsEta_,aPho->eta(),aPho->e1x5(),cut,type);
797  if (standAlone_) fill2DHistoVector(h_e1x5VsEt_, aPho->et(), aPho->e1x5(),cut,type);
798  fill2DHistoVector(p_e1x5VsEt_, aPho->et(), aPho->e1x5(),cut,type);
799 
800  if (standAlone_) fill2DHistoVector(h_e2x5VsEta_,aPho->eta(),aPho->e2x5(),cut,type);
801  fill2DHistoVector(p_e2x5VsEta_,aPho->eta(),aPho->e2x5(),cut,type);
802  if (standAlone_) fill2DHistoVector(h_e2x5VsEt_, aPho->et(), aPho->e2x5(),cut,type);
803  fill2DHistoVector(p_e2x5VsEt_, aPho->et(), aPho->e2x5(),cut,type);
804 
805  if (standAlone_) fill2DHistoVector(h_maxEXtalOver3x3VsEta_,aPho->eta(),aPho->maxEnergyXtal()/aPho->e3x3(),cut,type);
806  fill2DHistoVector(p_maxEXtalOver3x3VsEta_,aPho->eta(),aPho->maxEnergyXtal()/aPho->e3x3(),cut,type);
807  if (standAlone_) fill2DHistoVector(h_maxEXtalOver3x3VsEt_, aPho->et(), aPho->maxEnergyXtal()/aPho->e3x3(),cut,type);
808  fill2DHistoVector(p_maxEXtalOver3x3VsEt_, aPho->et(), aPho->maxEnergyXtal()/aPho->e3x3(),cut,type);
809 
810 
811  if (standAlone_) fill2DHistoVector(h_r1x5VsEta_,aPho->eta(),aPho->r1x5(),cut,type);
812  fill2DHistoVector(p_r1x5VsEta_,aPho->eta(),aPho->r1x5(),cut,type);
813  if (standAlone_) fill2DHistoVector(h_r1x5VsEt_, aPho->et(), aPho->r1x5(),cut,type);
814  fill2DHistoVector(p_r1x5VsEt_, aPho->et(), aPho->r1x5(),cut,type);
815 
816  if (standAlone_) fill2DHistoVector(h_r2x5VsEta_,aPho->eta(),aPho->r2x5(),cut,type);
817  fill2DHistoVector(p_r2x5VsEta_,aPho->eta(),aPho->r2x5(),cut,type);
818  if (standAlone_) fill2DHistoVector(h_r2x5VsEt_, aPho->et(), aPho->r2x5(),cut,type);
819  fill2DHistoVector(p_r2x5VsEt_, aPho->et(), aPho->r2x5(),cut,type);
820 
821  fill3DHistoVector(h_phoSigmaIetaIeta_,aPho->sigmaIetaIeta(),cut,type,part);
822  if (standAlone_) fill2DHistoVector(h_sigmaIetaIetaVsEta_,aPho->eta(),aPho->sigmaIetaIeta(),cut,type);
823  fill2DHistoVector(p_sigmaIetaIetaVsEta_,aPho->eta(),aPho->sigmaIetaIeta(),cut,type);
824 
825 
826 
827  //filling histograms for photons containing a bad ECAL channel
828  bool atLeastOneDeadChannel=false;
829  for(reco::CaloCluster_iterator bcIt = aPho->superCluster()->clustersBegin();bcIt != aPho->superCluster()->clustersEnd(); ++bcIt) { //loop over basic clusters in SC
830  for(vector< pair<DetId, float> >::const_iterator rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) { //loop over rec hits in basic cluster
831 
832  for(EcalRecHitCollection::const_iterator it = ecalRecHitCollection.begin(); it != ecalRecHitCollection.end(); ++it) { //loop over all rec hits to find the right ones
833  if (rhIt->first == (*it).id() ) { //found the matching rechit
834  if ( (*it).recoFlag() == 9 ) { //has a bad channel
835  atLeastOneDeadChannel=true;
836  break;
837  }
838  }
839  }
840  }
841  }
842  if ( atLeastOneDeadChannel ) {
843  fill2DHistoVector(h_phoPhi_BadChannels_,aPho->phi(),cut,type);
844  fill2DHistoVector(h_phoEta_BadChannels_,aPho->eta(),cut,type);
845  fill2DHistoVector(h_phoEt_BadChannels_, aPho->et(), cut,type);
846  }
847 
848 
849  // filling conversion-related histograms
850  if(aPho->hasConversionTracks()){
851  nConv[cut][0][0]++;
852  nConv[cut][0][part]++;
853  nConv[cut][type][0]++;
854  nConv[cut][type][part]++;
855  }
856 
857  //loop over conversions (don't forget, we're still inside the photon loop,
858  // i.e. these are all the conversions for this ONE photon, not for all the photons in the event)
859  reco::ConversionRefVector conversions = aPho->conversions();
860  for (unsigned int iConv=0; iConv<conversions.size(); iConv++) {
861 
862  reco::ConversionRef aConv=conversions[iConv];
863 
864  if ( aConv->nTracks() <2 ) continue;
865 
866  //fill histogram for denominator of vertex reconstruction efficiency plot
867  if(cut==0) h_phoEta_Vertex_->Fill(aConv->refittedPairMomentum().eta());
868 
869  if ( !(aConv->conversionVertex().isValid()) ) continue;
870 
871  float chi2Prob = ChiSquaredProbability( aConv->conversionVertex().chi2(), aConv->conversionVertex().ndof() );
872 
873  if(chi2Prob<0.0005) continue;
874 
875  fill2DHistoVector(h_vertexChi2Prob_,chi2Prob,cut,type);
876 
877 
878 
879  fill3DHistoVector(h_phoConvE_, aPho->energy(),cut,type,part);
880  fill3DHistoVector(h_phoConvEt_,aPho->et(),cut,type,part);
881  fill3DHistoVector(h_phoConvR9_,aPho->r9(),cut,type,part);
882 
883  if (cut==0 && isLoosePhoton){
884  h_convEta_Loose_->Fill(aPho->eta());
885  h_convEt_Loose_->Fill( aPho->et() );
886  }
887  if (cut==0 && isTightPhoton){
888  h_convEta_Tight_->Fill(aPho->eta());
889  h_convEt_Tight_->Fill( aPho->et() );
890  }
891 
892  fill2DHistoVector(h_phoConvEta_,aConv->refittedPairMomentum().eta(),cut,type);
893  fill3DHistoVector(h_phoConvPhi_,aConv->refittedPairMomentum().phi(),cut,type,part);
894 
895 
896  //we use the photon position because we'll be dividing it by a photon histogram (not a conversion histogram)
897  fill2DHistoVector(h_phoConvEtaForEfficiency_,aPho->eta(),cut,type);
898  fill3DHistoVector(h_phoConvPhiForEfficiency_,aPho->phi(),cut,type,part);
899 
900 
901  //vertex histograms
902  double convR= sqrt(aConv->conversionVertex().position().perp2());
903  double scalar = aConv->conversionVertex().position().x()*aConv->refittedPairMomentum().x() + aConv->conversionVertex().position().y()*aConv->refittedPairMomentum().y();
904  if ( scalar < 0 ) convR= -convR;
905 
906  fill2DHistoVector(h_convVtxRvsZ_,aConv->conversionVertex().position().z(), convR,cut,type);//trying to "see" R-Z view of tracker
907  fill2DHistoVector(h_convVtxZ_,aConv->conversionVertex().position().z(), cut,type);
908 
909 
910  if( fabs(aPho->eta()) > 1.5){//trying to "see" tracker endcaps
911  fill2DHistoVector(h_convVtxZEndcap_,aConv->conversionVertex().position().z(), cut,type);
912  }
913  else if(fabs(aPho->eta()) < 1){//trying to "see" tracker barrel
914  fill2DHistoVector(h_convVtxR_,convR,cut,type);
915  fill2DHistoVector(h_convVtxYvsX_,aConv->conversionVertex().position().x(),aConv->conversionVertex().position().y(),cut,type);
916  }
917 
918 
919  const std::vector<edm::RefToBase<reco::Track> > tracks = aConv->tracks();
920 
921 
922  for (unsigned int i=0; i<tracks.size(); i++) {
923  fill2DHistoVector(h_tkChi2_,tracks[i]->normalizedChi2(),cut,type);
924  fill2DHistoVector(p_tkChi2VsEta_,aPho->eta(),tracks[i]->normalizedChi2(),cut,type);
925  fill2DHistoVector(p_dCotTracksVsEta_,aPho->eta(),aConv->pairCotThetaSeparation(),cut,type);
926  fill2DHistoVector(p_nHitsVsEta_,aPho->eta(),float(tracks[i]->numberOfValidHits()),cut,type);
927  }
928 
929  //calculating delta eta and delta phi of the two tracks
930 
931  float DPhiTracksAtVtx = -99;
932  float dPhiTracksAtEcal= -99;
933  float dEtaTracksAtEcal= -99;
934 
935  float phiTk1= aConv->tracksPin()[0].phi();
936  float phiTk2= aConv->tracksPin()[1].phi();
937  DPhiTracksAtVtx = phiTk1-phiTk2;
938  DPhiTracksAtVtx = phiNormalization( DPhiTracksAtVtx );
939 
940  if (aConv->bcMatchingWithTracks().size() > 0 && aConv->bcMatchingWithTracks()[0].isNonnull() && aConv->bcMatchingWithTracks()[1].isNonnull() ) {
941  float recoPhi1 = aConv->ecalImpactPosition()[0].phi();
942  float recoPhi2 = aConv->ecalImpactPosition()[1].phi();
943  float recoEta1 = aConv->ecalImpactPosition()[0].eta();
944  float recoEta2 = aConv->ecalImpactPosition()[1].eta();
945 
946  recoPhi1 = phiNormalization(recoPhi1);
947  recoPhi2 = phiNormalization(recoPhi2);
948 
949  dPhiTracksAtEcal = recoPhi1 -recoPhi2;
950  dPhiTracksAtEcal = phiNormalization( dPhiTracksAtEcal );
951  dEtaTracksAtEcal = recoEta1 -recoEta2;
952 
953  }
954 
955 
956  fill3DHistoVector(h_dPhiTracksAtVtx_,DPhiTracksAtVtx,cut,type,part);
957  fill3DHistoVector(h_dPhiTracksAtEcal_,fabs(dPhiTracksAtEcal),cut,type,part);
958  fill3DHistoVector(h_dEtaTracksAtEcal_, dEtaTracksAtEcal,cut,type,part);
959  fill3DHistoVector(h_eOverPTracks_,aConv->EoverPrefittedTracks(),cut,type,part);
960  fill3DHistoVector(h_pOverETracks_,1./aConv->EoverPrefittedTracks(),cut,type,part);
961  fill3DHistoVector(h_dCotTracks_,aConv->pairCotThetaSeparation(),cut,type,part);
962 
963  }//end loop over conversions
964 
965  }//end loop over photons passing cuts
966  }//end loop over transverse energy cuts
967 
968 
969 
970 
971 
972  //make invariant mass plots
973 
974  if (isIsolated && aPho->et()>=invMassEtCut_){
975  for(unsigned int iPho2=iPho+1; iPho2 < photonHandle->size(); iPho2++) {
976  reco::PhotonRef aPho2(reco::PhotonRef(photonHandle, iPho2));
977 
978  // for (reco::PhotonCollection::const_iterator iPho2=iPho+1; iPho2!=photonCollection.end(); iPho2++){
979 
980  // edm::Ref<reco::PhotonCollection> photonref2(photonHandle, photonCounter); //note: it's correct to use photonCounter and not photonCounter+1
981  //since it has already been incremented earlier
982 
983  bool isTightPhoton2(false), isLoosePhoton2(false);
984  if ( photonSelection (aPho2) ) isLoosePhoton2=true;
985 
986  // Old if ( !isHeavyIon_ ) {
987  // isTightPhoton2 = (tightPhotonID)[aPho2];
988  // isLoosePhoton2 = (loosePhotonID)[aPho2];
989  // }
990 
991  bool isIsolated2=false;
992  if ( isolationStrength_ == 0) isIsolated2 = isLoosePhoton2;
993  if ( isolationStrength_ == 1) isIsolated2 = isTightPhoton2;
994 
995  reco::ConversionRefVector conversions = aPho->conversions();
996  reco::ConversionRefVector conversions2 = aPho2->conversions();
997 
998  if(isIsolated2 && aPho2->et()>=invMassEtCut_){
999 
1000  math::XYZTLorentzVector p12 = aPho->p4()+aPho2->p4();
1001  float gamgamMass2 = p12.Dot(p12);
1002 
1003 
1004  h_invMassAllPhotons_ -> Fill(sqrt( gamgamMass2 ));
1005  if(aPho->isEB() && aPho2->isEB()){h_invMassPhotonsEBarrel_ -> Fill(sqrt( gamgamMass2 ));}
1006  if(aPho->isEE() || aPho2->isEE()){h_invMassPhotonsEEndcap_ -> Fill(sqrt( gamgamMass2 ));}
1007 
1008  if(conversions.size()!=0 && conversions[0]->nTracks() >= 2){
1009  if(conversions2.size()!=0 && conversions2[0]->nTracks() >= 2) h_invMassTwoWithTracks_ -> Fill(sqrt( gamgamMass2 ));
1010  else h_invMassOneWithTracks_ -> Fill(sqrt( gamgamMass2 ));
1011  }
1012  else if(conversions2.size()!=0 && conversions2[0]->nTracks() >= 2) h_invMassOneWithTracks_ -> Fill(sqrt( gamgamMass2 ));
1013  else h_invMassZeroWithTracks_ -> Fill(sqrt( gamgamMass2 ));
1014  }
1015 
1016  }
1017 
1018  }
1019 
1020 
1021 
1022  }
1023 
1024 
1025  //filling number of photons/conversions per event histograms
1026  for (int cut = 0; cut !=numberOfSteps_; ++cut) {
1027  for(uint type=0;type!=types_.size();++type){
1028  for(uint part=0;part!=parts_.size();++part){
1029  h_nPho_[cut][type][part]-> Fill (float(nPho[cut][type][part]));
1030  h_nConv_[cut][type][part]-> Fill (float(nConv[cut][type][part]));
1031  }
1032  }
1033  }
1034 
1035 }//End of Analyze method
1036 
1037 void PhotonAnalyzer::endRun(const edm::Run& run, const edm::EventSetup& setup)
1038 {
1039  if(!standAlone_){
1040 
1041 
1042  dbe_->setCurrentFolder("Egamma/"+fName_+"/");
1043  //keep track of how many histos are in each folder
1044  totalNumberOfHistos_efficiencyFolder->Fill(histo_index_efficiency_);
1045  totalNumberOfHistos_invMassFolder->Fill(histo_index_invMass_);
1046  totalNumberOfHistos_photonsFolder->Fill(histo_index_photons_);
1047  totalNumberOfHistos_conversionsFolder->Fill(histo_index_conversions_);
1048 
1049  }
1050 
1051 }
1052 
1053 
1054 void PhotonAnalyzer::endJob()
1055 {
1056  //dbe_->showDirStructure();
1057  if(standAlone_){
1058  dbe_->setCurrentFolder("Egamma/"+fName_+"/");
1059 
1060  //keep track of how many histos are in each folder
1061  totalNumberOfHistos_efficiencyFolder->Fill(histo_index_efficiency_);
1062  totalNumberOfHistos_invMassFolder->Fill(histo_index_invMass_);
1063  totalNumberOfHistos_photonsFolder->Fill(histo_index_photons_);
1064  totalNumberOfHistos_conversionsFolder->Fill(histo_index_conversions_);
1065 
1066 
1067  dbe_->save(outputFileName_);
1068  }
1069 
1070 
1071 }
1072 
1074 
1075 
1076 
1077 float PhotonAnalyzer::phiNormalization(float & phi)
1078 {
1079  const float PI = 3.1415927;
1080  const float TWOPI = 2.0*PI;
1081 
1082  if(phi > PI) {phi = phi - TWOPI;}
1083  if(phi < -PI) {phi = phi + TWOPI;}
1084 
1085  return phi;
1086 }
1087 
1088 
1089 void PhotonAnalyzer::fill2DHistoVector(vector<vector<MonitorElement*> >& histoVector,double x, double y, int cut, int type){
1090 
1091  histoVector[cut][0]->Fill(x,y);
1092  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
1093 
1094 }
1095 
1096 void PhotonAnalyzer::fill2DHistoVector(vector<vector<MonitorElement*> >& histoVector, double x, int cut, int type){
1097 
1098  histoVector[cut][0]->Fill(x);
1099  histoVector[cut][type]->Fill(x);
1100 
1101 }
1102 
1103 void PhotonAnalyzer::fill3DHistoVector(vector<vector<vector<MonitorElement*> > >& histoVector,double x, int cut, int type, int part){
1104 
1105  histoVector[cut][0][0]->Fill(x);
1106  histoVector[cut][0][part]->Fill(x);
1107  histoVector[cut][type][0]->Fill(x);
1108  histoVector[cut][type][part]->Fill(x);
1109 
1110 }
1111 
1112 void PhotonAnalyzer::fill3DHistoVector(vector<vector<vector<MonitorElement*> > >& histoVector,double x, double y, int cut, int type, int part){
1113 
1114  histoVector[cut][0][0]->Fill(x,y);
1115  histoVector[cut][0][part]->Fill(x,y);
1116  histoVector[cut][type][0]->Fill(x,y);
1117  histoVector[cut][type][part]->Fill(x,y);
1118 
1119 }
1120 
1121 
1122 MonitorElement* PhotonAnalyzer::bookHisto(string histoName, string title, int bin, double min, double max)
1123 {
1124 
1125  int histo_index = 0;
1126  stringstream histo_number_stream;
1127 
1128  //determining which folder we're in
1129  if(dbe_->pwd().find( "InvMass" ) != string::npos){
1130  histo_index_invMass_++;
1131  histo_index = histo_index_invMass_;
1132  }
1133  if(dbe_->pwd().find( "Efficiencies" ) != string::npos){
1134  histo_index_efficiency_++;
1135  histo_index = histo_index_efficiency_;
1136  }
1137 
1138  histo_number_stream << "h_";
1139  if(histo_index<10) histo_number_stream << "0";
1140  histo_number_stream << histo_index;
1141 
1142  return dbe_->book1D(histo_number_stream.str()+"_"+histoName,title,bin,min,max);
1143 
1144 }
1145 
1146 
1147 void PhotonAnalyzer::book2DHistoVector(vector<vector<MonitorElement*> > &temp2DVector,
1148  string histoType, string histoName, string title,
1149  int xbin, double xmin,double xmax,
1150  int ybin, double ymin, double ymax)
1151 {
1152  int histo_index = 0;
1153 
1154  vector<MonitorElement*> temp1DVector;
1155 // vector<vector<MonitorElement*> > temp2DVector;
1156 
1157  //determining which folder we're in
1158  bool conversionPlot = false;
1159  if(dbe_->pwd().find( "Conversions" ) != string::npos) conversionPlot = true;
1160  bool TwoDPlot = false;
1161  if(histoName.find( "2D" ) != string::npos) TwoDPlot = true;
1162 
1163  if(conversionPlot){
1164  histo_index_conversions_++;
1165  histo_index = histo_index_conversions_;
1166  }
1167  else{
1168  histo_index_photons_++;
1169  histo_index = histo_index_photons_;
1170  }
1171 
1172  stringstream histo_number_stream;
1173  histo_number_stream << "h_";
1174  if(histo_index<10) histo_number_stream << "0";
1175  histo_number_stream << histo_index << "_";
1176 
1177 
1178 
1179  for(int cut = 0; cut != numberOfSteps_; ++cut){ //looping over Et cut values
1180 
1181  for(uint type=0;type!=types_.size();++type){ //looping over isolation type
1182 
1183  currentFolder_.str("");
1184  currentFolder_ << "Egamma/"+fName_+"/" << types_[type] << "Photons/Et above " << (cut+1)*cutStep_ << " GeV";
1185  if(conversionPlot) currentFolder_ << "/Conversions";
1186 
1187  dbe_->setCurrentFolder(currentFolder_.str());
1188 
1189  string kind;
1190  if(conversionPlot) kind = " Conversions: ";
1191  else kind = " Photons: ";
1192 
1193  if(histoType=="1D") temp1DVector.push_back(dbe_->book1D( histo_number_stream.str()+histoName,types_[type]+kind+title,xbin,xmin,xmax));
1194  else if(histoType=="2D"){
1195  if((TwoDPlot && type==0) || !TwoDPlot){//only book the 2D plots in the "AllPhotons" folder
1196  temp1DVector.push_back(dbe_->book2D( histo_number_stream.str()+histoName,types_[type]+kind+title,xbin,xmin,xmax,ybin,ymin,ymax));
1197  }
1198  }
1199  else if(histoType=="Profile") temp1DVector.push_back(dbe_->bookProfile( histo_number_stream.str()+histoName,types_[type]+kind+title,xbin,xmin,xmax,ybin,ymin,ymax,""));
1200  else cout << "bad histoType\n";
1201  }
1202 
1203  temp2DVector.push_back(temp1DVector);
1204  temp1DVector.clear();
1205  }
1206 
1207 // return temp2DVector;
1208 
1209 }
1210 
1211 
1212 void PhotonAnalyzer::book3DHistoVector(vector<vector<vector<MonitorElement*> > > &temp3DVector,
1213  string histoType, string histoName, string title,
1214  int xbin, double xmin,double xmax,
1215  int ybin, double ymin, double ymax)
1216 {
1217 
1218 
1219  int histo_index = 0;
1220 
1221  vector<MonitorElement*> temp1DVector;
1222  vector<vector<MonitorElement*> > temp2DVector;
1223 // vector<vector<vector<MonitorElement*> > > temp3DVector;
1224 
1225 
1226  //determining which folder we're in
1227  bool conversionPlot = false;
1228  if(dbe_->pwd().find( "Conversions" ) != string::npos) conversionPlot = true;
1229 
1230 
1231  if(conversionPlot){
1232  histo_index_conversions_++;
1233  histo_index = histo_index_conversions_;
1234  }
1235  else{
1236  histo_index_photons_++;
1237  histo_index = histo_index_photons_;
1238  }
1239 
1240 
1241 
1242  stringstream histo_number_stream;
1243  histo_number_stream << "h_";
1244  if(histo_index<10) histo_number_stream << "0";
1245  histo_number_stream << histo_index << "_";
1246 
1247  for(int cut = 0; cut != numberOfSteps_; ++cut){ //looping over Et cut values
1248 
1249  for(uint type=0;type!=types_.size();++type){ //looping over isolation type
1250 
1251  for(uint part=0;part!=parts_.size();++part){ //looping over different parts of the ecal
1252 
1253  currentFolder_.str("");
1254  currentFolder_ << "Egamma/"+fName_+"/" << types_[type] << "Photons/Et above " << (cut+1)*cutStep_ << " GeV";
1255  if(conversionPlot) currentFolder_ << "/Conversions";
1256 
1257  dbe_->setCurrentFolder(currentFolder_.str());
1258 
1259  string kind;
1260  if(conversionPlot) kind = " Conversions: ";
1261  else kind = " Photons: ";
1262 
1263  if(histoType=="1D") temp1DVector.push_back(dbe_->book1D( histo_number_stream.str()+histoName+parts_[part],types_[type]+kind+parts_[part]+": "+title,xbin,xmin,xmax));
1264  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));
1265  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,""));
1266  else cout << "bad histoType\n";
1267 
1268 
1269  }
1270 
1271  temp2DVector.push_back(temp1DVector);
1272  temp1DVector.clear();
1273  }
1274 
1275  temp3DVector.push_back(temp2DVector);
1276  temp2DVector.clear();
1277  }
1278 
1279  // return temp3DVector;
1280 }
1281 
1282 bool PhotonAnalyzer::photonSelection ( const reco::PhotonRef & pho) {
1283 
1284 
1285  bool result=true;
1286  if ( pho->pt() < minPhoEtCut_ ) result=false;
1287  if ( fabs(pho->eta()) > photonMaxEta_ ) result=false;
1288  if ( pho->isEBEEGap() ) result=false;
1289 
1290  double EtCorrHcalIso = pho->hcalTowerSumEtConeDR03() - 0.005*pho->pt();
1291  double EtCorrTrkIso = pho->trkSumPtHollowConeDR03() - 0.002*pho->pt();
1292 
1293  if (pho->r9() <=0.9) {
1294  if (pho->isEB() && (pho->hadTowOverEm()>0.075 || pho->sigmaIetaIeta() > 0.014)) result=false;
1295  if (pho->isEE() && (pho->hadTowOverEm()>0.075 || pho->sigmaIetaIeta() > 0.034)) result=false;
1297  if (EtCorrHcalIso>4.0) result=false;
1298  if (EtCorrTrkIso>4.0) result=false ;
1299  if ( pho->chargedHadronIso() > 4 ) result=false;
1300 
1301  } else {
1302  if (pho->isEB() && (pho->hadTowOverEm()>0.082 || pho->sigmaIetaIeta() > 0.014)) result=false;
1303  if (pho->isEE() && (pho->hadTowOverEm()>0.075 || pho->sigmaIetaIeta() > 0.034)) result=false;
1305  if (EtCorrHcalIso>50.0) result=false;
1306  if (EtCorrTrkIso>50.0) result=false;
1307  if ( pho->chargedHadronIso() > 4 ) result=false;
1308 
1309  }
1310 
1311 
1312 
1313  return result;
1314 }
PhotonAnalyzer::endJob
virtual void endJob()
Definition: PhotonAnalyzer.cc:1054
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
PhotonAnalyzer::beginJob
virtual void beginJob()
Definition: PhotonAnalyzer.cc:75
trigger::TriggerEvent
The single EDProduct to be saved for each event (AOD case)
Definition: TriggerEvent.h:25
DQMStore::book1D
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:954
PI
#define PI
Definition: HcalBeamMonitor.cc:14
PhotonAnalyzer::bookHisto
MonitorElement * bookHisto(std::string histoName, std::string title, int bin, double min, double max)
Definition: PhotonAnalyzer.cc:1122
deltaRMax
double deltaRMax
Definition: PFJetBenchmarkAnalyzer.cc:75
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
cuy.yBin
tuple yBin
Definition: cuy.py:891
edm::SortedCollection< EcalRecHit >::const_iterator
std::vector< EcalRecHit >::const_iterator const_iterator
Definition: SortedCollection.h:82
TWOPI
#define TWOPI
EgammaCoreTools.
Definition: HLTAlCaMonPi0.cc:42
PhotonAnalyzer::fill3DHistoVector
void fill3DHistoVector(std::vector< std::vector< std::vector< MonitorElement * > > > &histoVector, double x, int cut, int type, int part)
Definition: PhotonAnalyzer.cc:1103
SiStripMonitorClusterAlca_cfi.ymax
tuple ymax
Definition: SiStripMonitorClusterAlca_cfi.py:37
SiStripMonitorClusterAlca_cfi.xmin
tuple xmin
Definition: SiStripMonitorClusterAlca_cfi.py:34
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
PhotonAnalyzer::photonSelection
bool photonSelection(const reco::PhotonRef &p)
Definition: PhotonAnalyzer.cc:1282
PhotonAnalyzer::fill2DHistoVector
void fill2DHistoVector(std::vector< std::vector< MonitorElement * > > &histoVector, double x, int cut, int type)
Definition: PhotonAnalyzer.cc:1096
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
conversions_cfi.conversions
tuple conversions
Definition: conversions_cfi.py:8
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:611
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
query.result
tuple result
Definition: query.py:137
reco::deltaR
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
ChiSquaredProbability
float ChiSquaredProbability(double chiSquared, double nrDOF)
Definition: ChiSquaredProbability.cc:13
DQMStore::bookProfile
MonitorElement * bookProfile(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, const char *option="s")
Definition: DQMStore.cc:1268
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:76
SiStripMonitorClusterAlca_cfi.ymin
tuple ymin
Definition: SiStripMonitorClusterAlca_cfi.py:36
PhotonAnalyzer::phiNormalization
float phiNormalization(float &a)
Definition: PhotonAnalyzer.cc:1077
PhotonAnalyzer::book3DHistoVector
void book3DHistoVector(std::vector< std::vector< std::vector< MonitorElement * > > > &toFill, std::string histoType, std::string histoName, std::string title, int xbin, double xmin, double xmax, int ybin=1, double ymin=1, double ymax=2)
Definition: PhotonAnalyzer.cc:1212
dbe_
DQMStore * dbe_
Definition: PFJetBenchmarkAnalyzer.cc:79
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:310
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
PhotonAnalyzer::~PhotonAnalyzer
virtual ~PhotonAnalyzer()
Definition: PhotonAnalyzer.cc:72
PhotonAnalyzer::endRun
virtual void endRun(const edm::Run &, const edm::EventSetup &)
Definition: PhotonAnalyzer.cc:1037
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
trigger::Keys
std::vector< size_type > Keys
Definition: TriggerTypeDefs.h:20
part
part
Definition: HCALResponse.h:20
PhotonAnalyzer::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PhotonAnalyzer.cc:451
reco::PhotonCollection
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
edm::Handle::product
T const * product() const
Definition: Handle.h:81
SiStripMonitorClusterAlca_cfi.xmax
tuple xmax
Definition: SiStripMonitorClusterAlca_cfi.py:35
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
DQMStore::save
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, const uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE", const bool resetMEsAfterWriting=false)
Definition: DQMStore.cc:2540
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:89
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
x
Definition: DDAxes.h:10
DQMStore::bookInt
MonitorElement * bookInt(const char *name)
Book int.
Definition: DQMStore.cc:861
DQMStore::book2D
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:1082
HcalObjRepresent::setup
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Definition: HcalObjRepresent.h:628
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
DQMStore::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:667
PhotonAnalyzer::PhotonAnalyzer
PhotonAnalyzer(const edm::ParameterSet &)
Definition: PhotonAnalyzer.cc:21
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:303
edm::Run
Definition: Run.h:41
GOODCOLL_filter_cfg.cut
tuple cut
Definition: GOODCOLL_filter_cfg.py:65
DQMStore::pwd
const std::string & pwd(void) const
Definition: DQMStore.cc:639
hitfit::scalar
double scalar(const CLHEP::HepGenMatrix &m)
Return the matrix as a scalar. Raise an assertion if the matris is not .
Definition: matutil.cc:183
PhotonAnalyzer::book2DHistoVector
void book2DHistoVector(std::vector< std::vector< MonitorElement * > > &toFill, std::string histoType, std::string histoName, std::string title, int xbin, double xmin, double xmax, int ybin=1, double ymin=1, double ymax=2)
Definition: PhotonAnalyzer.cc:1147
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5