CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ValidIsoTrkCalib.cc
Go to the documentation of this file.
1 //
2 // Package: ValidIsoTrkCalib
3 // Class: ValidIsoTrkCalib
4 //
5 
6 /*
7  Description: Validation for Isolated tracks Calibration
8 
9  Implementation:
10 See the twiki page for details:
11 https://twiki.cern.ch/twiki/bin/view/CMS/ValidIsoTrkCalib
12 
13 */
14 
15 //
16 // Original Author: Andrey Pozdnyakov
17 // Created: Tue Nov 4 01:16:05 CET 2008
18 // $Id: ValidIsoTrkCalib.cc,v 1.14 2012/11/12 21:08:18 dlange Exp $
19 //
20 
21 // system include files
22 
23 #include <memory>
24 
25 // user include files
29 
30 //#include "FWCore/Framework/interface/Event.h"
31 //#include "FWCore/Framework/interface/MakerMacros.h"
32 //#include "FWCore/ParameterSet/interface/ParameterSet.h"
33 
42 
51 
53 
54 #include "TROOT.h"
55 #include "TFile.h"
56 #include "TTree.h"
57 #include "TH1F.h"
58 
59 
60 #include <fstream>
61 #include <map>
62 
63 using namespace edm;
64 using namespace std;
65 using namespace reco;
66 
68 public:
69  explicit ValidIsoTrkCalib(const edm::ParameterSet&);
71 
72  // double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint);
73 
74 private:
75 
76  virtual void beginJob() ;
77  virtual void analyze(const edm::Event&, const edm::EventSetup&);
78  virtual void endJob() ;
79 
80 
81 
82  // ----------member data ---------------------------
83 
84  //Variables from HcalIsotrackAnalyzer
85 
88  double taECALCone_;
89  double taHCALCone_;
90 
91  const CaloGeometry* geo;
95  std::vector<edm::InputTag> genecalLabel_;
100 
101  //std::string m_inputTrackLabel;
102  //std::string m_hcalLabel;
103 
105  string AxB_;
107 
110  //string calibFactorsFileName_;
111  // string corrfile;
112 
114  //bool takeAllRecHits_, takeGenTracks_;
115  int gen, iso, pix;
116  float genPt[500], genPhi[500], genEta[500];
117  float isoPt[500], isoPhi[500], isoEta[500];
118  float pixPt[500], pixPhi[500], pixEta[500];
119  //int hbheiEta[5000],hbheiPhi[5000],hbheDepth[5000];
120  //float hbheEnergy[5000];
121 
122 
123  int NisoTrk;
124  float trackPt, trackE, trackEta, trackPhi;
125  float ptNear;
126  float ptrack, rvert;
127  //float eecal, ehcal;
128 
129  int MinNTrackHitsBarrel, MinNTECHitsEndcap;
130  double energyECALmip, maxPNear;
131  double energyMinIso, energyMaxIso;
132 
133 
134  Float_t emEnergy;
135  Float_t targetE;
136 
137  //TFile* rootFile;
138  // TTree* tree;
139  TTree *tTree, *fTree;
140 
141  Float_t xTrkEcal;
142  Float_t yTrkEcal;
143  Float_t zTrkEcal;
144 
145  Float_t xTrkHcal;
146  Float_t yTrkHcal;
147  Float_t zTrkHcal;
148 
149  int Nhits;
150  float eClustBefore; //Calo energy before calibration
151  float eClustAfter; //After calibration
152  float eTrack; //Track energy
153  float etaTrack;
154  float phiTrack;
155  float eECAL; // Energy deposited in ECAL
156  int numHits; //number of rechits
157  //float e3x3;
158 
166  int iEta;
167  int iPhi;
168  int iEtaTr;
169  int iPhiTr;
170  float iDr, delR;
171  int dietatr;
172  int diphitr;
173 
174  float iTime;
175  float HTime[100];
176  float e3x3Before;
177  float e3x3After;
178  float e5x5Before;
179  float e5x5After;
182  float PtNearBy;
183  float numVH, numVS, numValidTrkHits, numValidTrkStrips;
184 
186  // map<UInt_t, Float_t> CalibFactors;
187  // Bool_t ReadCalibFactors(string);
188 
189  TH1F *nTracks;
190 
192  // int Lumi_n;
193 
194 
195 };
196 
197 
199 
200 {
201 
202  //takeAllRecHits_=iConfig.getUntrackedParameter<bool>("takeAllRecHits");
203  takeGenTracks_=iConfig.getUntrackedParameter<bool>("takeGenTracks");
204 
205  genTracksLabel_ = iConfig.getParameter<edm::InputTag>("genTracksLabel");
206  genhbheLabel_= iConfig.getParameter<edm::InputTag>("genHBHE");
207  //genhoLabel_=iConfig.getParameter<edm::InputTag>("genHO");
208  //genecalLabel_=iConfig.getParameter<std::vector<edm::InputTag> >("genECAL");
209 
210  // m_hcalLabel = iConfig.getUntrackedParameter<std::string> ("hcalRecHitsLabel","hbhereco");
211 
212  hbheLabel_= iConfig.getParameter<edm::InputTag>("hbheInput");
213  hoLabel_=iConfig.getParameter<edm::InputTag>("hoInput");
214  //eLabel_=iConfig.getParameter<edm::InputTag>("eInput");
215  trackLabel_ = iConfig.getParameter<edm::InputTag>("HcalIsolTrackInput");
216  trackLabel1_ = iConfig.getParameter<edm::InputTag>("trackInput");
217 
218  associationConeSize_=iConfig.getParameter<double>("associationConeSize");
219  allowMissingInputs_=iConfig.getUntrackedParameter<bool>("allowMissingInputs", true);
220 // outputFileName_=iConfig.getParameter<std::string>("outputFileName");
221  // calibFactorsFileName_=iConfig.getParameter<std::string>("calibFactorsFileName");
222 
223 
224  AxB_=iConfig.getParameter<std::string>("AxB");
225  calibrationConeSize_=iConfig.getParameter<double>("calibrationConeSize");
226 
227  MinNTrackHitsBarrel = iConfig.getParameter<int>("MinNTrackHitsBarrel");
228  MinNTECHitsEndcap = iConfig.getParameter<int>("MinNTECHitsEndcap");
229  energyECALmip = iConfig.getParameter<double>("energyECALmip");
230  energyMinIso = iConfig.getParameter<double>("energyMinIso");
231  energyMaxIso = iConfig.getParameter<double>("energyMaxIso");
232  maxPNear = iConfig.getParameter<double>("maxPNear");
233 
234  edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
235  parameters_.loadParameters( parameters );
236  trackAssociator_.useDefaultPropagator();
237 
238  // taECALCone_=iConfig.getUntrackedParameter<double>("TrackAssociatorECALCone",0.5);
239  //taHCALCone_=iConfig.getUntrackedParameter<double>("TrackAssociatorHCALCone",0.6);
240 
241 }
242 
243 
245 {
246 
247  // do anything here that needs to be done at desctruction time
248  // (e.g. close files, deallocate resources etc.)
249 
250 }
251 
252 
253 // ------------ method called to for each event ------------
254 void
256 {
257  using namespace edm;
258 
259  try{
260  edm::ESHandle <HcalRespCorrs> recalibCorrs;
261  iSetup.get<HcalRespCorrsRcd>().get("recalibrate",recalibCorrs);
262  respRecalib = recalibCorrs.product();
263 
264  LogInfo("CalibConstants")<<" Loaded: OK ";
265 
266  }catch(const cms::Exception & e) {
267  LogWarning("CalibConstants")<<" Not Found!! ";
268  }
269 
271  iEvent.getByLabel(genTracksLabel_, generalTracks);
272 
274  iEvent.getByLabel(trackLabel1_,isoProdTracks);
275 
276 
278  //edm::Handle<reco::TrackCollection> isoPixelTracks;
279  iEvent.getByLabel(trackLabel_,isoPixelTracks);
280 
281  /*
282  edm::Handle<EcalRecHitCollection> ecal;
283  iEvent.getByLabel(eLabel_,ecal);
284  const EcalRecHitCollection Hitecal = *(ecal.product());
285  */
286 
288  iEvent.getByLabel(hbheLabel_,hbhe);
289  const HBHERecHitCollection Hithbhe = *(hbhe.product());
290 
292  iSetup.get<CaloGeometryRecord>().get(pG);
293  geo = pG.product();
294 
295  const CaloSubdetectorGeometry* gHcal = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
296  //Note: even though it says HcalBarrel, we actually get the whole Hcal detector geometry!
297 
298  // Lumi_n=iEvent.luminosityBlock();
299  parameters_.useEcal = true;
300  parameters_.useHcal = true;
301  parameters_.useCalo = false;
302  parameters_.useMuon = false;
303  //parameters_.dREcal = taECALCone_;
304  //parameters_.dRHcal = taHCALCone_;
305 
306  //cout<<"Hello World. TrackCollectionSize: "<< isoPixelTracks->size()<<endl;
307 
308 if (isoPixelTracks->size()==0) return;
309 
310 
311 for (reco::TrackCollection::const_iterator trit=isoProdTracks->begin(); trit!=isoProdTracks->end(); trit++)
312  {
313 
314  reco::IsolatedPixelTrackCandidateCollection::const_iterator isoMatched=isoPixelTracks->begin();
315  //reco::TrackCollection::const_iterator isoMatched=isoPixelTracks->begin();
316  bool matched=false;
317 
318  //for (reco::IsolatedPixelTrackCandidateCollection::const_iterator trit = isoPixelTracks->begin(); trit!=isoPixelTracks->end(); trit++)
319  for (reco::IsolatedPixelTrackCandidateCollection::const_iterator it = isoPixelTracks->begin(); it!=isoPixelTracks->end(); it++)
320  //for (reco::TrackCollection::const_iterator it = isoPixelTracks->begin(); it!=isoPixelTracks->end(); it++)
321  {
322 
323  if (abs((trit->pt() - it->pt())/it->pt()) < 0.005 && abs(trit->eta() - it->eta()) < 0.01)
324  {
325  isoMatched=it;
326  matched=true;
327  break;
328  }
329  }
330  // CUT
331 
332  if (!matched) continue;
333  if(isoMatched->maxPtPxl() > maxPNear) continue;
334 
335  ptNear = isoMatched->maxPtPxl();
336  //cout<<"Point 0.1 isoMatch. ptnear: "<<ptNear<<endl;
337 
338 
339  // CUT
340  if (trit->hitPattern().numberOfValidHits()<MinNTrackHitsBarrel) continue;
341  if (fabs(trit->eta())>1.47&&trit->hitPattern().numberOfValidStripTECHits()<MinNTECHitsEndcap) continue;
342 
343  //cout<<"Point 0.2.1 after numofvalidhits HB: "<<trit->hitPattern().numberOfValidHits()<<endl;
344  //cout<<"Point 0.2.2 after numofvalidstrips HE: "<<trit->hitPattern().numberOfValidStripTECHits()<<endl;
345 
346  numVH = trit->hitPattern().numberOfValidHits();
347  numVS = trit->hitPattern().numberOfValidStripTECHits();
348 
349 
350 
351  trackE = sqrt(trit->px()*trit->px()+trit->py()*trit->py()+trit->pz()*trit->pz()+0.14*0.14);
352  trackPt = trit->pt();
353  trackEta = trit->eta();
354  trackPhi = trit->phi();
355 
356  emEnergy = isoMatched->energyIn();
357 
358  //cout<<"Point 0.3. Matched :: pt: "<<trit->pt()<<" wholeEnergy: "<<trackE<<" emEnergy: "<<emEnergy<<" eta: "<<etahcal<<" phi: "<<phihcal<<endl;
359  //cout<<"Point 0.4. EM energy in cone: "<<emEnergy<<" EtaHcal: "<<etahcal<<" PhiHcal: "<<phihcal<<endl;
360 
361  TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iSetup,trackAssociator_.getFreeTrajectoryState(iSetup, *trit), parameters_);
362 
363  //float etaecal=info.trkGlobPosAtEcal.eta();
364  //float phiecal=info.trkGlobPosAtEcal.phi();
365  // float etahcal=info.trkGlobPosAtHcal.eta();
366  // float phihcal=info.trkGlobPosAtHcal.phi();
367 
368 
369  xTrkEcal=info.trkGlobPosAtEcal.x();
370  yTrkEcal=info.trkGlobPosAtEcal.y();
371  zTrkEcal=info.trkGlobPosAtEcal.z();
372 
373  xTrkHcal=info.trkGlobPosAtHcal.x();
374  yTrkHcal=info.trkGlobPosAtHcal.y();
375  zTrkHcal=info.trkGlobPosAtHcal.z();
376 
377  if (xTrkEcal==0 && yTrkEcal==0&& zTrkEcal==0) {cout<<"zero point at Ecal"<<endl; continue;}
378  if (xTrkHcal==0 && yTrkHcal==0&& zTrkHcal==0) {cout<<"zero point at Hcal"<<endl; continue;}
379 
380  /*GlobalVector trackMomAtEcal = info.trkMomAtEcal;
381  GlobalVector trackMomAtHcal = info.trkMomAtHcal;
382 
383  PxTrkHcal = trackMomAtHcal.x();
384  PyTrkHcal = trackMomAtHcal.y();
385  PzTrkHcal = trackMomAtHcal.z();
386  */
387 
388  GlobalPoint gPointEcal(xTrkEcal,yTrkEcal,zTrkEcal);
389  GlobalPoint gPointHcal(xTrkHcal,yTrkHcal,zTrkHcal);
390 
391  int iphitrue = -10;
392  int ietatrue = 100;
393  const HcalDetId tempId = gHcal->getClosestCell(gPointHcal);
394  ietatrue = tempId.ieta();
395  iphitrue = tempId.iphi();
396 
397 
398  MaxHit_struct MaxHit;
399 
400  MaxHit.hitenergy=-100.;
401 
402  //container for used recHits
403  std::vector<DetId> usedHits;
404  //
405  usedHits.clear();
406  //cout <<"Point 1. Entrance to HBHECollection"<<endl;
407 
408  //float dddeta = 1000.;
409  //float dddphi = 1000.;
410  //int iphitrue = 1234;
411  //int ietatrue = 1234;
412 
413  GlobalPoint gPhot;
414 
415 
416  for (HBHERecHitCollection::const_iterator hhit=Hithbhe.begin(); hhit!=Hithbhe.end(); hhit++)
417  {
418 
419  //check that this hit was not considered before and push it into usedHits
420  bool hitIsUsed=false;
421  for (uint32_t i=0; i<usedHits.size(); i++)
422  {
423  if (usedHits[i]==hhit->id()) hitIsUsed=true;
424  }
425  if (hitIsUsed) continue;
426  usedHits.push_back(hhit->id());
427  //
428 
429  // rof 16.05.2008 start: include the possibility for recalibration
430  float recal = 1;
431  // rof end
432 
433  GlobalPoint pos = geo->getPosition(hhit->detid());
434  //float phihit = pos.phi();
435  //float etahit = pos.eta();
436 
437  int iphihitm = (hhit->id()).iphi();
438  int ietahitm = (hhit->id()).ieta();
439  int depthhit = (hhit->id()).depth();
440  float enehit = hhit->energy() * recal;
441 
442  if (depthhit!=1) continue;
443 
444  /*
445  float dphi = fabs(phihcal - phihit);
446  if(dphi > 4.*atan(1.)) dphi = 8.*atan(1.) - dphi;
447  float deta = fabs(etahcal - etahit);
448  float dr = sqrt(dphi*dphi + deta*deta);
449  */
450 
451 
452  //double distAtHcal = getDistInPlaneTrackDir(gPointHcal, trackMomAtHcal, pos);
453  double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
454 
455  if(distAtHcal < associationConeSize_)
456  {
457 
458  for (HBHERecHitCollection::const_iterator hhit2=Hithbhe.begin(); hhit2!=Hithbhe.end(); hhit2++)
459  {
460  int iphihitm2 = (hhit2->id()).iphi();
461  int ietahitm2 = (hhit2->id()).ieta();
462  int depthhit2 = (hhit2->id()).depth();
463  float enehit2 = hhit2->energy() * recal;
464 
465  if (iphihitm==iphihitm2 && ietahitm==ietahitm2 && depthhit!=depthhit2) enehit = enehit+enehit2;
466 
467  }
468 
469 
470  //cout<<"IN CONE ieta: "<<ietahitm<<" iphi: "<<iphihitm<<" depthhit: "<<depthhit<<" dr: "<<dr<<" energy: "<<enehit<<endl;
471 
472  //Find a Hit with Maximum Energy
473 
474  if(enehit > MaxHit.hitenergy)
475  {
476  MaxHit.hitenergy = enehit;
477  MaxHit.ietahitm = (hhit->id()).ieta();
478  MaxHit.iphihitm = (hhit->id()).iphi();
479  MaxHit.dr = distAtHcal;
480  //MaxHit.depthhit = (hhit->id()).depth();
481  MaxHit.depthhit = 1;
482 
483  //gPhot = geo->getPosition(hhit->detid());
484  }
485 
486  }
487  } //end of all HBHE hits cycle
488 
489  usedHits.clear();
490 
491  //cout<<"Hottest ieta: "<<MaxHit.ietahitm<<" iphi: "<<MaxHit.iphihitm<<" dr: "<<MaxHit.dr<<endl;
492  //cout<<"Track ieta: "<<ietatrue<<" iphi: "<<iphitrue<<endl;
493 
494  //cout<<"Point 3. MaxHit :::En "<<MaxHit.hitenergy<<" ieta: "<<MaxHit.ietahitm<<" iphi: "<<MaxHit.iphihitm<<endl;
495 
496  Bool_t passCuts = kFALSE;
497  if(trackE > energyMinIso && trackE < energyMaxIso && emEnergy < energyECALmip && MaxHit.hitenergy > 0. && abs(MaxHit.ietahitm)<29) passCuts = kTRUE;
498 
499  //cout<<"Pont 0.1.1. trackE:"<<trackE<<" emEn: "<<emEnergy<<endl;
500 
501 
502  numHits=0;
503 
504  eClustBefore = 0.0;
505  eClustAfter = 0.0;
506  eBeforeDepth1 = 0.0;
507  eAfterDepth1 = 0.0;
508  eBeforeDepth2 = 0.0;
509  eAfterDepth2 = 0.0;
510  CentHitFactor = 0.0;
511  e3x3After = 0.0;
512  e3x3Before = 0.0;
513  e5x5After = 0.0;
514  e5x5Before = 0.0;
515 
516 
517  for (HBHERecHitCollection::const_iterator hhit=Hithbhe.begin(); hhit!=Hithbhe.end(); hhit++)
518  {
519 
520  //check that this hit was not considered before and push it into usedHits
521  bool hitIsUsed=false;
522  for (uint32_t i=0; i<usedHits.size(); i++)
523  {
524  if (usedHits[i]==hhit->id()) hitIsUsed=true;
525  }
526  if (hitIsUsed) continue;
527  usedHits.push_back(hhit->id());
528 
529  int DIETA = 100;
530  if(MaxHit.ietahitm*(hhit->id()).ieta()>0)
531  {
532  DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
533  }
534  if(MaxHit.ietahitm*(hhit->id()).ieta()<0)
535  {
536  DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
537  DIETA = DIETA>0 ? DIETA-1 : DIETA+1;
538  }
539 
540  int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
541  DIPHI = DIPHI>36 ? 72-DIPHI : DIPHI;
542 
543 
544 
545  int numbercell=100; //always collect Wide clastor!
546 
547  //if(AxB_=="5x5" || AxB_=="3x3" || AxB_=="7x7"|| AxB_=="Cone")
548 
549  //if(AxB_=="3x3") numbercell = 1;
550  //if(AxB_=="5x5") numbercell = 2;
551  //if(AxB_=="Cone") numbercell = 1000;
552 
553  if( abs(DIETA)<=numbercell && (abs(DIPHI)<=numbercell || ( abs(MaxHit.ietahitm)>=20 && abs(DIPHI)<=numbercell+1)) )
554  {
555  const GlobalPoint pos2 = geo->getPosition(hhit->detid());
556 
557  if(passCuts && hhit->energy()>0)
558  {
559 
560 
561  float factor = 0.0;
562  // factor = CalibFactors[hhit->id()];
563  factor = respRecalib -> getValues(hhit->id())->getValue();
564 
565 
566 
567  //if(i<5){cout<<" calib factors: "<<factor<<" ij "<<100*i+j<<endl;}
568 
569  if (hhit->id().ieta() == MaxHit.ietahitm && hhit->id().iphi() == MaxHit.iphihitm) CentHitFactor=factor;
570 
571  if (hhit->id().ieta() == MaxHit.ietahitm && hhit->id().iphi() == MaxHit.iphihitm) iTime = hhit->time();
572 
573  if(AxB_!="3x3" && AxB_!="5x5" && AxB_!="Cone") LogWarning(" AxB ")<<" Not supported: "<< AxB_;
574 
575 
576  if (abs(DIETA)<=2 && (abs(DIPHI)<=2 || ((abs(MaxHit.ietahitm)>20 && abs(DIPHI)<=4) && !( (abs(MaxHit.ietahitm)==21 || abs(MaxHit.ietahitm)==22) && abs((hhit->id()).ieta())<=20 && abs(DIPHI)>2) )) )
577  {
578 
579  e5x5Before += hhit->energy();
580  e5x5After += hhit->energy()*factor;
581  }
582 
583 
584  if (abs(DIETA)<=1 && (abs(DIPHI)<=1 || ((abs(MaxHit.ietahitm)>20 && abs(DIPHI)<=2) && !(abs(MaxHit.ietahitm)==21 && abs((hhit->id()).ieta())<=20 && abs(DIPHI)>1) )) )
585  {
586 
587  e3x3Before += hhit->energy();
588  e3x3After += hhit->energy()*factor;
589  }
590 
591 
592  if(AxB_=="5x5")
593  {
594 
595  if (abs(DIETA)<=2 && (abs(DIPHI)<=2 || ( abs(MaxHit.ietahitm)>20 && abs(DIPHI)<=4) ) )
596  {
597  if (abs(MaxHit.ietahitm)==21 && abs((hhit->id()).ieta())<=20 && abs(DIPHI)>3) continue;
598 
599  HTime[numHits]= hhit->time();
600  numHits++;
601 
602  eClustBefore += hhit->energy();
603  eClustAfter += hhit->energy()*factor;
604 
605  if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29))
606  {
607  eBeforeDepth1 += hhit->energy();
608  eAfterDepth1 += hhit->energy()*factor;
609  }
610  else if((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29))
611  {
612  eBeforeDepth2 += hhit->energy();
613  eAfterDepth2 += hhit->energy()*factor;
614  }
615  }
616  }//end of 5x5
617 
618 
619 
620  if(AxB_=="3x3")
621  {
622  if (abs(DIETA)<=1 && (abs(DIPHI)<=1 || ( abs(MaxHit.ietahitm)>20 && abs(DIPHI)<=2) ) )
623  {
624  if (abs(MaxHit.ietahitm)==21 && abs((hhit->id()).ieta())<=20 && abs(DIPHI)>2) continue;
625 
626  HTime[numHits]= hhit->time();
627  numHits++;
628 
629  eClustBefore += hhit->energy();
630  eClustAfter += hhit->energy() * factor;
631 
632  if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29))
633  {
634  eBeforeDepth1 += hhit->energy();
635  eAfterDepth1 += hhit->energy() * factor;
636  }
637  else if((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29))
638  {
639  eBeforeDepth2 += hhit->energy();
640  eAfterDepth2 += hhit->energy() * factor;
641  }
642 
643  }
644  }//end of 3x3
645 
646 
647  if (AxB_=="Cone" && getDistInPlaneSimple(gPointHcal, pos2) < calibrationConeSize_) {
648 
649  HTime[numHits]= hhit->time();
650  numHits++;
651 
652  eClustBefore += hhit->energy();
653  eClustAfter += hhit->energy() * factor;
654 
655  if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29))
656  {
657  eBeforeDepth1 += hhit->energy();
658  eAfterDepth1 += hhit->energy() * factor;
659  }
660  else if((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29))
661  {
662  eBeforeDepth2 += hhit->energy();
663  eAfterDepth2 += hhit->energy() * factor;
664  }
665 
666 
667  }//end of Cone
668 
669  }//end of passCuts
670 
671  }//end of DIETA DIPHI
672 
673  } //end of associatedcone HBHE hits cycle
674 
675 
676  int dieta_M_P = 100;
677  int diphi_M_P = 100;
678  if(MaxHit.ietahitm*ietatrue>0) {dieta_M_P = abs (MaxHit.ietahitm-ietatrue);}
679  if(MaxHit.ietahitm*ietatrue<0) {dieta_M_P = abs(MaxHit.ietahitm-ietatrue)-1;}
680  diphi_M_P = abs(MaxHit.iphihitm-iphitrue);
681  diphi_M_P = diphi_M_P>36 ? 72-diphi_M_P : diphi_M_P;
682 
683 
684  if(passCuts)
685 
686  {
687  eventNumber = iEvent.id().event();
688  runNumber = iEvent.id().run();
689 
690  eCentHitBefore = MaxHit.hitenergy;
691  eCentHitAfter = MaxHit.hitenergy*CentHitFactor;
692  eECAL = emEnergy;
693  numValidTrkHits = numVH;
694  numValidTrkStrips = numVS;
695  PtNearBy = ptNear;
696 
697 
698  eTrack = trackE;
699  phiTrack = trackPhi;
700  etaTrack = trackEta;
701 
702  iEta = MaxHit.ietahitm;
703  iPhi = MaxHit.iphihitm;
704 
705  iEtaTr = ietatrue;
706  iPhiTr = iphitrue;
707  iDr = sqrt(diphi_M_P*diphi_M_P+dieta_M_P*dieta_M_P);
708  delR = MaxHit.dr;
709  dietatr = dieta_M_P;
710  diphitr = diphi_M_P;
711 
712  fTree -> Fill();
713  }
714 
715  } //end of isoProdTracks cycle
716 
717 
718 
719 /* ------------------ Some stuff for general tracks ---------- ----*/
720  //cout<<" generalTracks Size: "<< generalTracks->size()<<endl;
721  int n = generalTracks->size();
722  nTracks->Fill(n);
723 
724  if(takeGenTracks_ && iEvent.id().event()%10==1)
725  {
726  gen = generalTracks->size();
727  iso = isoProdTracks->size();
728  pix = isoPixelTracks->size();
729 
730  genPt[0] = -33;
731  genPhi[0] = -33;
732  genEta[0] = -33;
733 
734  isoPt[0] = -33;
735  isoPhi[0] = -33;
736  isoEta[0] = -33;
737 
738  pixPt[0] = -33;
739  pixPhi[0] = -33;
740  pixEta[0] = -33;
741 
742  Int_t gencount=0, isocount=0, pixcount=0;
743  for (reco::TrackCollection::const_iterator gentr=generalTracks->begin(); gentr!=generalTracks->end(); gentr++)
744  {
745  genPt[gencount] = gentr->pt();
746  genPhi[gencount] = gentr->phi();
747  genEta[gencount] = gentr->eta();
748  gencount++;
749  }
750 
751  for (reco::TrackCollection::const_iterator isotr=isoProdTracks->begin(); isotr!=isoProdTracks->end(); isotr++)
752  {
753  isoPt[isocount] = isotr->pt();
754  isoPhi[isocount] = isotr->phi();
755  isoEta[isocount] = isotr->eta();
756  isocount++;
757  }
758 
759  for (reco::IsolatedPixelTrackCandidateCollection::const_iterator pixtr=isoPixelTracks->begin(); pixtr!=isoPixelTracks->end(); pixtr++)
760  {
761  pixPt[pixcount] = pixtr->pt();
762  pixPhi[pixcount] = pixtr->phi();
763  pixEta[pixcount] = pixtr->eta();
764  pixcount++;
765  }
766  }
767 
768  tTree -> Fill();
769 
770 }
771 
772 // ------------ method called once each job just before starting event loop ------------
773 void
775 {
776 
777  // if(!ReadCalibFactors(calibFactorsFileName_.c_str() )) {cout<<"Cant read file with cailib coefficients!! ---"<<endl;}
778 
779 // try{
780 // edm::ESHandle <HcalRespCorrs> recalibCorrs;
781 // iSetup.get<HcalRespCorrsRcd>().get("recalibrate",recalibCorrs);
782 // respRecalib = recalibCorrs.product();
783 //
784 // LogInfo("CalibConstants")<<" Loaded: OK ";
785 //
786 // }catch(const cms::Exception & e) {
787 // LogWarning("CalibConstants")<<" Not Found!! ";
788 // }
789 
790 
791  // rootFile = new TFile(outputFileName_.c_str(),"RECREATE");
792 
793  //@@@@@@@@@@@@@
794  //TFileDirectory ValDir = fs->mkdir("Validation");
795 
796  nTracks = fs->make<TH1F>("nTracks","general;number of general tracks",11,-0.5,10.5);
797 
798 
799  tTree = fs->make<TTree>("tTree", "Tree for gen info");
800 
801  fTree = fs->make<TTree>("fTree", "Tree for IsoTrack Calibration");
802 
803  fTree->Branch("eventNumber", &eventNumber, "eventNumber/I");
804  fTree->Branch("runNumber", &runNumber, "runNumber/I");
805 
806  fTree->Branch("eClustBefore", &eClustBefore,"eClustBefore/F");
807  fTree->Branch("eClustAfter", &eClustAfter,"eClustAfter/F");
808  fTree->Branch("eTrack", &eTrack, "eTrack/F");
809  fTree->Branch("etaTrack", &etaTrack, "etaTrack/F");
810  fTree->Branch("phiTrack", &phiTrack, "phiTrack/F");
811 
812  fTree->Branch("numHits", &numHits, "numHits/I");
813  fTree->Branch("eECAL", &eECAL, "eECAL/F");
814  fTree->Branch("PtNearBy", &PtNearBy, "PtNearBy/F");
815  fTree->Branch("numValidTrkHits", &numValidTrkHits, "numValidTrkHits/F");
816  fTree->Branch("numValidTrkStrips", &numValidTrkStrips, "numValidTrkStrips/F");
817 
818  fTree->Branch("eBeforeDepth1", &eBeforeDepth1, "eBeforeDepth1/F");
819  fTree->Branch("eBeforeDepth2", &eBeforeDepth2, "eBeforeDepth2/F");
820  fTree->Branch("eAfterDepth1", &eAfterDepth1, "eAfterDepth1/F");
821  fTree->Branch("eAfterDepth2", &eAfterDepth2, "eAfterDepth2/F");
822 
823  fTree->Branch("e3x3Before", &e3x3Before, "e3x3Before/F");
824  fTree->Branch("e3x3After", &e3x3After, "e3x3After/F");
825  fTree->Branch("e5x5Before", &e5x5Before, "e5x5Before/F");
826  fTree->Branch("e5x5After", &e5x5After, "e5x5After/F");
827 
828  fTree->Branch("eCentHitAfter", &eCentHitAfter, "eCentHitAfter/F");
829  fTree->Branch("eCentHitBefore", &eCentHitBefore, "eCentHitBefore/F");
830  fTree->Branch("iEta", &iEta, "iEta/I");
831  fTree->Branch("iPhi", &iPhi, "iPhi/I");
832 
833 
834  fTree->Branch("iEtaTr", &iEtaTr, "iEtaTr/I");
835  fTree->Branch("iPhiTr", &iPhiTr, "iPhiTr/I");
836  fTree->Branch("dietatr", &dietatr, "dietatr/I");
837  fTree->Branch("diphitr", &diphitr, "diphitr/I");
838  fTree->Branch("iDr", &iDr, "iDr/F");
839  fTree->Branch("delR", &delR, "delR/F");
840 
841  fTree->Branch("iTime", &iTime, "iTime/F");
842  fTree->Branch("HTime", HTime, "HTime[numHits]/F");
843 
844  fTree->Branch("xTrkEcal", &xTrkEcal, "xTrkEcal/F");
845  fTree->Branch("yTrkEcal", &yTrkEcal, "yTrkEcal/F");
846  fTree->Branch("zTrkEcal", &zTrkEcal, "zTrkEcal/F");
847  fTree->Branch("xTrkHcal", &xTrkHcal, "xTrkHcal/F");
848  fTree->Branch("yTrkHcal", &yTrkHcal, "yTrkHcal/F");
849  fTree->Branch("zTrkHcal", &zTrkHcal, "zTrkHcal/F");
850 
851  if(takeGenTracks_) {
852  tTree->Branch("gen", &gen, "gen/I");
853  tTree->Branch("iso", &iso, "iso/I");
854  tTree->Branch("pix", &pix, "pix/I");
855  tTree->Branch("genPt", genPt, "genPt[gen]/F");
856  tTree->Branch("genPhi", genPhi, "genPhi[gen]/F");
857  tTree->Branch("genEta", genEta, "genEta[gen]/F");
858 
859  tTree->Branch("isoPt", isoPt, "isoPt[iso]/F");
860  tTree->Branch("isoPhi", isoPhi, "isoPhi[iso]/F");
861  tTree->Branch("isoEta", isoEta, "isoEta[iso]/F");
862 
863  tTree->Branch("pixPt", pixPt, "pixPt[pix]/F");
864  tTree->Branch("pixPhi", pixPhi, "pixPhi[pix]/F");
865  tTree->Branch("pixEta", pixEta, "pixEta[pix]/F");
866  }
867 
868 
869 }
870 
871 
872 
873 // ------------ method called once each job just after ending the event loop ------------
874 void
876 {
877  // rootFile->Write();
878  //rootFile->Close();
879 }
880 
881 //define this as a plug-in
RunNumber_t run() const
Definition: EventID.h:42
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:44
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
dictionary parameters
Definition: Parameters.py:2
tuple trackPt
Definition: listHistos.py:120
virtual void beginJob()
std::vector< T >::const_iterator const_iterator
virtual void endJob()
#define abs(x)
Definition: mlp_lapack.h:159
TrackAssociatorParameters parameters_
DEFINE_FWK_MODULE(HiMixingModule)
math::XYZPoint trkGlobPosAtHcal
void beginJob()
Definition: Breakpoints.cc:15
int iEvent
Definition: GenABIO.cc:243
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
const CaloGeometry * geo
T sqrt(T t)
Definition: SSEVec.h:48
int ieta() const
get the cell ieta
Definition: HcalDetId.h:38
TrackDetectorAssociator trackAssociator_
const HcalRespCorrs * respRecalib
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
ValidIsoTrkCalib(const edm::ParameterSet &)
std::vector< edm::InputTag > genecalLabel_
int iphi() const
get the cell iphi
Definition: HcalDetId.h:40
const T & get() const
Definition: EventSetup.h:55
edm::Service< TFileService > fs
T const * product() const
Definition: ESHandle.h:62
edm::EventID id() const
Definition: EventBase.h:56
virtual void analyze(const edm::Event &, const edm::EventSetup &)
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
const JetExtendedData & getValue(const Container &, const reco::JetBaseRef &)
get value for the association. Throw exception if no association found
double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint)
Definition: ConeDefinition.h:9
tuple cout
Definition: gather_cfg.py:121