CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ggPFClusters.cc
Go to the documentation of this file.
1 #include "RecoEgamma/EgammaTools/interface/ggPFClusters.h"
2 #include "DataFormats/EcalDetId/interface/EBDetId.h"
3 #include "DataFormats/EcalDetId/interface/EEDetId.h"
4 #include "DataFormats/Math/interface/deltaR.h"
5 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
6 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
7 #include "TMath.h"
8 /*
9 class by Rishi Patel rpatel@ cern.ch
10 */
11 
12 
13 ggPFClusters::ggPFClusters(
14  // reco::Photon PFPhoton,
15  edm::Handle<EcalRecHitCollection>& EBReducedRecHits,
16  edm::Handle<EcalRecHitCollection>& EEReducedRecHits,
17  const CaloSubdetectorGeometry* geomBar,
18  const CaloSubdetectorGeometry* geomEnd
19 
20  ):
21  EBReducedRecHits_(EBReducedRecHits),
22  EEReducedRecHits_(EEReducedRecHits),
23  geomBar_(geomBar),
24  geomEnd_(geomEnd)
25 {
26 
27 }
28 
29 
30 ggPFClusters::~ggPFClusters(){
31 
32 }
33 
34 
35 vector<reco::CaloCluster> ggPFClusters::getPFClusters(reco::SuperCluster sc){
36  vector<reco::CaloCluster> PFClust;
37  reco::CaloCluster_iterator it=sc.clustersBegin();
38  //get PFCluster Position from Basic Clusters, get Energy from RecHits
39  for(;it!=sc.clustersEnd();++it){
40  std::vector< std::pair<DetId, float> >bcCells=(*it)->hitsAndFractions();
41  DetId seedXtalId = bcCells[0].first ;
42  int detector = seedXtalId.subdetId(); //use Seed to check if EB or EE
43  bool isEb;
44  if(detector==1)isEb=true;
45  else isEb=false;
46 
47  float ClusSum=SumPFRecHits(bcCells, isEb); //return PFCluster Energy by Matching to RecHit and using the fractions from the bcCells
48  CaloCluster calo(ClusSum, (*it)->position());//store in CaloCluster (parent of PFCluster)
49  for(unsigned int i=0; i<bcCells.size(); ++i){
50  calo.addHitAndFraction(bcCells[i].first, bcCells[i].second);//Store DetIds and Fractions
51  }
52  PFClust.push_back(calo);//store in the Vector
53  }
54  return PFClust;
55 }
56 
57 float ggPFClusters::SumPFRecHits(std::vector< std::pair<DetId, float> >& bcCells, bool isEB){
58  float ClustSum=0;
59  for(unsigned int i=0; i<bcCells.size(); ++i){//loop over the Basic Cluster Cells
60  EcalRecHitCollection::const_iterator eb;
61  EcalRecHitCollection::const_iterator ee;
62 
63  if(isEB){
64  for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){//loop over RecHits
65  if(eb->id().rawId()==bcCells[i].first.rawId()){//match
66  float cellE=eb->energy()* bcCells[i].second;//Energy times fraction
67  ClustSum=ClustSum+cellE;
68  }
69  }
70  }
71  else{
72  for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){//loop over RecHits
73  if(ee->id().rawId()==bcCells[i].first.rawId()){//match
74  float cellE=ee->energy()* bcCells[i].second;//Energy times fraction
75  ClustSum=ClustSum+cellE;
76  break;
77  }
78  }
79  }
80  }
81  return ClustSum;
82 }
83 
84 float ggPFClusters::getPFSuperclusterOverlap(reco::CaloCluster PFClust, reco::SuperCluster sc){
85  float overlap=0;
86  reco::CaloCluster_iterator pit=sc.clustersBegin();
87  std::vector< std::pair<DetId, float> >bcCellsPF=PFClust.hitsAndFractions();
88 
89  DetId seedXtalId = bcCellsPF[0].first ;
90  int detector = seedXtalId.subdetId() ;
91  bool isEb;
92  std::vector< std::pair<DetId, float> >bcCellsreco;
93  if(detector==1)isEb=true;
94  else isEb=false;
95  for(;pit!=sc.clustersEnd();++pit){//fill vector of basic Clusters from SuperCluster
96  std::vector< std::pair<DetId, float> >bcCells2=(*pit)->hitsAndFractions();
97  for(unsigned int h=0; h<bcCells2.size(); ++h)bcCellsreco.push_back(bcCells2[h]);
98  }
99  float clustOverlap=0;
100  clustOverlap=PFRecHitsSCOverlap(//find overlap of a PFCluster with SuperCluster
101  bcCellsPF,
102  bcCellsreco,
103  isEb);
104  overlap=clustOverlap;
105  return overlap;
106 
107 
108 }
109 
110 
111 float ggPFClusters::getPFSuperclusterOverlap(reco::CaloCluster PFClust, reco::Photon phot){
112  float overlap=0;
113  SuperClusterRef recoSC=phot.superCluster();
114  reco::CaloCluster_iterator pit=recoSC->clustersBegin();
115 
116  std::vector< std::pair<DetId, float> >bcCellsPF=PFClust.hitsAndFractions();
117 
118  DetId seedXtalId = bcCellsPF[0].first ;
119  int detector = seedXtalId.subdetId() ;
120  bool isEb;
121  std::vector< std::pair<DetId, float> >bcCellsreco;
122  if(detector==1)isEb=true;
123  else isEb=false;
124  for(;pit!=recoSC->clustersEnd();++pit){//fill vector of basic Clusters from SuperCluster
125  std::vector< std::pair<DetId, float> >bcCells2=(*pit)->hitsAndFractions();
126  for(unsigned int h=0; h<bcCells2.size(); ++h)bcCellsreco.push_back(bcCells2[h]);
127  }
128  float clustOverlap=0;
129  clustOverlap=PFRecHitsSCOverlap(//find overlap of a PFCluster with SuperCluster
130  bcCellsPF,
131  bcCellsreco,
132  isEb);
133  overlap=clustOverlap;
134  return overlap;
135 }
136 
137 float ggPFClusters::PFRecHitsSCOverlap(
138  std::vector< std::pair<DetId, float> >& bcCells1,
139  std::vector< std::pair<DetId, float> >& bcCells2,
140  bool isEB){
141  float OverlapSum=0;
142  multimap <DetId, float> pfDID;
143  multimap <DetId, float> recoDID;
144  multimap<DetId, float>::iterator pfit;
145  multimap<DetId, float>::iterator pit;
146  vector<DetId>matched;
147  vector<float>matchedfrac;
148  //fill Multimap of DetId, fraction for PFCluster
149  for(unsigned int i=0; i<bcCells1.size(); ++i){
150  pfDID.insert(make_pair(bcCells1[i].first, bcCells1[i].second));
151  }
152  //fill Multimap of DetId, fraction for SuperCluster
153  for(unsigned int i=0; i<bcCells2.size(); ++i){
154  recoDID.insert(make_pair(bcCells2[i].first, bcCells2[i].second));
155  }
156  pit=recoDID.begin();
157  pfit=pfDID.begin();
158 
159  for(; pfit!=pfDID.end(); ++pfit){
160  pit=recoDID.find(pfit->first);//match DetId from PFCluster to Supercluster
161  if(pit!=recoDID.end()){
162  // cout<<"Found Match "<<endl;
163  matched.push_back(pfit->first);//store detId
164  matchedfrac.push_back(pfit->second);//store fraction
165  }
166  }
167 
168  for(unsigned int m=0; m<matched.size(); ++m){ //loop over matched cells
169  DetId det=matched[m];
170  EcalRecHitCollection::const_iterator eb;
171  EcalRecHitCollection::const_iterator ee;
172  if(isEB){
173  for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
174  if(eb->id().rawId()==det.rawId()){
175  float cellE=eb->energy()* matchedfrac[m]; //compute overlap
176  OverlapSum=OverlapSum+cellE;
177  break;
178  }
179  }
180  }
181  else{
182  for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
183  if(ee->id().rawId()==det.rawId()){
184  float cellE=ee->energy() * matchedfrac[m];//compute overlap
185  OverlapSum=OverlapSum+cellE;
186  break;
187  }
188  }
189  }
190 
191  }
192  return OverlapSum;
193 }
194 
195 void ggPFClusters::localCoordsEB( reco::CaloCluster clus, float &etacry, float &phicry, int &ieta, int &iphi, float &thetatilt, float &phitilt){
196  const math::XYZPoint position_ = clus.position();
197  double Theta = -position_.theta()+0.5*TMath::Pi();
198  double Eta = position_.eta();
199  double Phi = TVector2::Phi_mpi_pi(position_.phi());
200 
201  //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
202  // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
203  const float X0 = 0.89; const float T0 = 7.4;
204  double depth = X0 * (T0 + log(clus.energy()));
205 
206  //find max energy crystal
207  std::vector< std::pair<DetId, float> > crystals_vector = clus.hitsAndFractions();
208  float drmin = 999.;
209  EBDetId crystalseed(crystals_vector[0].first);
210  //printf("starting loop over crystals, etot = %5f:\n",clus.energy());
211  for (unsigned int icry=0; icry!=crystals_vector.size(); ++icry) {
212 
213  EBDetId crystal(crystals_vector[icry].first);
214 
215  const CaloCellGeometry* cell=geomBar_->getGeometry(crystal);
216  GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
217  double EtaCentr = center_pos.eta();
218  double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
219 
220  float dr = reco::deltaR(Eta,Phi,EtaCentr,PhiCentr);
221  if (dr<drmin) {
222  drmin = dr;
223  crystalseed = crystal;
224  }
225 
226  }
227 
228  ieta = crystalseed.ieta();
229  iphi = crystalseed.iphi();
230 
231  // Get center cell position from shower depth
232  const CaloCellGeometry* cell=geomBar_->getGeometry(crystalseed);
233  const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid*>(cell);
234 
235  thetatilt = cpyr->getThetaAxis();
236  phitilt = cpyr->getPhiAxis();
237 
238  GlobalPoint center_pos = cpyr->getPosition(depth);
239 
240  double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
241  double PhiWidth = (TMath::Pi()/180.);
242  phicry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth;
243  //Some flips to take into account ECAL barrel symmetries:
244  if (ieta<0) phicry *= -1.;
245 
246  double ThetaCentr = -center_pos.theta()+0.5*TMath::Pi();
247  double ThetaWidth = (TMath::Pi()/180.)*TMath::Cos(ThetaCentr);
248  etacry = (Theta-ThetaCentr)/ThetaWidth;
249  //cout<<"eta, phi raw w/o widths "<<(Theta-ThetaCentr)<<", "<<(TVector2::Phi_mpi_pi(Phi-PhiCentr))<<endl;
250  //flip to take into account ECAL barrel symmetries:
251  if (ieta<0) etacry *= -1.;
252  return;
253 
254 }
255 
256 void ggPFClusters::localCoordsEE(reco::CaloCluster clus, float &xcry, float &ycry, int &ix, int &iy, float &thetatilt, float &phitilt){
257  const math::XYZPoint position_ = clus.position();
258  //double Theta = -position_.theta()+0.5*TMath::Pi();
259  double Eta = position_.eta();
260  double Phi = TVector2::Phi_mpi_pi(position_.phi());
261  double X = position_.x();
262  double Y = position_.y();
263 
264  //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
265  // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
266  const float X0 = 0.89; float T0 = 1.2;
267  //different T0 value if outside of preshower coverage
268  if (TMath::Abs(clus.eta())<1.653) T0 = 3.1;
269 
270  double depth = X0 * (T0 + log(clus.energy()));
271 
272  //find max energy crystal
273  std::vector< std::pair<DetId, float> > crystals_vector = clus.hitsAndFractions();
274  float drmin = 999.;
275  EEDetId crystalseed(crystals_vector[0].first);
276  //printf("starting loop over crystals, etot = %5f:\n",bclus.energy());
277  for (unsigned int icry=0; icry!=crystals_vector.size(); ++icry) {
278 
279  EEDetId crystal(crystals_vector[icry].first);
280 
281  const CaloCellGeometry* cell=geomEnd_->getGeometry(crystal);
282  GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
283  double EtaCentr = center_pos.eta();
284  double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
285 
286  float dr = reco::deltaR(Eta,Phi,EtaCentr,PhiCentr);
287  if (dr<drmin) {
288  drmin = dr;
289  crystalseed = crystal;
290  }
291 
292  }
293 
294  ix = crystalseed.ix();
295  iy = crystalseed.iy();
296 
297  // Get center cell position from shower depth
298  const CaloCellGeometry* cell=geomEnd_->getGeometry(crystalseed);
299  const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid*>(cell);
300 
301  thetatilt = cpyr->getThetaAxis();
302  phitilt = cpyr->getPhiAxis();
303 
304  GlobalPoint center_pos = cpyr->getPosition(depth);
305 
306  double XCentr = center_pos.x();
307  double XWidth = 2.59;
308  xcry = (X-XCentr)/XWidth;
309 
310 
311  double YCentr = center_pos.y();
312  double YWidth = 2.59;
313  ycry = (Y-YCentr)/YWidth;
314  //cout<<"x, y raw w/o widths "<<(X-XCentr)<<", "<<(Y-YCentr)<<endl;
315 }
316 
317 float ggPFClusters::get5x5Element(int i, int j,
318  std::vector< std::pair<DetId, float> >&
319  bcCells,
320  bool isEB){
321 
322  Fill5x5Map(bcCells,isEB);
323  if(abs(i)>2 ||abs(j)>2)return 0; //outside 5x5
324  int ind1=i+2;
325  int ind2=j+2;
326  float E=e5x5_[ind1][ind2]; //return element from 5x5 cells
327  return E;
328 }
329 
330 void ggPFClusters::Fill5x5Map(std::vector< std::pair<DetId, float> >& bcCells,
331  bool isEB){
332 
333  for(int i=0; i<5; ++i)
334  for(int j=0; j<5; ++j)e5x5_[i][j]=0;
335  //cout<<"here 2 "<<endl;
336  EcalRecHitCollection::const_iterator eb;
337  EcalRecHitCollection::const_iterator ee;
338 
339  DetId idseed=FindSeed(bcCells, isEB);//return seed crystal
340 
341  for(unsigned int i=0; i<bcCells.size(); ++i){
342  DetId id=bcCells[i].first;
343  if(isEB){
344  EBDetId EBidSeed=EBDetId(idseed.rawId());
345  int deta=EBDetId::distanceEta(id,idseed);
346  int dphi=EBDetId::distancePhi(id,idseed);
347  if(abs(dphi)<=2 && abs(deta)<=2){
348  for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
349  if(eb->id().rawId()==bcCells[i].first.rawId()){
350  EBDetId EBid=EBDetId(id.rawId());
351  int ind1=EBidSeed.ieta()-EBid.ieta();
352  int ind2=EBid.iphi()-EBidSeed.iphi();
353  if(EBidSeed.ieta() * EBid.ieta() > 0){
354  ind1=EBid.ieta()-EBidSeed.ieta();
355  }
356  else{ //near EB+ EB-
357  //bugged
358  //ind1=(1-(EBidSeed.ieta()-EBid.ieta()));
359  ind1 = (EBid.ieta()-EBidSeed.ieta())*(abs(EBid.ieta()-EBidSeed.ieta())-1)/abs(EBid.ieta()-EBidSeed.ieta());
360  }
361  int iEta=ind1+2;
362  int iPhi=ind2+2;
363  e5x5_[iEta][iPhi]=eb->energy()* bcCells[i].second;
364  }
365  }
366  }
367  }
368  else{
369  EEDetId EEidSeed=EEDetId(idseed.rawId());
370  int dx=EEDetId::distanceX(id,idseed);
371  int dy=EEDetId::distanceY(id,idseed);
372  if(abs(dx)<=2 && abs(dy)<=2){
373  for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
374  if(ee->id().rawId()==bcCells[i].first.rawId()){
375  EEDetId EEid=EEDetId(id.rawId());
376  int ind1=EEid.ix()-EEidSeed.ix();
377  int ind2=EEid.iy()-EEidSeed.iy();
378  int ix=ind1+2;
379  int iy=ind2+2;
380  e5x5_[ix][iy]=ee->energy()* bcCells[i].second;
381  }
382  }
383  }
384  }
385  }
386 }
387 
388 DetId ggPFClusters::FindSeed(std::vector< std::pair<DetId, float> >& bcCells, bool isEB){
389  //first find seed:
390  EcalRecHitCollection::const_iterator eb;
391  EcalRecHitCollection::const_iterator ee;
392  DetId idseed=bcCells[0].first;
393  float PFSeedE=0;
394  //find seed by largest energy matching
395  for(unsigned int i=0; i<bcCells.size(); ++i){
396  if(isEB){
397  for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
398  if(eb->id().rawId()==bcCells[i].first.rawId()){
399  float cellE=eb->energy()* bcCells[i].second;
400  if(PFSeedE<cellE){
401  PFSeedE=cellE;
402  idseed=bcCells[i].first;
403  }
404  break;
405  }
406  }
407  }
408  else{
409  for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
410 
411  if(ee->id().rawId()==bcCells[i].first.rawId()){
412  float cellE=ee->energy()* bcCells[i].second;
413  if(PFSeedE<cellE){
414  PFSeedE=cellE;
415  idseed=bcCells[i].first;
416  }
417  break;
418  }
419  }
420  }
421  }
422  return idseed;
423 }
424 
425 std::pair<float, float> ggPFClusters::ClusterWidth(
426  vector<reco::CaloCluster>& PFClust){
427  pair<float, float> widths(0,0);
428  multimap<float, unsigned int>ClusterMap;
429  //fill in Multimap to order by Energy (ascending order
430  for(unsigned int i=0; i<PFClust.size();++i)ClusterMap.insert(make_pair(PFClust[i].energy(), i));
431  //reverse iterator to get cluster with largest first
432  std::multimap<float, unsigned int>::reverse_iterator it;
433  it=ClusterMap.rbegin();
434  unsigned int max_c=(*it).second;
435  std::vector< std::pair<DetId, float> >bcCells=PFClust[max_c].hitsAndFractions();
436  EcalRecHitCollection::const_iterator eb;
437  EcalRecHitCollection::const_iterator ee;
438  float numeratorEtaWidth=0;
439  float numeratorPhiWidth=0;
440  float denominator=0;
441 
442  DetId seedXtalId = bcCells[0].first ;
443  int detector = seedXtalId.subdetId() ;
444  bool isEb;
445  if(detector==1)isEb=true;
446  else isEb=false;
447 
448 
449  if(isEb){
450  for(unsigned int i=0; i<bcCells.size(); ++i){
451  for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
452  if(eb->id().rawId()==bcCells[i].first.rawId()){
453 
454  double energyHit = eb->energy()*bcCells[i].second;
455  DetId id=bcCells[i].first;
456  float eta=geomBar_->getGeometry(id)->getPosition().eta();
457  float phi=geomBar_->getGeometry(id)->getPosition().phi();
458  float dEta = eta - PFClust[max_c].eta();
459  float dPhi = phi - PFClust[max_c].phi();
460  if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; }
461  if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; }
462  numeratorEtaWidth=dEta*dEta*energyHit+numeratorEtaWidth;
463  numeratorPhiWidth=dPhi*dPhi*energyHit+numeratorPhiWidth;
464  denominator=energyHit+denominator;
465  break;
466  }
467  }
468 
469  }
470  }
471  else{
472  for(unsigned int i=0; i<bcCells.size(); ++i){
473  for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
474  if(ee->id().rawId()==bcCells[i].first.rawId()){
475  double energyHit = ee->energy()*bcCells[i].second;
476  DetId id=bcCells[i].first;
477  float eta=geomEnd_->getGeometry(id)->getPosition().eta();
478  float phi=geomEnd_->getGeometry(id)->getPosition().phi();
479  float dEta = eta - PFClust[max_c].eta();
480  float dPhi = phi - PFClust[max_c].phi();
481  if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; }
482  if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; }
483  numeratorEtaWidth=dEta*dEta*energyHit+numeratorEtaWidth;
484  numeratorPhiWidth=dPhi*dPhi*energyHit+numeratorPhiWidth;
485  denominator=energyHit+denominator;
486  break;
487  }
488  }
489  }
490  }
491  float etaWidth=sqrt(numeratorEtaWidth/denominator);
492  float phiWidth=sqrt(numeratorPhiWidth/denominator);
493  widths=make_pair(etaWidth, phiWidth);
494  return widths;
495 }
496 
497 double ggPFClusters::LocalEnergyCorrection(const GBRForest *ReaderLCEB, const GBRForest *ReaderLCEE,reco::CaloCluster PFClust,float beamspotZ ){
498 
499  //double Coor=1;
500  double PFClustCorr=PFClust.energy();
501  float ClusEta=PFClust.eta();
502  float ClusE=PFClust.energy();
503  float ClusPhi=PFClust.phi();
504 
505  std::vector<float>inputs;
506  std::vector< std::pair<DetId, float> >bcCells=PFClust.hitsAndFractions();
507  bool isEb;
508  DetId seedXtalId = bcCells[0].first ;
509  int detector = seedXtalId.subdetId();
510  if(detector==1)isEb=true;
511  else isEb=false;
512  //Shower Shape Variables:
513  Fill5x5Map(bcCells, isEb);
514  float Eseed=get5x5Element(0, 0, bcCells, isEb);
515 
516  float Etop=get5x5Element(0, 1, bcCells, isEb);
517  float Ebottom=get5x5Element(0, -1, bcCells, isEb);
518  float Eleft=get5x5Element(-1, 0, bcCells, isEb);
519  float Eright=get5x5Element(1, 0, bcCells, isEb);
520  float E5x5=0; float E3x3=0;
521  for(int i=-2; i<3; ++i)
522  for(int j=-2; j<3; ++j){
523  float e=get5x5Element(i, j, bcCells, isEb);
524  if(abs(i)<2)E3x3=E3x3+e;
525  E5x5=e+E5x5;
526  }
527  //Fill5x5Map(bcCells, isEb);
528  float E3x1=e5x5_[2][2]+e5x5_[1][2]+e5x5_[3][2];
529  float E1x3=e5x5_[2][2]+e5x5_[2][1]+e5x5_[2][3];
530  float E1x5=e5x5_[2][2]+e5x5_[2][0]+e5x5_[2][1]+e5x5_[2][3]+e5x5_[2][4];
531 
532  float E2x5Top=e5x5_[0][4]+e5x5_[1][4]+e5x5_[2][4]+e5x5_[3][4]+e5x5_[4][4]
533  +e5x5_[0][3]+e5x5_[1][3]+e5x5_[2][3]+e5x5_[3][3]+e5x5_[4][3];
534  //add up bottom edge of 5x5 2 rows
535  float E2x5Bottom=e5x5_[0][0]+e5x5_[1][0]+e5x5_[2][0]+e5x5_[3][0]+e5x5_[4][0]
536  +e5x5_[0][1]+e5x5_[1][1]+e5x5_[2][1]+e5x5_[3][1]+e5x5_[4][1];
537  //add up left edge of 5x5 2 rows
538  float E2x5Left=e5x5_[0][1]+e5x5_[0][1]+e5x5_[0][2]+e5x5_[0][3]+e5x5_[0][4]
539  +e5x5_[1][0]+e5x5_[1][1]+e5x5_[1][2]+e5x5_[1][3]+e5x5_[1][4];
540  //add up right edge of 5x5 2 rows
541  float E2x5Right=e5x5_[4][0]+e5x5_[4][1]+e5x5_[4][2]+e5x5_[4][3]+e5x5_[4][4]
542  +e5x5_[3][0]+e5x5_[3][1]+e5x5_[3][2]+e5x5_[3][3]+e5x5_[3][4];
543  //find max 2x5 from the center
544  float centerstrip=e5x5_[2][2]+e5x5_[2][0]+e5x5_[2][1]+e5x5_[2][3]+e5x5_[2][4];
545  float rightstrip=e5x5_[3][2]+e5x5_[3][0]+e5x5_[3][1]+e5x5_[3][3]+e5x5_[3][4];
546  float leftstrip=e5x5_[1][2]+e5x5_[1][0]+e5x5_[1][1]+e5x5_[1][3]+e5x5_[1][4];
547  float E2x5Max=0;
548  if(rightstrip>leftstrip)E2x5Max=rightstrip+centerstrip;
549  else E2x5Max=leftstrip+centerstrip;
550  //get Local Coordinates
551  if(isEb){
552  float etacry; float phicry; int ieta; int iphi; float thetatilt; float phitilt;
553  int iEtaCrack=3;
554  if(abs(ieta)==1 || abs(ieta)==2 )
555  iEtaCrack=abs(ieta);
556  if(abs(ieta)>2 && abs(ieta)<24)
557  iEtaCrack=3;
558  if(abs(ieta)==24)
559  iEtaCrack=4;
560  if(abs(ieta)==25)
561  iEtaCrack=5;
562  if(abs(ieta)==26)
563  iEtaCrack=6;
564  if(abs(ieta)==27)
565  iEtaCrack=7;
566  if(abs(ieta)>27 && abs(ieta)<44)
567  iEtaCrack=8;
568  if(abs(ieta)==44)
569  iEtaCrack=9;
570  if(abs(ieta)==45)
571  iEtaCrack=10;
572  if(abs(ieta)==46)
573  iEtaCrack=11;
574  if(abs(ieta)==47)
575  iEtaCrack=12;
576  if(abs(ieta)>47 && abs(ieta)<64)
577  iEtaCrack=13;
578  if(abs(ieta)==64)
579  iEtaCrack=14;
580  if(abs(ieta)==65)
581  iEtaCrack=15;
582  if(abs(ieta)==66)
583  iEtaCrack=16;
584  if(abs(ieta)==67)
585  iEtaCrack=17;
586  if(abs(ieta)>67 && abs(ieta)<84)
587  iEtaCrack=18;
588  if(abs(ieta)==84)
589  iEtaCrack=19;
590  if(abs(ieta)==85)
591  iEtaCrack=20;
592  localCoordsEB(PFClust, etacry, phicry, ieta, iphi, thetatilt, phitilt);
593  inputs.push_back(beamspotZ);
594  inputs.push_back(ClusEta/fabs(ClusEta));
595  inputs.push_back(fabs(ClusEta));
596  inputs.push_back(fabs(ClusPhi));
597  inputs.push_back(log(ClusE));
598  inputs.push_back((Eseed/ClusE));
599  inputs.push_back((Etop/ClusE));
600  inputs.push_back((Ebottom/ClusE));
601  inputs.push_back((Eleft/ClusE));
602  inputs.push_back((Eright/ClusE));
603  inputs.push_back(E3x3/ClusE);
604  inputs.push_back(E1x3/ClusE);
605  inputs.push_back(E3x1/ClusE);
606  inputs.push_back(E5x5/ClusE);
607  inputs.push_back(E1x5/ClusE);
608  inputs.push_back(E2x5Max/ClusE);
609  inputs.push_back(E2x5Top/ClusE);
610  inputs.push_back(E2x5Bottom/ClusE);
611  inputs.push_back(E2x5Left/ClusE);
612  inputs.push_back(E2x5Right/ClusE);
613  inputs.push_back(etacry);
614  inputs.push_back(phicry);
615  inputs.push_back(iphi%2);
616  inputs.push_back(ieta%5);
617  inputs.push_back(iphi%20);
618  inputs.push_back(iEtaCrack);
619  int size=inputs.size();
620  float PFInputs[26];
621  for(int i=0; i<size; ++i)PFInputs[i]=inputs[i];
622  PFClustCorr= ReaderLCEB->GetResponse(PFInputs)*ClusE;
623  }
624  else{
625  float xcry; float ycry; int ix; int iy; float thetatilt; float phitilt;
626  localCoordsEE(PFClust, xcry, ycry, ix, iy, thetatilt, phitilt);
627  inputs.push_back(beamspotZ);
628  inputs.push_back(ClusEta/fabs(ClusEta));
629  inputs.push_back(fabs(ClusEta));
630  inputs.push_back(fabs(ClusPhi));
631  inputs.push_back(log(ClusE));
632  inputs.push_back((Eseed/ClusE));
633  inputs.push_back((Etop/ClusE));
634  inputs.push_back((Ebottom/ClusE));
635  inputs.push_back((Eleft/ClusE));
636  inputs.push_back((Eright/ClusE));
637  inputs.push_back(E3x3/ClusE);
638  inputs.push_back(E1x3/ClusE);
639  inputs.push_back(E3x1/ClusE);
640  inputs.push_back(E5x5/ClusE);
641  inputs.push_back(E1x5/ClusE);
642  inputs.push_back(E2x5Max/ClusE);
643  inputs.push_back(E2x5Top/ClusE);
644  inputs.push_back(E2x5Bottom/ClusE);
645  inputs.push_back(E2x5Left/ClusE);
646  inputs.push_back(E2x5Right/ClusE);
647  inputs.push_back(xcry);
648  inputs.push_back(ycry);
649  int size=inputs.size();
650  float PFInputs[22];
651  for(int i=0; i<size; ++i)PFInputs[i]=inputs[i];
652  PFClustCorr= ReaderLCEE->GetResponse(PFInputs) *ClusE;
653  }
654 
655  return PFClustCorr;
656 }
657 void ggPFClusters::BasicClusterPFCandLink(
658  reco::SuperCluster sc,
659  std::vector<reco::PFCandidatePtr>&insideBox,
660  std::vector<DetId>& MatchedRH
661  ){
662  std::vector<reco::PFCandidatePtr>Linked;
663  for(unsigned int p=0; p<insideBox.size(); ++p){
664  math::XYZPointF position_ = insideBox[p]->positionAtECALEntrance();
665  //math::XYZPointF positionvtx(position_.x()+insideBox[p]->vx(),
666  // position_.y()+insideBox[p]->vy(),
667  // position_.z()+insideBox[p]->vz());
668  //position_=positionvtx;
669  //math::XYZVector position_=insideBox[p]->momentum();
670  if(insideBox[p]->pdgId()==22){
671  double Theta = -position_.theta()+0.5*TMath::Pi();
672  double Eta = position_.eta();
673  double Phi = TVector2::Phi_mpi_pi(position_.phi());
674  double X = position_.x();
675  double Y = position_.y();
676  reco::CaloCluster_iterator cit=sc.clustersBegin();
677  std::vector< std::pair<DetId, float> > crystals_vector = (*cit)->hitsAndFractions();
678  DetId seedXtalId = crystals_vector[0].first;
679  int detector = seedXtalId.subdetId();
680  bool isEb=false;
681  float X0 = 0.89; float T0 = 7.4;
682  double depth = X0 * (T0 + log((*cit)->energy()));
683  if(detector==1){
684  X0 = 0.89; T0 = 7.4;
685  depth = X0 * (T0 + log((*cit)->energy()));
686  isEb=true;
687  }
688  else{
689  X0 = 0.89; T0 = 1.2;
690  if(fabs(Eta)<1.653)T0=3.1;
691  depth = X0 * (T0 + log((*cit)->energy()));
692 
693  isEb=false;
694  }
695  crystals_vector.clear();
696  for(; cit!=sc.clustersEnd(); ++cit){
697  bool matchBC=false;
698  crystals_vector = (*cit)->hitsAndFractions();
699 
700 
701  for (unsigned int icry=0; icry<crystals_vector.size(); ++icry){
702 
703  if(isEb){
704 
705  EBDetId crystal(crystals_vector[icry].first);
706  const CaloCellGeometry* cell=geomBar_->getGeometry(crystal);
707 
708  GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
709  //GlobalPoint vtx(insideBox[p]->vx(), insideBox[p]->vy(), insideBox[p]->vz());
710  //GlobalPoint center_pos(oldcenter_pos.x()-vtx.x(),
711  // oldcenter_pos.y()-vtx.y(),
712  // oldcenter_pos.z()-vtx.z());
713  //double EtaCentr = center_pos.eta();
714  double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
715  double PhiWidth = (TMath::Pi()/180.);
716  double ThetaCentr = -center_pos.theta()+0.5*TMath::Pi();
717  double ThetaWidth = (TMath::Pi()/180.)*TMath::Cos(ThetaCentr);
718  double phicry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth;
719  double etacry=(Theta-ThetaCentr)/ThetaWidth;
720  if(fabs(etacry)<0.6 && fabs(phicry)<0.6){
721  Linked.push_back(insideBox[p]);
722  matchBC=true;
723  MatchedRH.push_back(crystals_vector[icry].first);
724  break;
725 
726  }
727 
728  }//Barrel
729  else{
730 
731  EEDetId crystal(crystals_vector[icry].first);
732  const CaloCellGeometry* cell=geomEnd_->getGeometry(crystal);
733 
734  GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
735  //GlobalPoint vtx(insideBox[p]->vx(), insideBox[p]->vy(), insideBox[p]->vz());
736  //GlobalPoint center_pos(oldcenter_pos.x()-vtx.x(), oldcenter_pos.y()-vtx.y(), oldcenter_pos.z()-vtx.z());
737  double XCentr = center_pos.x();
738  double XWidth = 2.59;
739  double xcry = (X-XCentr)/XWidth;
740  double YCentr = center_pos.y();
741  double YWidth = 2.59;
742  double ycry = (Y-YCentr)/YWidth;
743  if(fabs(xcry)<0.6 && fabs(ycry)<0.6){
744  Linked.push_back(insideBox[p]);
745  MatchedRH.push_back(crystals_vector[icry].first);
746  matchBC=true;
747  break;
748  }
749 
750  }//Endcap
751 
752  }
753  if(matchBC)break;
754  }
755  }
756  if(abs(insideBox[p]->pdgId())==211){
757  float drmin = 999.;
758  float deta=999;
759  float dphi=999;
760  reco::CaloCluster_iterator cit=sc.clustersBegin();
761  DetId seedCrystal=(*cit)->hitsAndFractions()[0].first;
762  for(; cit!=sc.clustersEnd(); ++cit){
763  DetId crystals_vector_seed=(*cit)->hitsAndFractions()[0].first;
764  math::XYZVector photon_directionWrtVtx((*cit)->position().x()-insideBox[p]->vx(),
765  (*cit)->position().y()-insideBox[p]->vy(),
766  (*cit)->position().z()-insideBox[p]->vz()
767  );
768  float dR=deltaR(photon_directionWrtVtx.eta(), photon_directionWrtVtx.phi(), position_.eta(), position_.phi());
769  //float dR=deltaR((*cit)->eta(), (*cit)->phi(), position_.eta(), position_.phi());
770  if(dR<drmin){
771  drmin=dR;
772  deta=fabs(photon_directionWrtVtx.eta()-position_.eta());
773  dphi=acos(cos(photon_directionWrtVtx.phi()- position_.phi()));
774  seedCrystal=crystals_vector_seed;
775  }
776  }
777  if(deta<0.05 && dphi<0.07){
778  Linked.push_back(insideBox[p]);
779  MatchedRH.push_back(seedCrystal);
780  }
781  }
782 
783  }
784  insideBox.clear();
785  insideBox=Linked;
786 }
GBRForest::GetResponse
double GetResponse(const float *vector) const
Definition: GBRForest.h:51
TwoPi
const double TwoPi
Definition: CosmicMuonParameters.h:19
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
reco::CaloCluster::position
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:123
i
int i
Definition: DBlmapReader.cc:9
EEDetId::ix
int ix() const
Definition: EEDetId.h:71
reco::CaloCluster::addHitAndFraction
void addHitAndFraction(DetId id, float fraction)
Definition: CaloCluster.h:182
RecoTauValidation_cfi.denominator
tuple denominator
Definition: RecoTauValidation_cfi.py:74
EEDetId::distanceX
static int distanceX(const EEDetId &a, const EEDetId &b)
Definition: EEDetId.cc:602
ggPFClusters::ClusterWidth
virtual std::pair< float, float > ClusterWidth(vector< reco::CaloCluster > &PFClust)
Definition: ggPFClusters.cc:425
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:68
SiPixelRawToDigiRegional_cfi.inputs
tuple inputs
Definition: SiPixelRawToDigiRegional_cfi.py:8
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1108
edm::SortedCollection< EcalRecHit >::const_iterator
std::vector< EcalRecHit >::const_iterator const_iterator
Definition: SortedCollection.h:80
ggPFClusters::localCoordsEE
virtual void localCoordsEE(reco::CaloCluster clus, float &xcry, float &ycry, int &ix, int &iy, float &thetatilt, float &phitilt)
Definition: ggPFClusters.cc:256
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
ggPFClusters::geomBar_
const CaloSubdetectorGeometry * geomBar_
Definition: ggPFClusters.h:69
X
#define X(str)
Definition: MuonsGrabber.cc:49
abs
#define abs(x)
Definition: mlp_lapack.h:159
reco::Photon::superCluster
reco::SuperClusterRef superCluster() const
Ref to SuperCluster.
Definition: Photon.cc:59
reco::CaloCluster::hitsAndFractions
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:189
Phi_mpi_pi
double Phi_mpi_pi(double x)
Definition: JetUtil.h:24
math::XYZPointF
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:11
TruncatedPyramid::getPhiAxis
CCGFloat getPhiAxis() const
Definition: TruncatedPyramid.cc:77
eta
T eta() const
Definition: Basic3DVectorLD.h:172
CaloCellGeometry::getPosition
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.cc:46
ggPFClusters::BasicClusterPFCandLink
void BasicClusterPFCandLink(reco::SuperCluster sc, std::vector< reco::PFCandidatePtr > &insideBox, std::vector< DetId > &MatchedRH)
Definition: ggPFClusters.cc:657
ggPFClusters::SumPFRecHits
virtual float SumPFRecHits(std::vector< std::pair< DetId, float > > &bcCells, bool isEB)
Definition: ggPFClusters.cc:57
ggPFClusters::e5x5_
float e5x5_[5][5]
Definition: ggPFClusters.h:71
reco::CaloCluster::eta
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:160
ggPFClusters::get5x5Element
virtual float get5x5Element(int i, int j, std::vector< std::pair< DetId, float > > &bcCells, bool isEB)
Definition: ggPFClusters.cc:317
ggPFClusters::~ggPFClusters
virtual ~ggPFClusters()
Definition: ggPFClusters.cc:30
CaloSubdetectorGeometry::getGeometry
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
Definition: CaloSubdetectorGeometry.cc:43
reco::deltaR
double deltaR(double eta1, double phi1, double eta2, double phi2)
Definition: deltaR.h:19
muon::overlap
bool overlap(const reco::Muon &muon1, const reco::Muon &muon2, double pullX=1.0, double pullY=1.0, bool checkAdjacentChambers=false)
Definition: MuonSelectors.cc:640
ggPFClusters::EBReducedRecHits_
Handle< EcalRecHitCollection > EBReducedRecHits_
Definition: ggPFClusters.h:67
PV3DBase::theta
Geom::Theta< T > theta() const
Definition: PV3DBase.h:74
ggPFClusters::getPFSuperclusterOverlap
float getPFSuperclusterOverlap(reco::CaloCluster PFClust, reco::SuperCluster sc)
Definition: ggPFClusters.cc:84
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:46
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
ggPFClusters::ggPFClusters
ggPFClusters(edm::Handle< EcalRecHitCollection > &EBReducedRecHits, edm::Handle< EcalRecHitCollection > &EEReducedRecHits, const CaloSubdetectorGeometry *geomBar, const CaloSubdetectorGeometry *geomEnd)
Definition: ggPFClusters.cc:13
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:234
ggPFClusters::Fill5x5Map
virtual void Fill5x5Map(std::vector< std::pair< DetId, float > > &bcCells, bool isEB)
Definition: ggPFClusters.cc:330
EBDetId::distanceEta
static int distanceEta(const EBDetId &a, const EBDetId &b)
Definition: EBDetId.cc:185
ggPFClusters::localCoordsEB
virtual void localCoordsEB(reco::CaloCluster clus, float &etacry, float &phicry, int &ieta, int &iphi, float &thetatilt, float &phitilt)
Definition: ggPFClusters.cc:195
dPhi
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
EBDetId::distancePhi
static int distancePhi(const EBDetId &a, const EBDetId &b)
Definition: EBDetId.cc:193
ggPFClusters::getPFClusters
virtual vector< reco::CaloCluster > getPFClusters(reco::SuperCluster)
Definition: ggPFClusters.cc:35
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
ggPFClusters::LocalEnergyCorrection
double LocalEnergyCorrection(const GBRForest *ReaderLCEB, const GBRForest *ReaderLCEE, reco::CaloCluster PFClust, float beamspotZ)
Definition: ggPFClusters.cc:497
j
int j
Definition: DBlmapReader.cc:9
reco::CaloCluster::energy
double energy() const
cluster energy
Definition: CaloCluster.h:120
EEDetId::iy
int iy() const
Definition: EEDetId.h:77
ggPFClusters::FindSeed
virtual DetId FindSeed(std::vector< std::pair< DetId, float > > &bcCells, bool isEB)
Definition: ggPFClusters.cc:388
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:44
EEDetId::distanceY
static int distanceY(const EEDetId &a, const EEDetId &b)
Definition: EEDetId.cc:607
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
TruncatedPyramid::getThetaAxis
CCGFloat getThetaAxis() const
Definition: TruncatedPyramid.cc:71
first
bool first
Definition: L1TdeRCT.cc:94
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
DetId
Definition: DetId.h:20
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
ggPFClusters::geomEnd_
const CaloSubdetectorGeometry * geomEnd_
Definition: ggPFClusters.h:70
TruncatedPyramid::getPosition
const GlobalPoint getPosition(CCGFloat depth) const
Definition: TruncatedPyramid.cc:65
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
TruncatedPyramid
A base class to handle the particular shape of Ecal Xtals. Taken from ORCA Calorimetry Code...
Definition: TruncatedPyramid.h:19
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:75
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
reco::SuperCluster::clustersBegin
CaloCluster_iterator clustersBegin() const
fist iterator over BasicCluster constituents
Definition: SuperCluster.h:65
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
benchmark_cfg.pdgId
tuple pdgId
Definition: benchmark_cfg.py:18
reco::CaloCluster::phi
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:163
X0
static const double X0
Definition: CastorGeometryData.h:31
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
ggPFClusters::EEReducedRecHits_
Handle< EcalRecHitCollection > EEReducedRecHits_
Definition: ggPFClusters.h:68
ggPFClusters::PFRecHitsSCOverlap
virtual float PFRecHitsSCOverlap(std::vector< std::pair< DetId, float > > &bcCells1, std::vector< std::pair< DetId, float > > &bcCells2, bool isEB)
Definition: ggPFClusters.cc:137
reco::SuperCluster::clustersEnd
CaloCluster_iterator clustersEnd() const
last iterator over BasicCluster constituents
Definition: SuperCluster.h:68
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5