CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalCorrPFCalculation.cc
Go to the documentation of this file.
1 
2 
3 #include "FWCore/Framework/interface/EDAnalyzer.h"
4 #include "FWCore/Utilities/interface/InputTag.h"
5 
6 #include "Geometry/Records/interface/CaloGeometryRecord.h"
7 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
8 
9 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
10 #include "DataFormats/HcalRecHit/interface/HcalSourcePositionData.h"
11 #include <DataFormats/EcalDetId/interface/EBDetId.h>
12 #include <DataFormats/EcalDetId/interface/EEDetId.h>
13 
14 #include "DataFormats/HcalDetId/interface/HcalElectronicsId.h"
15 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
16 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
17 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
18 
19 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
20 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
21 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
22 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
23 #include "CondFormats/HcalObjects/interface/HcalPFCorrs.h"
24 #include "CondFormats/DataRecord/interface/HcalPFCorrsRcd.h"
25 #include "MagneticField/Engine/interface/MagneticField.h"
26 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
27 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
28 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
29 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
30 #include "DataFormats/GeometrySurface/interface/Plane.h"
31 
32 
33 #include "Calibration/HcalCalibAlgos/interface/CommonUsefulStuff.h"
34 
35 #include <iostream>
36 #include "TProfile.h"
37 
38 using namespace edm;
39 using namespace std;
40 using namespace reco;
41 
42 class HcalCorrPFCalculation : public edm::EDAnalyzer {
43  public:
44  HcalCorrPFCalculation(edm::ParameterSet const& conf);
45  ~HcalCorrPFCalculation();
46  virtual void analyze(edm::Event const& ev, edm::EventSetup const& c) override;
47  virtual void beginJob() override ;
48  virtual void endJob() override ;
49  private:
50 
51  double RecalibFactor(HcalDetId id);
52 
53  bool Respcorr_;
54  bool PFcorr_;
55  bool Conecorr_;
56  double radius_;
57 
58  double clusterConeSize_,associationConeSize_, ecalCone_, neutralIsolationCone_, trackIsolationCone_;
59  float eECAL, eECAL09cm, eECAL40cm;
60 // double energyECALmip;
61 
62  float eCentHit, eTrack, eParticle;
63  float CentHitFactor;
64  int iEta, iPhi;
65  //int iEtaTr, iPhiTr;
66  float iDr, delR;
67  float e3x3, e5x5;
68  float phiTrack, etaTrack;
69  float phiGPoint, etaGPoint;
70 
71  Bool_t doHF;
72  Bool_t AddRecalib;
73  int nevtot;
74 
75  const HcalRespCorrs* respRecalib;
76  const HcalPFCorrs* pfRecalib;
77 
78  SteppingHelixPropagator* stepPropF;
79  MagneticField *theMagField;
80 
81  TrackDetectorAssociator trackAssociator_;
82  TrackAssociatorParameters parameters_;
83 
84  const CaloGeometry* geo;
85  const HcalGeometry* geoHcal;
86 
87  Float_t xTrkHcal,yTrkHcal,zTrkHcal;
88  Float_t xTrkEcal, yTrkEcal, zTrkEcal;
89  Float_t xAtHcal, yAtHcal, zAtHcal;
90 
91  float eEcalCone, eHcalCone, eHcalConeNoise;
92  int UsedCells, UsedCellsNoise;
93  float phiParticle, etaParticle;
94 
95  Int_t numValidTrkHits[50], numValidTrkStrips[50], numLayers[50];
96 
97  Bool_t trkQual[50];
98 
99  TProfile *nCells, *nCellsNoise, *enHcal, *enHcalNoise;
100  // TH1F *enEcalB, *enEcalE;
101  TTree *pfTree, *tracksTree;
102  // TFile *rootFile;
103  edm::Service<TFileService> fs;
104  Int_t nTracks;
105  Float_t genEta,genPhi, trackEta[50],trackPhi[50], trackP[50] , delRmc[50];
106 
107  edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
108  edm::EDGetTokenT<HORecHitCollection> tok_ho_;
109  edm::EDGetTokenT<HFRecHitCollection> tok_hf_;
110 
111  edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
112  edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
113  edm::EDGetTokenT<reco::TrackCollection> tok_tracks_;
114  edm::EDGetTokenT<edm::HepMCProduct> tok_gen_;
115 };
116 
117 
118 HcalCorrPFCalculation::HcalCorrPFCalculation(edm::ParameterSet const& iConfig) {
119 
120  tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheRecHitCollectionTag"));
121  tok_hf_ = consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfRecHitCollectionTag"));
122  tok_ho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoRecHitCollectionTag"));
123 
124  // should maybe add these options to configuration - cowden
125  tok_EE_ = consumes<EcalRecHitCollection>( edm::InputTag("ecalRecHit","EcalRecHitsEE") );
126  tok_EB_ = consumes<EcalRecHitCollection>( edm::InputTag("ecalRecHit","EcalRecHitsEB") );
127  tok_tracks_ = consumes<reco::TrackCollection>( edm::InputTag("generalTracks") );
128  tok_gen_ = consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"));
129 
130  // outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile", "myfile.root");
131 
132  Respcorr_ = iConfig.getUntrackedParameter<bool>("RespcorrAdd", false);
133  PFcorr_ = iConfig.getUntrackedParameter<bool>("PFcorrAdd", false);
134  Conecorr_ = iConfig.getUntrackedParameter<bool>("ConeCorrAdd", true);
135  //radius_ = iConfig.getUntrackedParameter<double>("ConeRadiusCm", 40.);
136  //energyECALmip = iConfig.getParameter<double>("energyECALmip");
137 
138  edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
139  edm::ConsumesCollector iC = consumesCollector();
140  parameters_.loadParameters( parameters, iC );
141  trackAssociator_.useDefaultPropagator();
142 
143  associationConeSize_=iConfig.getParameter<double>("associationConeSize");
144  clusterConeSize_=iConfig.getParameter<double>("clusterConeSize");
145  ecalCone_=iConfig.getParameter<double>("ecalCone");
146  trackIsolationCone_ = iConfig.getParameter<double>("trackIsolationCone");
147  neutralIsolationCone_ = iConfig.getParameter<double>("neutralIsolationCone");
148 
149  // AxB_=iConfig.getParameter<std::string>("AxB");
150 
151 }
152 
153 
154 double HcalCorrPFCalculation::RecalibFactor(HcalDetId id)
155 {
156  Float_t resprecal = 1.;
157  Float_t pfrecal = 1.;
158  if(AddRecalib) {
159  if(Respcorr_) resprecal = respRecalib -> getValues(id)->getValue();
160  if(PFcorr_) pfrecal = pfRecalib -> getValues(id)->getValue();
161  }
162  Float_t factor = resprecal*pfrecal;
163  return factor;
164 }
165 
166 HcalCorrPFCalculation::~HcalCorrPFCalculation() {}
167 
168 void HcalCorrPFCalculation::analyze(edm::Event const& ev, edm::EventSetup const& c) {
169 
170  AddRecalib=kFALSE;
171 
172  try{
173 
174  edm::ESHandle <HcalRespCorrs> recalibCorrs;
175  c.get<HcalRespCorrsRcd>().get("recalibrate",recalibCorrs);
176  respRecalib = recalibCorrs.product();
177 
178  edm::ESHandle <HcalPFCorrs> pfCorrs;
179  c.get<HcalPFCorrsRcd>().get("recalibrate",pfCorrs);
180  pfRecalib = pfCorrs.product();
181 
182  AddRecalib = kTRUE;
183  // LogMessage("CalibConstants")<<" OK ";
184 
185  }catch(const cms::Exception & e) {
186  LogWarning("CalibConstants")<<" Not Found!! ";
187  }
188 
189  edm::Handle<HBHERecHitCollection> hbhe;
190  ev.getByToken(tok_hbhe_, hbhe);
191  const HBHERecHitCollection Hithbhe = *(hbhe.product());
192 
193  edm::Handle<HFRecHitCollection> hfcoll;
194  ev.getByToken(tok_hf_, hfcoll);
195  const HFRecHitCollection Hithf = *(hfcoll.product());
196 
197  edm::Handle<HORecHitCollection> hocoll;
198  ev.getByToken(tok_ho_, hocoll);
199  const HORecHitCollection Hitho = *(hocoll.product());
200 
201  edm::Handle<EERecHitCollection> ecalEE;
202  ev.getByToken(tok_EE_,ecalEE);
203  const EERecHitCollection HitecalEE = *(ecalEE.product());
204 
205  edm::Handle<EBRecHitCollection> ecalEB;
206  ev.getByToken(tok_EB_,ecalEB);
207  const EBRecHitCollection HitecalEB = *(ecalEB.product());
208 
209  // temporary collection of EB+EE recHits
210  std::auto_ptr<EcalRecHitCollection> tmpEcalRecHitCollection(new EcalRecHitCollection);
211  for(EcalRecHitCollection::const_iterator recHit = (*ecalEB).begin(); recHit != (*ecalEB).end(); ++recHit)
212  {tmpEcalRecHitCollection->push_back(*recHit);}
213  for(EcalRecHitCollection::const_iterator recHit = (*ecalEE).begin(); recHit != (*ecalEE).end(); ++recHit)
214  {tmpEcalRecHitCollection->push_back(*recHit);}
215  const EcalRecHitCollection Hitecal = *tmpEcalRecHitCollection;
216 
217 
218  edm::Handle<reco::TrackCollection> generalTracks;
219  ev.getByToken(tok_tracks_, generalTracks);
220 
221  edm::ESHandle<CaloGeometry> pG;
222  c.get<CaloGeometryRecord>().get(pG);
223  geo = pG.product();
224 
225  /*
226  edm::ESHandle<HcalGeometry> hcalG;
227  c.get<HcalGeometryRecord>().get(hcalG);
228  geoHcal = hcalG.product();
229  */
230 
231  const CaloSubdetectorGeometry* gHcal = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
232 
233  parameters_.useEcal = true;
234  parameters_.useHcal = true;
235  parameters_.useCalo = false;
236  parameters_.useMuon = false;
237  parameters_.dREcal = 0.5;
238  parameters_.dRHcal = 0.6;
239 
240  //parameters_.dREcal = taECALCone_;
241  //parameters_.dRHcal = taHCALCone_;
242 
243 
244  edm::ESHandle<MagneticField> bField;
245  c.get<IdealMagneticFieldRecord>().get(bField);
246  stepPropF = new SteppingHelixPropagator(&*bField,alongMomentum);
247  stepPropF->setMaterialMode(false);
248  stepPropF->applyRadX0Correction(true);
249 
250 
251  // double eta_bin[42]={0.,.087,.174,.261,.348,.435,.522,.609,.696,.783,
252  //.870,.957,1.044,1.131,1.218,1.305,1.392,1.479,1.566,1.653,1.740,1.830,1.930,2.043,2.172,
253  //2.322,2.500,2.650,2.853,3.000,3.139,3.314,3.489,3.664,3.839,4.013,4.191,4.363,4.538,4.716,4.889,5.191};
254 
255  // MC info
256  //double phi_MC = -999999.; // phi of initial particle from HepMC
257  //double eta_MC = -999999.; // eta of initial particle from HepMC
258  double mom_MC = 50.; // P of initial particle from HepMC
259  //bool MC = false;
260 
261  // MC information
262 
263 
264  edm::Handle<edm::HepMCProduct> evtMC;
265  // ev.getByLabel("generatorSmeared",evtMC);
266  ev.getByToken(tok_gen_,evtMC);
267  if (!evtMC.isValid())
268  {
269  std::cout << "no HepMCProduct found" << std::endl;
270  }
271  else
272  {
273  //MC=true;
274  // std::cout << "*** source HepMCProduct found"<< std::endl;
275  }
276 
277  // MC particle with highest pt is taken as a direction reference
278  double maxPt = -99999.;
279  int npart = 0;
280 
281  GlobalPoint initpos (0,0,0);
282 
283  HepMC::GenEvent * myGenEvent = new HepMC::GenEvent(*(evtMC->GetEvent()));
284  for ( HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin();
285  p != myGenEvent->particles_end(); ++p )
286  {
287  phiParticle = (*p)->momentum().phi();
288  etaParticle = (*p)->momentum().eta();
289  double pt = (*p)->momentum().perp();
290  mom_MC = (*p)->momentum().rho();
291  if(pt > maxPt) { npart++; maxPt = pt; /*phi_MC = phiParticle; eta_MC = etaParticle;*/ }
292  GlobalVector mom ((*p)->momentum().x(),(*p)->momentum().y(),(*p)->momentum().z());
293  int charge = -1;
294 
295  if(abs((*p)->pdg_id())==211) charge = (*p)->pdg_id()/abs((*p)->pdg_id()); // pions only !!!
296  else continue;
297 
298 
299  /* Propogate the partoicle to Hcal */
300 
301 
302  const FreeTrajectoryState *freetrajectorystate_ =
303  new FreeTrajectoryState(initpos, mom ,charge , &(*theMagField));
304 
305 
306  TrackDetMatchInfo info = trackAssociator_.associate(ev, c, *freetrajectorystate_ , parameters_);
307 
308  GlobalPoint barrelMC(0,0,0), endcapMC(0,0,0), forwardMC1(0,0,0), forwardMC2(0,0,0);
309 
310  GlobalPoint gPointHcal(0.,0.,0.);
311 
312  /*
313  xTrkHcal=info.trkGlobPosAtHcal.x();
314  yTrkHcal=info.trkGlobPosAtHcal.y();
315  zTrkHcal=info.trkGlobPosAtHcal.z();
316  //GlobalPoint gPointHcal(xTrkHcal,yTrkHcal,zTrkHcal);
317 
318  GlobalVector trackMomAtHcal = info.trkMomAtHcal;
319 
320  if (xTrkHcal==0 && yTrkHcal==0 && zTrkHcal==0) continue;
321  */
322 
323 
324  if(fabs(etaParticle)<1.392) {
325  Cylinder *cylinder = new Cylinder(181.1, Surface::PositionType(0,0,0),
326  Surface::RotationType());
327 
328  TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*cylinder));
329  if(steppingHelixstateinfo_.isValid() )
330  { barrelMC = steppingHelixstateinfo_.freeState()->position(); }
331 
332  }
333 
334 
335  doHF = kFALSE;
336  if(fabs(etaParticle)>=1.392 && fabs(etaParticle)<5.191) {
337 
338  Surface::RotationType rot(GlobalVector(1,0,0),GlobalVector(0,1,0));
339 
340  Surface::PositionType pos(0., 0.,400.5*fabs(etaParticle)/etaParticle);
341  PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos,rot);
342 
343 
344  TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane));
345  if(steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta())<3.0 )
346  endcapMC = steppingHelixstateinfo_.freeState()->position();
347  if(steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta())>3.0 )
348  doHF=kTRUE;
349  }
350  if(doHF) {
351 
352  if (abs(etaParticle)>5.191) continue;
353 
354  if(abs(etaParticle)>2.99) {
355  Surface::RotationType rot(GlobalVector(1,0,0),GlobalVector(0,1,0));
356 
357  Surface::PositionType pos1(0., 0.,1125*fabs(etaParticle)/etaParticle);
358  // Surface::PositionType pos1(0., 0.,1115*fabs(etaParticle)/etaParticle);
359  Surface::PositionType pos2(0., 0.,1137*fabs(etaParticle)/etaParticle);
360  PlaneBuilder::ReturnType aPlane1 = PlaneBuilder().plane(pos1,rot);
361  PlaneBuilder::ReturnType aPlane2 = PlaneBuilder().plane(pos2,rot);
362 
363 
364  TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane1));
365  if(steppingHelixstateinfo_.isValid() )
366  forwardMC1 = steppingHelixstateinfo_.freeState()->position();
367 
368  steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane2));
369  if(steppingHelixstateinfo_.isValid() )
370  forwardMC2 = steppingHelixstateinfo_.freeState()->position();
371  }
372  }
373  /* ------------ ------------ ----------------------------- */
374 
375 
376 
377  /* Finding the closest cell at Hcal */
378  Int_t iphitrue = -10;
379  Int_t ietatrue = 100;
380  HcalDetId tempId, tempId1, tempId2;
381 
382 
383  if (abs(etaParticle)<1.392)
384  {
385  gPointHcal = barrelMC;
386  tempId = gHcal->getClosestCell(gPointHcal);
387  }
388  if (abs(etaParticle)>=1.392 && abs(etaParticle)<3.0)
389  {
390  gPointHcal = endcapMC;
391  tempId = gHcal->getClosestCell(gPointHcal);
392  }
393  if (abs(etaParticle)>=3.0 && abs(etaParticle)<5.191)
394  {
395  /*
396  tempId1 = gHcal->getClosestCell(forwardMC1);
397  tempId2 = gHcal->getClosestCell(forwardMC2);
398  if (deltaR(tempId1.eta(), tempId1.phi(), etaParticle, phiParticle) < deltaR(tempId2.eta(), tempId2.phi(), etaParticle, phiParticle))
399  gPointHcal = forwardMC1;
400 
401  else gPointHcal = forwardMC2;
402  */
403  gPointHcal = forwardMC1;
404  tempId = gHcal->getClosestCell(gPointHcal);
405  //tempId = gHcal->CaloSubdetectorGeometry::getClosestCell(gPointHcal);
406  }
407 
408 
409 
410  tempId = gHcal->getClosestCell(gPointHcal);
411 
412  ietatrue = tempId.ieta();
413  iphitrue = tempId.iphi();
414 
415  etaGPoint = gPointHcal.eta();
416  phiGPoint = gPointHcal.phi();
417 
418  //xAtHcal = gPointHcal.x();
419  //yAtHcal = gPointHcal.y();
420  //zAtHcal = gPointHcal.z();
421  /* ----------------- ------------------------ */
422 
423 
424  if (gPointHcal.x()==0 && gPointHcal.y()==0 && gPointHcal.z()==0)
425  {/*cout <<"gPointHcal is Zero!"<<endl;*/ continue;}
426 
427 
428  float etahcal=gPointHcal.eta();
429  // float phihcal=gPointHcal.phi();
430  if (abs(etahcal)>5.192) continue;
431  //if (abs(etahcal)>3.0 && abs(etahcal)<5.191)
432 
433  //cout <<gPointHcal.x() <<" "<<gPointHcal.y() <<" "<<gPointHcal.z()<<" "<<gPointHcal.eta()<<" "<<gPointHcal.phi()<<" "<<ietatrue<<" "<<iphitrue <<endl;
434 
435  // if (ietatrue==100 || iphitrue==-10) {cout<<"ietatrue: "<<ietatrue<<" iphitrue: "<<iphitrue<<" etahcal: "<<etahcal<<" phihcal: "<<phihcal<<endl;}
436 
437 
438 
439 
440  /* ------------- Calculate Ecal Energy using TrackAssociator ---------------------- */
441 
442  //float etaecal=info.trkGlobPosAtEcal.eta();
443 
444  xTrkEcal=info.trkGlobPosAtEcal.x();
445  yTrkEcal=info.trkGlobPosAtEcal.y();
446  zTrkEcal=info.trkGlobPosAtEcal.z();
447 
448  GlobalPoint gPointEcal(xTrkEcal,yTrkEcal,zTrkEcal);
449 
450  eECAL = ecalEnergyInCone(gPointEcal, ecalCone_, Hitecal, geo);
451  eECAL09cm = ecalEnergyInCone(gPointEcal, 9, Hitecal, geo);
452  eECAL40cm = ecalEnergyInCone(gPointEcal, 40, Hitecal, geo);
453 
454  eEcalCone = eECAL;
455  //if(abs(etaecal)<1.5) enEcalB -> Fill(eEcalCone);
456  //if(abs(etaecal)>1.5 && abs(etaecal)<3.1) enEcalE -> Fill(eEcalCone);
457 
458  /* ------------------------------ -------------------------- ----------------- */
459 
460 
461  /* ----------------- Find the Hottest Hcal RecHit --------------------------- */
462  MaxHit_struct MaxHit;
463 
464  MaxHit.hitenergy=-100.;
465  Float_t recal = 1.0;
466 
467  //Hcal:
468  eHcalCone = 0.;
469  eHcalConeNoise = 0.;
470  UsedCells = 0;
471  UsedCellsNoise = 0;
472  e3x3 = 0.;
473  e5x5 = 0.;
474 
475  for (HBHERecHitCollection::const_iterator hhit=Hithbhe.begin(); hhit!=Hithbhe.end(); hhit++)
476  //for (HcalRecHitCollection::const_iterator hhit=Hithcal.begin(); hhit!=Hithcal.end(); hhit++)
477  {
478  recal = RecalibFactor(hhit->detid());
479  GlobalPoint pos = geo->getPosition(hhit->detid());
480 
481  int iphihit = (hhit->id()).iphi();
482  int ietahit = (hhit->id()).ieta();
483  int depthhit = (hhit->id()).depth();
484  float enehit = hhit->energy()* recal;
485 
486  if (depthhit!=1) continue;
487 
488  double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
489 
490  if(distAtHcal < clusterConeSize_)
491  {
492  for (HBHERecHitCollection::const_iterator hhit2=Hithbhe.begin(); hhit2!=Hithbhe.end(); hhit2++)
493  //for (HcalRecHitCollection::const_iterator hhit2=Hithcal.begin(); hhit2!=Hithcal.end(); hhit2++)
494  {
495  int iphihit2 = (hhit2->id()).iphi();
496  int ietahit2 = (hhit2->id()).ieta();
497  int depthhit2 = (hhit2->id()).depth();
498  float enehit2 = hhit2->energy() * recal;
499 
500  if (iphihit==iphihit2 && ietahit==ietahit2 && depthhit!=depthhit2) enehit = enehit+enehit2;
501 
502  }
503 
504 
505  if(enehit > MaxHit.hitenergy)
506  {
507  MaxHit.hitenergy = enehit;
508  MaxHit.ietahitm = (hhit->id()).ieta();
509  MaxHit.iphihitm = (hhit->id()).iphi();
510  MaxHit.dr = distAtHcal;
511  //MaxHit.depthhit = (hhit->id()).depth();
512  MaxHit.depthhit = 1;
513  }
514  }
515  }
516 
517 
518  for (HFRecHitCollection::const_iterator hhit=Hithf.begin(); hhit!=Hithf.end(); hhit++)
519  {
520 
521  recal = RecalibFactor(hhit->detid());
522 
523  GlobalPoint pos = geo->getPosition(hhit->detid());
524 
525  int iphihit = (hhit->id()).iphi();
526  int ietahit = (hhit->id()).ieta();
527  int depthhit = (hhit->id()).depth();
528  float enehit = hhit->energy()* recal;
529 
530  double distAtHcal = getDistInPlaneSimple(gPointHcal,pos);
531 
532  if(distAtHcal < associationConeSize_)
533  {
534  for (HFRecHitCollection::const_iterator hhit2=Hithf.begin(); hhit2!=Hithf.end(); hhit2++)
535  {
536  int iphihit2 = (hhit2->id()).iphi();
537  int ietahit2 = (hhit2->id()).ieta();
538  int depthhit2 = (hhit2->id()).depth();
539  float enehit2 = hhit2->energy() * recal;
540 
541  if (iphihit==iphihit2 && ietahit==ietahit2 && depthhit!=depthhit2) enehit = enehit+enehit2;
542 
543  }
544 
545  if(enehit > MaxHit.hitenergy)
546  {
547  MaxHit.hitenergy = enehit;
548  MaxHit.ietahitm = (hhit->id()).ieta();
549  MaxHit.iphihitm = (hhit->id()).iphi();
550  MaxHit.dr = distAtHcal;
551  MaxHit.depthhit = 1;
552  }
553  }
554  }
555 
556  /* ---------------------- ---------------- -------------------------------------------- */
557 
558 
559  /* ----------- Collect Hcal Energy in a Cone (and also 3x3 and 5x5 around the hottest cell)*/
560 
561  for (HBHERecHitCollection::const_iterator hhit=Hithbhe.begin(); hhit!=Hithbhe.end(); hhit++)
562  // for (HcalRecHitCollection::const_iterator hhit=Hithcal.begin(); hhit!=Hithcal.end(); hhit++)
563  {
564 
565  recal = RecalibFactor(hhit->detid());
566  //cout<<"recal: "<<recal<<endl;
567 
568  GlobalPoint pos = geo->getPosition(hhit->detid());
569 
570  int iphihit = (hhit->id()).iphi();
571  int ietahit = (hhit->id()).ieta();
572  int depthhit = (hhit->id()).depth();
573  float enehit = hhit->energy()* recal;
574 
575  //if (depthhit!=1) continue;
576 
577  //Set noise RecHit opposite to track hits
578  int iphihitNoise = iphihit >36 ? iphihit-36 : iphihit+36;
579  int ietahitNoise = -ietahit;
580  int depthhitNoise = depthhit;
581 
582  double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
583  if(distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
584  {
585  eHcalCone += enehit;
586  UsedCells++;
587 
588 
589  int DIETA = 100;
590  if(MaxHit.ietahitm*(hhit->id()).ieta()>0)
591  { DIETA = MaxHit.ietahitm - (hhit->id()).ieta();}
592  if(MaxHit.ietahitm*(hhit->id()).ieta()<0)
593  { DIETA = MaxHit.ietahitm - (hhit->id()).ieta(); DIETA = DIETA>0 ? DIETA-1 : DIETA+1;}
594  int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
595  DIPHI = DIPHI>36 ? 72-DIPHI : DIPHI;
596 
597 
598  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) )) )
599  {e5x5 += hhit->energy();}
600  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) )) )
601  {e3x3 += hhit->energy();}
602 
603  // cout<<"track: ieta "<<ietahit<<" iphi: "<<iphihit<<" depth: "<<depthhit<<" energydepos: "<<enehit<<endl;
604 
605  for (HBHERecHitCollection::const_iterator hhit2=Hithbhe.begin(); hhit2!=Hithbhe.end(); hhit2++)
606  {
607  recal = RecalibFactor(hhit2->detid());
608  int iphihit2 = (hhit2->id()).iphi();
609  int ietahit2 = (hhit2->id()).ieta();
610  int depthhit2 = (hhit2->id()).depth();
611  float enehit2 = hhit2->energy()* recal;
612 
613  if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2>0.)
614  {
615  eHcalConeNoise += hhit2->energy()*recal;
616  UsedCellsNoise++;
617  //cout<<"Noise: ieta "<<ietahit2<<" iphi: "<<iphihit2<<" depth: "<<depthhit2<<" energydepos: "<<enehit2<<endl;
618  }
619  }
620  }
621  } //end of all HBHE hits loop
622 
623 
624  for (HFRecHitCollection::const_iterator hhit=Hithf.begin(); hhit!=Hithf.end(); hhit++)
625  {
626 
627  recal = RecalibFactor(hhit->detid());
628 
629  GlobalPoint pos = geo->getPosition(hhit->detid());
630  //float phihit = pos.phi();
631  //float etahit = pos.eta();
632 
633  int iphihit = (hhit->id()).iphi();
634  int ietahit = (hhit->id()).ieta();
635  int depthhit = (hhit->id()).depth();
636  float enehit = hhit->energy()* recal;
637 
638  //Set noise RecHit opposite to track hits
639  int iphihitNoise = iphihit >36 ? iphihit-36 : iphihit+36;
640  int ietahitNoise = -ietahit;
641  int depthhitNoise = depthhit;
642 
643  double distAtHcal = getDistInPlaneSimple(gPointHcal,pos);
644 
645  if(distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
646  //if(dr<radius_ && enehit>0.)
647  {
648 
649  eHcalCone += enehit;
650  UsedCells++;
651 
652  int DIETA = 100;
653  if(MaxHit.ietahitm*(hhit->id()).ieta()>0)
654  { DIETA = MaxHit.ietahitm - (hhit->id()).ieta();}
655  if(MaxHit.ietahitm*(hhit->id()).ieta()<0)
656  { DIETA = MaxHit.ietahitm - (hhit->id()).ieta(); DIETA = DIETA>0 ? DIETA-1 : DIETA+1;}
657  int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
658  DIPHI = DIPHI>36 ? 72-DIPHI : DIPHI;
659 
660 
661  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) )) )
662  {e5x5 += hhit->energy();}
663  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) )) )
664  {e3x3 += hhit->energy();}
665 
666  for (HFRecHitCollection::const_iterator hhit2=Hithf.begin(); hhit2!=Hithf.end(); hhit2++)
667  {
668  recal = RecalibFactor(hhit2->detid());
669 
670  int iphihit2 = (hhit2->id()).iphi();
671  int ietahit2 = (hhit2->id()).ieta();
672  int depthhit2 = (hhit2->id()).depth();
673  float enehit2 = hhit2->energy()* recal;
674 
675  if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2>0.01)
676  {
677  eHcalConeNoise += hhit2->energy()*recal;
678  UsedCellsNoise++;
679  }
680  }
681 
682  }
683  } //end of all HF hits loop
684 
685  /* ---------------- -------------------- ---------------------------------------------- -------- */
686 
687 
688  /* ------------- - Track-MC matching (if any tracks are in event) ------------ - */
689 
690  nTracks=0;
691 
692  delRmc[0] = 5;
693 
694  float delR_track_particle = 100;
695 
696  for (reco::TrackCollection::const_iterator track1=generalTracks->begin(); track1!=generalTracks->end(); track1++)
697  {
698  delR_track_particle = deltaR(etaParticle, phiParticle, track1->eta(), track1->phi());
699 
700  trackEta[nTracks] = track1 -> eta();
701  trackPhi[nTracks] = track1 -> phi();
702  trackP[nTracks] = sqrt(track1->px()*track1->px() + track1->py()*track1->py() + track1->pz()*track1->pz());
703 
704  delRmc[nTracks] = delR_track_particle;
705  numValidTrkHits[nTracks] = track1->hitPattern().numberOfValidHits();
706  numValidTrkStrips[nTracks] = track1->hitPattern().numberOfValidStripTECHits();
707  numLayers[nTracks] = track1->hitPattern().trackerLayersWithMeasurement(); //layers crossed
708  trkQual[nTracks] = track1->quality(reco::TrackBase::highPurity);
709 
710  nTracks++;
711  }
712 
713 
714  /* ------------------ ------------------------------ ------- */
715 
716 
717  int dieta_M_P = 100;
718  int diphi_M_P = 100;
719  if(MaxHit.ietahitm*ietatrue>0) {dieta_M_P = abs (MaxHit.ietahitm-ietatrue);}
720  if(MaxHit.ietahitm*ietatrue<0) {dieta_M_P = abs(MaxHit.ietahitm-ietatrue)-1;}
721  diphi_M_P = abs(MaxHit.iphihitm-iphitrue);
722 
723  diphi_M_P = diphi_M_P>36 ? 72-diphi_M_P : diphi_M_P;
724  iDr = sqrt(diphi_M_P*diphi_M_P+dieta_M_P*dieta_M_P);
725 
726  /* if (iDr>15)
727  {
728 cout<<"diphi: "<<diphi_M_P<<" dieta: "<<dieta_M_P<<" iDr: "<<iDr<<" ietatrue:"<<ietatrue<<" iphitrue:"<<iphitrue<<endl;
729 cout<<"M ieta: "<<MaxHit.ietahitm<<" M iphi: "<<MaxHit.iphihitm<<endl;
730 
731 }*/
732 
733 
734  Bool_t passCuts = kFALSE;
735  passCuts=kTRUE;
736  //if(eEcalCone < energyECALmip && iDr<2.) passCuts = kTRUE;
737  //if(MaxHit.hitenergy>0.) passCuts = kTRUE;
738 
739 
740  if(passCuts)
741  {
742  /*
743  enHcal -> Fill(ietatrue, eHcalCone);
744  nCells -> Fill(ietatrue, UsedCells);
745  enHcalNoise -> Fill(ietatrue, eHcalConeNoise);
746  nCellsNoise -> Fill(ietatrue, UsedCellsNoise);
747  */
748 
749 
750  //e3x3=0; e5x5=0;
751 
752  iEta = ietatrue;
753  iPhi = iphitrue;
754 
755  //iEta = MaxHit.ietahitm;
756  //iPhi = MaxHit.iphihitm;
757  delR = MaxHit.dr;
758  eCentHit = MaxHit.hitenergy;
759 
760  eParticle = mom_MC;
761  //eTrack = mom_MC;
762  //phiTrack = phiParticle;
763  //etaTrack = etaParticle;
764 
765  pfTree->Fill();
766  }
767 
768  } //Hep:MC
769 
770 
771 }
772 
773 void HcalCorrPFCalculation::beginJob(){
774 
775  pfTree = fs -> make<TTree>("pfTree", "Tree for pf info");
776 
777 
778  pfTree->Branch("nTracks", &nTracks, "nTracks/I");
779  pfTree->Branch("trackEta", trackEta, "trackEta[nTracks]/F");
780  pfTree->Branch("trackPhi", trackPhi, "trackPhi[nTracks]/F");
781  pfTree->Branch("trackP", trackP, "trackP[nTracks]/F");
782 
783  pfTree->Branch("delRmc", delRmc, "delRmc[nTracks]/F");
784  pfTree->Branch("numValidTrkHits", numValidTrkHits, "numValidTrkHits[nTracks]/I");
785  pfTree->Branch("numValidTrkStrips", numValidTrkStrips, "numValidTrkStrips[nTracks]/I");
786  pfTree->Branch("numLayers", numLayers, "numLayers[nTracks]/I");
787  pfTree->Branch("trkQual", trkQual, "trkQual[nTracks]/O");
788 
789 
790  pfTree->Branch("eEcalCone", &eEcalCone, "eEcalCone/F");
791  pfTree->Branch("eHcalCone", &eHcalCone, "eHcalCone/F");
792  pfTree->Branch("eHcalConeNoise", &eHcalConeNoise, "eHcalConeNoise/F");
793 
794  pfTree->Branch("UsedCellsNoise", &UsedCellsNoise, "UsedCellsNoise/I");
795  pfTree->Branch("UsedCells", &UsedCells, "UsedCells/I");
796 
797  pfTree->Branch("eCentHit", &eCentHit , "eCentHit/F");
798 
799  pfTree->Branch("eParticle", &eParticle, "eParticle/F");
800  pfTree->Branch("etaParticle", &etaParticle, "etaParticle/F");
801  pfTree->Branch("phiParticle", &phiParticle, "phiParticle/F");
802 
803  pfTree->Branch("etaGPoint", &etaGPoint, "etaGPoint/F");
804  pfTree->Branch("phiGPoint", &phiGPoint, "phiGPoint/F");
805 
806  pfTree->Branch("xAtHcal", &xAtHcal, "xAtHcal/F");
807  pfTree->Branch("yAtHcal", &yAtHcal, "yAtHcal/F");
808  pfTree->Branch("zAtHcal", &zAtHcal, "zAtHcal/F");
809 
810  pfTree->Branch("eECAL09cm", &eECAL09cm, "eECAL09cm/F");
811  pfTree->Branch("eECAL40cm", &eECAL40cm, "eECAL40cm/F");
812  pfTree->Branch("eECAL", &eECAL, "eECAL/F");
813 
814  pfTree->Branch("e3x3 ", &e3x3 , "e3x3/F");
815  pfTree->Branch("e5x5", &e5x5 , "e5x5/F");
816 
817  pfTree->Branch("iDr", &iDr, "iDr/F");
818  pfTree->Branch("delR", &delR, "delR/F");
819 
820  pfTree->Branch("iEta", &iEta, "iEta/I");
821  pfTree->Branch("iPhi", &iPhi, "iPhi/I");
822 
823  // pfTree->Branch("numValidTrkHits", &numValidTrkHits, "numValidTrkHits/I");
824  // pfTree->Branch("numValidTrkStrips", &numValidTrkStrips, "numValidTrkStrips/I");
825  // pfTree->Branch("trkQual", &trkQual, "trkQual/");
826  // pfTree->Branch("numLayers", &numLayers, "numLayers/I");
827 
828 }
829 
830 
831 void HcalCorrPFCalculation::endJob()
832 {}
833 
834 
835 #include "FWCore/PluginManager/interface/ModuleDef.h"
836 #include "FWCore/Framework/interface/MakerMacros.h"
837 
838 
839 //define this as a plug-in
840 DEFINE_FWK_MODULE(HcalCorrPFCalculation);
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HcalCorrPFCalculation::tok_hbhe_
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
Definition: HcalCorrPFCalculation.cc:107
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
HcalCorrPFCalculation::geoHcal
const HcalGeometry * geoHcal
Definition: HcalCorrPFCalculation.cc:85
Parameters.parameters
dictionary parameters
Definition: Parameters.py:2
HcalCorrPFCalculation::pfRecalib
const HcalPFCorrs * pfRecalib
Definition: HcalCorrPFCalculation.cc:76
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
ecalEnergyInCone
double ecalEnergyInCone(const GlobalPoint center, double radius, const EcalRecHitCollection ecalCol, const CaloGeometry *geo)
Definition: CommonUsefulStuff.h:220
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
HcalCorrPFCalculation::beginJob
virtual void beginJob() override
Definition: HcalCorrPFCalculation.cc:773
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
PlaneBuilder::plane
ReturnType plane(const PositionType &pos, const RotationType &rot) const
Definition: PlaneBuilder.h:22
HcalCorrPFCalculation::parameters_
TrackAssociatorParameters parameters_
Definition: HcalCorrPFCalculation.cc:82
edm::SortedCollection< EcalRecHit >::const_iterator
std::vector< EcalRecHit >::const_iterator const_iterator
Definition: SortedCollection.h:82
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
npart
double npart
Definition: HydjetWrapper.h:49
ev
bool ev
Definition: Hydjet2Hadronizer.cc:95
MaxHit_struct::hitenergy
float hitenergy
Definition: CommonUsefulStuff.h:50
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:106
HLT_25ns14e33_v1_cff.InputTag
tuple InputTag
Definition: HLT_25ns14e33_v1_cff.py:45333
bk::beginJob
void beginJob()
Definition: Breakpoints.cc:15
MaxHit_struct::ietahitm
int ietahitm
Definition: CommonUsefulStuff.h:48
DataMixerPreMix_cff.generalTracks
tuple generalTracks
Definition: DataMixerPreMix_cff.py:16
HcalCorrPFCalculation::analyze
virtual void analyze(edm::Event const &ev, edm::EventSetup const &c) override
Definition: HcalCorrPFCalculation.cc:168
SteppingHelixPropagator_cfi.SteppingHelixPropagator
tuple SteppingHelixPropagator
Definition: SteppingHelixPropagator_cfi.py:3
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
TrajectoryStateOnSurface::freeState
FreeTrajectoryState const * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:76
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
HcalDetId::ieta
int ieta() const
get the cell ieta
Definition: HcalDetId.h:51
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ecalTB2006H4_GenSimDigiReco_cfg.bField
tuple bField
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:235
HcalCorrPFCalculation::geo
const CaloGeometry * geo
Definition: HcalCorrPFCalculation.cc:84
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:75
dbtoconf.conf
tuple conf
Definition: dbtoconf.py:185
HcalCorrPFCalculation::respRecalib
const HcalRespCorrs * respRecalib
Definition: HcalCorrPFCalculation.cc:75
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
CaloSubdetectorGeometry::getClosestCell
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: CaloSubdetectorGeometry.cc:56
HLT_25ns14e33_v1_cff.depth
tuple depth
Definition: HLT_25ns14e33_v1_cff.py:6940
HcalDetId::iphi
int iphi() const
get the cell iphi
Definition: HcalDetId.h:53
HcalCorrPFCalculation::endJob
virtual void endJob() override
Definition: HcalCorrPFCalculation.cc:831
FreeTrajectoryState::position
GlobalPoint position() const
Definition: FreeTrajectoryState.h:84
HcalCorrPFCalculation::stepPropF
SteppingHelixPropagator * stepPropF
Definition: HcalCorrPFCalculation.cc:78
MaxHit_struct::iphihitm
int iphihitm
Definition: CommonUsefulStuff.h:47
edm::Handle::product
T const * product() const
Definition: Handle.h:81
HcalCorrPFCalculation::tok_hf_
edm::EDGetTokenT< HFRecHitCollection > tok_hf_
Definition: HcalCorrPFCalculation.cc:109
HcalCorrPFCalculation::tok_EE_
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
Definition: HcalCorrPFCalculation.cc:111
HcalCorrPFCalculation::tok_EB_
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
Definition: HcalCorrPFCalculation.cc:112
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:56
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
phi
Geom::Phi< T > phi() const
Definition: Basic3DVectorLD.h:163
HcalCorrPFCalculation::tok_ho_
edm::EDGetTokenT< HORecHitCollection > tok_ho_
Definition: HcalCorrPFCalculation.cc:108
HcalCorrPFCalculation::trackAssociator_
TrackDetectorAssociator trackAssociator_
Definition: HcalCorrPFCalculation.cc:81
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:76
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
MaxHit_struct::dr
float dr
Definition: CommonUsefulStuff.h:51
HcalCorrPFCalculation::tok_tracks_
edm::EDGetTokenT< reco::TrackCollection > tok_tracks_
Definition: HcalCorrPFCalculation.cc:113
HcalCorrPFCalculation::RecalibFactor
double RecalibFactor(HcalDetId id)
Definition: HcalCorrPFCalculation.cc:154
MaxHit_struct::depthhit
int depthhit
Definition: CommonUsefulStuff.h:49
TrackDetMatchInfo::trkGlobPosAtEcal
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
Definition: TrackDetMatchInfo.h:39
getDistInPlaneSimple
double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint)
Definition: ConeDefinition.h:9
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
HcalCorrPFCalculation::fs
edm::Service< TFileService > fs
Definition: HcalCorrPFCalculation.cc:103
HcalCorrPFCalculation::tok_gen_
edm::EDGetTokenT< edm::HepMCProduct > tok_gen_
Definition: HcalCorrPFCalculation.cc:114
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
MaxHit_struct
Definition: CommonUsefulStuff.h:45
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:303
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
HcalCorrPFCalculation::HcalCorrPFCalculation
HcalCorrPFCalculation(edm::ParameterSet const &conf)
Definition: HcalCorrPFCalculation.cc:118
Generated for CMSSW Reference Manual by doxygen 1.8.5