CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ValidationHcalIsoTrackAlCaReco.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Calibration/HcalCalibAlgos/plugins
4 // Class: ValidationHcalIsoTrackAlCaReco
5 //
13 //
14 // Original Author: Grigory SAFRONOV, Sergey PETRUSHANKO
15 // Created: Tue Oct 14 16:10:31 CEST 2008
16 //
17 //
18 
19 
20 // system include files
21 #include <memory>
22 
23 // user include files
24 
26 
29 
32 
34 
38 
44 
47 
51 
55 
58 
59 
60 // Sergey +
61 
64 
65 // Sergey -
66 
68 #include <fstream>
69 
70 #include "TH1F.h"
71 
72 
73 double ValidationHcalIsoTrackAlCaReco::getDist(double eta1, double phi1, double eta2, double phi2)
74 {
75  double dphi = fabs(phi1 - phi2);
76  if(dphi>acos(-1)) dphi = 2*acos(-1)-dphi;
77  double dr = sqrt(dphi*dphi + pow(eta1-eta2,2));
78  return dr;
79 }
80 
81 // Sergey +
82 
83 double ValidationHcalIsoTrackAlCaReco::getDistInCM(double eta1, double phi1, double eta2, double phi2)
84 {
85  double dR, Rec;
86  double theta1=2*atan(exp(-eta1));
87  double theta2=2*atan(exp(-eta2));
88  if (fabs(eta1)<1.479) Rec=129;
89  else Rec=275;
90  //|vect| times tg of acos(scalar product)
91  dR=fabs((Rec/sin(theta1))*tan(acos(sin(theta1)*sin(theta2)*(sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2))+cos(theta1)*cos(theta2))));
92  return dR;
93 }
94 
95 // Sergey -
96 
97 std::pair<int,int> ValidationHcalIsoTrackAlCaReco::towerIndex(double eta, double phi)
98 {
99  int ieta=0;
100  int iphi=0;
101  for (int i=1; i<21; i++)
102  {
103  if (fabs(eta)<(i*0.087)&&fabs(eta)>(i-1)*0.087) ieta=int(fabs(eta)/eta)*i;
104  }
105  if (fabs(eta)>1.740&&fabs(eta)<1.830) ieta=int(fabs(eta)/eta)*21;
106  if (fabs(eta)>1.830&&fabs(eta)<1.930) ieta=int(fabs(eta)/eta)*22;
107  if (fabs(eta)>1.930&&fabs(eta)<2.043) ieta=int(fabs(eta)/eta)*23;
108 
109  double delta=phi+0.174532925;
110  if (delta<0) delta=delta+2*acos(-1);
111  if (fabs(eta)<1.740)
112  {
113  for (int i=0; i<72; i++)
114  {
115  if (delta<(i+1)*0.087266462&&delta>i*0.087266462) iphi=i;
116  }
117  }
118  else
119  {
120  for (int i=0; i<36; i++)
121  {
122  if (delta<2*(i+1)*0.087266462&&delta>2*i*0.087266462) iphi=2*i;
123  }
124  }
125 
126  return std::pair<int,int>(ieta,iphi);
127 }
128 
129 
131 
132  tok_simTrack_ = consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("simTracksTag"));
133 
134  folderName_ = iConfig.getParameter<std::string>("folderName");
135  saveToFile_=iConfig.getParameter<bool>("saveToFile");
136  outRootFileName_=iConfig.getParameter<std::string>("outputRootFileName");
137  tok_hlt_ = consumes<trigger::TriggerEvent>(iConfig.getParameter<edm::InputTag>("hltTriggerEventLabel"));
138  hltFilterTag_=iConfig.getParameter<edm::InputTag>("hltL3FilterLabel");
139  tok_arITr_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(iConfig.getParameter<edm::InputTag>("alcarecoIsoTracksLabel"));
140  recoTrLabel_=iConfig.getParameter<edm::InputTag>("recoTracksLabel");
141  pThr_=iConfig.getUntrackedParameter<double>("pThrL3",0);
142  heLow_=iConfig.getUntrackedParameter<double>("lowerHighEnergyCut",40);
143  heLow_=iConfig.getUntrackedParameter<double>("upperHighEnergyCut",60);
144 
145  nTotal=0;
146  nHLTL3accepts=0;
147 }
148 
150 {}
151 
153 {
154  nTotal++;
155 
157  iEvent.getByToken(tok_hlt_,trEv);
158 
160  iEvent.getByToken(tok_arITr_,recoIsoTracks);
161 
162  const trigger::TriggerObjectCollection& TOCol(trEv->getObjects());
163 
164  trigger::Keys KEYS;
165  const trigger::size_type nFilt(trEv->sizeFilters());
166  for (trigger::size_type iFilt=0; iFilt!=nFilt; iFilt++)
167  {
168  if (trEv->filterTag(iFilt)==hltFilterTag_)
169  {
170  KEYS=trEv->filterKeys(iFilt);
171  }
172  }
173 
174  trigger::size_type nReg=KEYS.size();
175 
176  std::vector<double> trigEta;
177  std::vector<double> trigPhi;
178  bool trig=false;
179 
180  //checks with IsoTrack trigger results
181  for (trigger::size_type iReg=0; iReg<nReg; iReg++)
182  {
183  const trigger::TriggerObject& TObj(TOCol[KEYS[iReg]]);
184  if (TObj.p()<pThr_) continue;
185  hl3eta->Fill(TObj.eta(),1);
186  hl3AbsEta->Fill(fabs(TObj.eta()),1);
187  hl3phi->Fill(TObj.phi(),1);
188 
189  if (recoIsoTracks->size()>0)
190  {
191  double minRecoL3dist=1000;
192  reco::IsolatedPixelTrackCandidateCollection::const_iterator mrtr;
193  for (reco::IsolatedPixelTrackCandidateCollection::const_iterator rtrit=recoIsoTracks->begin(); rtrit!=recoIsoTracks->end(); rtrit++)
194  {
195  double R=getDist(rtrit->eta(),rtrit->phi(),TObj.eta(),TObj.phi());
196  if (R<minRecoL3dist)
197  {
198  mrtr=rtrit;
199  minRecoL3dist=R;
200  }
201  }
202  hOffL3TrackMatch->Fill(minRecoL3dist,1);
203  hOffL3TrackPtRat->Fill(TObj.pt()/mrtr->pt(),1);
204  }
205 
206  hl3Pt->Fill(TObj.pt(),1);
207  trig=true;
208  trigEta.push_back(TObj.eta());
209  trigPhi.push_back(TObj.phi());
210  }
211 
212  //general distributions
213  for (reco::IsolatedPixelTrackCandidateCollection::const_iterator itr=recoIsoTracks->begin(); itr!=recoIsoTracks->end(); itr++)
214  {
215  bool match=false;
216  for (unsigned int l=0; l<trigEta.size(); l++)
217  {
218  if (getDist(itr->eta(),itr->phi(),trigEta[l],trigPhi[l])<0.2) match=true;
219  }
220  if (match&&trig)
221  {
222  hOffEtaFP->Fill(itr->eta(),1);
223  hOffPhiFP->Fill(itr->phi(),1);
224  }
225 
226  hOffEta->Fill(itr->eta(),1);
227  hOffPhi->Fill(itr->phi(),1);
228 
229  hOffAbsEta->Fill(fabs(itr->eta()),1);
230 
231  hDeposEcalInner->Fill(itr->energyIn(),1);
232  hDeposEcalOuter->Fill(itr->energyOut(),1);
233 
234  hTracksSumP->Fill(itr->sumPtPxl(),1);
235  hTracksMaxP->Fill(itr->maxPtPxl(),1);
236 
237  if (fabs(itr->eta())<0.5) hOffP_0005->Fill(itr->p(),1);
238  if (fabs(itr->eta())>0.5&&fabs(itr->eta())<1.0) hOffP_0510->Fill(itr->p(),1);
239  if (fabs(itr->eta())>1.0&&fabs(itr->eta())<1.5) hOffP_1015->Fill(itr->p(),1);
240  if (fabs(itr->eta())<1.5&&fabs(itr->eta())<2.0) hOffP_1520->Fill(itr->p(),1);
241 
242  hOffP->Fill(itr->p(),1);
243 
244  std::pair<int,int> TI=towerIndex(itr->eta(),itr->phi());
245  hOccupancyFull->Fill(TI.first,TI.second,1);
246  if (itr->p()>heLow_&&itr->p()<heUp_) hOccupancyHighEn->Fill(TI.first,TI.second,1);
247  }
248 
249 // Sergey +
250 
251  std::cout << std::endl << " End / Start " << std::endl;
252 
254  iEvent.getByToken<edm::SimTrackContainer>(tok_simTrack_, simTracks);
255 
256  for (reco::IsolatedPixelTrackCandidateCollection::const_iterator bll=recoIsoTracks->begin(); bll!=recoIsoTracks->end(); bll++)
257  {
258 
259  std::cout<<"ISO Pt " << bll->pt() << " P " << bll->p() << " Eta "<< bll->eta() << " Phi "<< bll->phi()<< std::endl;
260 
261  double distanceMin = 1.;
262  double SimPtMatched = 1.;
263  double SimPhiMatched = 1.;
264  double SimEtaMatched = 1.;
265  double SimDistMatched = 1.;
266  double SimPMatched = 1.;
267  double neuen = 0.;
268  double neuenm = 0.;
269  int neun = 0;
270 
271  for(edm::SimTrackContainer::const_iterator tracksCI = simTracks->begin();
272  tracksCI != simTracks->end(); tracksCI++){
273 
274  int partIndex = tracksCI->genpartIndex();
275  if (tracksCI->momentum().eta() > (bll->eta()-0.1)
276  && tracksCI->momentum().eta() < (bll->eta()+0.1)
277  && tracksCI->momentum().phi() > (bll->phi()-0.1)
278  && tracksCI->momentum().phi() < (bll->phi()+0.1)
279 // && tracksCI->momentum().e() > (0.5*bll->p())
280 // && tracksCI->momentum().e() < (2.*bll->p())
281  && tracksCI->momentum().e() > 2.
282  && fabs(tracksCI->charge()) == 1 && partIndex >0)
283  {
284 
285  double distance=getDist(tracksCI->momentum().eta(),tracksCI->momentum().phi(),bll->eta(),bll->phi());
286  double distanceCM=getDistInCM(tracksCI->momentum().eta(),tracksCI->momentum().phi(),bll->eta(),bll->phi());
287 
288  if (distanceMin > distance) {
289  distanceMin = distance;
290  SimPtMatched = tracksCI->momentum().pt();
291  SimPhiMatched = tracksCI->momentum().phi();
292  SimEtaMatched = tracksCI->momentum().eta();
293  SimDistMatched = distance;
294  SimPMatched = sqrt(tracksCI->momentum().pt()*tracksCI->momentum().pt() + tracksCI->momentum().pz()*tracksCI->momentum().pz());
295  }
296 
297  std::cout<<" Pt "<<tracksCI->momentum().pt()
298  << " Energy " << tracksCI->momentum().e()
299  << " Eta "<< tracksCI->momentum().eta()
300  << " Phi "<< tracksCI->momentum().phi()
301  << " Ind " << partIndex
302  << " Cha " << tracksCI->charge()
303  << " Dis " << distance
304  << " DCM " << distanceCM
305  << std::endl;
306 
307  }
308 
309 
310  if (
311  tracksCI->momentum().eta() > (bll->eta()-0.5)
312  && tracksCI->momentum().eta() < (bll->eta()+0.5)
313  && tracksCI->momentum().phi() > (bll->phi()-0.5)
314  && tracksCI->momentum().phi() < (bll->phi()+0.5)
315 // && tracksCI->momentum().e() > 2.
316  && tracksCI->charge() == 0 && partIndex >0)
317  {
318 
319  double distance=getDist(tracksCI->momentum().eta(),tracksCI->momentum().phi(),bll->eta(),bll->phi());
320  double distanceCM=getDistInCM(tracksCI->momentum().eta(),tracksCI->momentum().phi(),bll->eta(),bll->phi());
321 
322  std::cout<<"NEU Pt "<<tracksCI->momentum().pt()
323  << " Energy " << tracksCI->momentum().e()
324  << " Eta "<< tracksCI->momentum().eta()
325  << " Phi "<< tracksCI->momentum().phi()
326  << " Ind " << partIndex
327  << " Cha " << tracksCI->charge()
328  << " Dis " << distance
329  << " DCM " << distanceCM
330  << std::endl;
331 
332  if (distanceCM < 40.){
333 
334  neuen = neuen + tracksCI->momentum().e();
335  neun = neun + 1;
336  if (neuenm < tracksCI->momentum().e()) neuenm = tracksCI->momentum().e();
337 
338  }
339 
340  }
341 
342  }
343 
344  hSimNN->Fill(neun,1);
345  hSimNE->Fill(neuen,1);
346  hSimNM->Fill(neuenm,1);
347 
348  if (distanceMin < 0.1) {
349 
350  hSimPt->Fill(SimPtMatched,1);
351  hSimPhi->Fill(SimPhiMatched,1);
352  hSimEta->Fill(SimEtaMatched,1);
353  hSimAbsEta->Fill(fabs(SimEtaMatched),1);
354  hSimDist->Fill(SimDistMatched,1);
355  hSimPtRatOff->Fill(SimPtMatched/bll->pt(),1);
356  hSimP->Fill(SimPMatched,1);
357  hSimN->Fill(1,1);
358 
359  std::cout<<"S Pt "<< SimPtMatched
360  << std::endl;
361  }
362 
363  if (distanceMin > 0.1) {
364  hSimN->Fill(0,1);
365  }
366 
367 
368 
369 
370  }
371 
372 // Sergey -
373 
374 
375 
376 }
377 
379 {
382 
383  hl3Pt=dbe_->book1D("hl3Pt","pT of hlt L3 objects",1000,0,1000);
384  hl3Pt->setAxisTitle("pT(GeV)",1);
385  hl3eta=dbe_->book1D("hl3eta","eta of hlt L3 objects",16,-2,2);
386  hl3eta->setAxisTitle("eta",1);
387  hl3AbsEta=dbe_->book1D("hl3AbsEta","|eta| of hlt L3 objects",8,0,2);
388  hl3AbsEta->setAxisTitle("eta",1);
389  hl3phi=dbe_->book1D("hl3phi","phi of hlt L3 objects",16,-3.2,3.2);
390  hl3phi->setAxisTitle("phi",1);
391 
392  hOffEta=dbe_->book1D("hOffEta","eta of alcareco objects",100,-2,2);
393  hOffEta->setAxisTitle("eta",1);
394  hOffPhi=dbe_->book1D("hOffPhi","phi of alcareco objects",100,-3.2,3.2);
395  hOffPhi->setAxisTitle("phi",1);
396  hOffP=dbe_->book1D("hOffP","p of alcareco objects",1000,0,1000);
397  hOffP->setAxisTitle("E(GeV)",1);
398  hOffP_0005=dbe_->book1D("hOffP_0005","p of alcareco objects, |eta|<0.5",1000,0,1000);
399  hOffP_0005->setAxisTitle("E(GeV)",1);
400  hOffP_0510=dbe_->book1D("hOffP_0510","p of alcareco objects, 0.5<|eta|<1.0",1000,0,1000);
401  hOffP_0510->setAxisTitle("E(GeV)",1);
402  hOffP_1015=dbe_->book1D("hOffP_1015","p of alcareco objects, 1.0<|eta|<1.5",1000,0,1000);
403  hOffP_1015->setAxisTitle("E(GeV)",1);
404  hOffP_1520=dbe_->book1D("hOffP_1520","p of alcareco objects, 1.5<|eta|<2.0",1000,0,1000);
405  hOffP_1520->setAxisTitle("E(GeV)",1);
406  hOffEtaFP=dbe_->book1D("hOffEtaFP","eta of alcareco objects, FP",16,-2,2);
407  hOffEtaFP->setAxisTitle("eta",1);
408  hOffAbsEta=dbe_->book1D("hOffAbsEta","|eta| of alcareco objects",8,0,2);
409  hOffAbsEta->setAxisTitle("|eta|",1);
410  hOffPhiFP=dbe_->book1D("hOffPhiFP","phi of alcareco objects, FP",16,-3.2,3.2);
411  hOffPhiFP->setAxisTitle("phi",1);
412  hTracksSumP=dbe_->book1D("hTracksSumP","summary p of tracks in the isolation cone",100,0,20);
413  hTracksSumP->setAxisTitle("E(GeV)");
414  hTracksMaxP=dbe_->book1D("hTracksMaxP","maximum p among tracks in the isolation cone",100,0,20);
415  hTracksMaxP->setAxisTitle("E(GeV)");
416 
417  hDeposEcalInner=dbe_->book1D("hDeposEcalInner","ecal energy deposition in inner cone around track",1000,0,1000);
418  hDeposEcalInner->setAxisTitle("E(GeV)");
419  hDeposEcalOuter=dbe_->book1D("hDeposEcalOuter","ecal energy deposition in outer cone around track",1000,0,1000);
420  hDeposEcalOuter->setAxisTitle("E(GeV)");
421 
422  hOccupancyFull=dbe_->book2D("hOccupancyFull","number of tracks per tower, full energy range",48,-25,25,73,0,73);
423  hOccupancyFull->setAxisTitle("ieta",1);
424  hOccupancyFull->setAxisTitle("iphi",2);
425  hOccupancyFull->getTH2F()->SetOption("colz");
426  hOccupancyFull->getTH2F()->SetStats(kFALSE);
427  hOccupancyHighEn=dbe_->book2D("hOccupancyHighEn","number of tracks per tower, high energy tracks",48,-25,25,73,0,73);
428  hOccupancyHighEn->setAxisTitle("ieta",1);
429  hOccupancyHighEn->setAxisTitle("iphi",2);
430  hOccupancyHighEn->getTH2F()->SetOption("colz");
431  hOccupancyHighEn->getTH2F()->SetStats(kFALSE);
432 
433  hOffL3TrackMatch=dbe_->book1D("hOffL3TrackMatch","Distance from L3 object to offline track",200,0,0.5);
434  hOffL3TrackMatch->setAxisTitle("R(eta,phi)",1);
435  hOffL3TrackPtRat=dbe_->book1D("hOffL3TrackPtRat","Ratio of pT: L3/offline",100,0,10);
436  hOffL3TrackPtRat->setAxisTitle("ratio L3/offline",1);
437 
438 
439 // Sergey +
440 
441  hSimPt=dbe_->book1D("hSimPt","pT of matched SimTrack",1000,0,1000);
442  hSimPt->setAxisTitle("pT(GeV)",1);
443 
444  hSimPhi=dbe_->book1D("hSimPhi","Phi of matched SimTrack",100,-3.2,3.2);
445  hSimPhi->setAxisTitle("phi",1);
446 
447  hSimEta=dbe_->book1D("hSimEta","Eta of matched SimTrack",100,-2.,2.);
448  hSimEta->setAxisTitle("eta",1);
449 
450  hSimAbsEta=dbe_->book1D("hSimAbsEta","|eta| of matched SimTrack",8,0.,2.);
451  hSimAbsEta->setAxisTitle("|eta|",1);
452 
453  hSimDist=dbe_->book1D("hSimDist","Distance from matched SimTrack to Offline Track",200,0,0.1);
454  hSimDist->setAxisTitle("R(eta,phi)",1);
455 
456  hSimPtRatOff=dbe_->book1D("hSimPtRatOff","pT Sim / pT Offline",100,0,10);
457  hSimPtRatOff->setAxisTitle("pT Sim / pT Offline",1);
458 
459  hSimP=dbe_->book1D("hSimP","p of matched SimTrack",1000,0,1000);
460  hSimP->setAxisTitle("p(GeV)",1);
461 
462  hSimN=dbe_->book1D("hSimN","Number matched",2,0,2);
463  hSimN->setAxisTitle("Offline/SimTRack - Matched or Not",1);
464 
465  hSimNN=dbe_->book1D("hSimNN","Number of the neutral particles in cone on ECAL",100,0,100);
466  hSimNN->setAxisTitle("Number",1);
467 
468  hSimNE=dbe_->book1D("hSimNE","Total energy of the neutral particles in cone on ECAL",100,0,100);
469  hSimNE->setAxisTitle("Energy",1);
470 
471  hSimNM=dbe_->book1D("hSimNM","Maximum energy of the neutral particles in cone on ECAL",100,0,100);
472  hSimNM->setAxisTitle("Energy",1);
473 
474 // Sergey -
475 
476 
477 }
478 
480 
481 if(dbe_)
482  {
484  }
485 }
486 
dbl * delta
Definition: mlp_gen.cc:36
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
virtual void analyze(const edm::Event &, const edm::EventSetup &)
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:954
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
float phi() const
Definition: TriggerObject.h:58
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double getDistInCM(double eta1, double phi1, double eta2, double phi2)
edm::EDGetTokenT< edm::SimTrackContainer > tok_simTrack_
double getDist(double, double, double, double)
float eta() const
Definition: TriggerObject.h:57
T eta() const
uint16_t size_type
edm::EDGetTokenT< reco::IsolatedPixelTrackCandidateCollection > tok_arITr_
void Fill(long long x)
Single trigger physics object (e.g., an isolated muon)
Definition: TriggerObject.h:22
int iEvent
Definition: GenABIO.cc:230
T sqrt(T t)
Definition: SSEVec.h:48
std::pair< int, int > towerIndex(double eta, double phi)
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
std::vector< TriggerObject > TriggerObjectCollection
collection of trigger physics objects (e.g., all isolated muons)
Definition: TriggerObject.h:81
std::vector< size_type > Keys
edm::EDGetTokenT< trigger::TriggerEvent > tok_hlt_
ValidationHcalIsoTrackAlCaReco(const edm::ParameterSet &)
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
tuple cout
Definition: gather_cfg.py:121
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
TH2F * getTH2F(void) const
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
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
std::vector< SimTrack > SimTrackContainer
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:667
Definition: DDAxes.h:10