CMS 3D CMS Logo

DetIdAssociator.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HTrackAssociator
4 // Class: HDetIdAssociator
5 //
6 /*
7 
8  Description: <one line class summary>
9 
10  Implementation:
11  <Notes on implementation>
12 */
13 //
14 // Original Author: Dmytro Kovalskyi
15 // Modified for ECAL+HCAL by: Michal Szleper
16 // Created: Fri Apr 21 10:59:41 PDT 2006
17 //
18 //
19 
20 
22 
23 #include <memory>
24 
25 
26 // surfaces is a vector of GlobalPoint representing outermost point on a cylinder
27 std::vector<GlobalPoint> HDetIdAssociator::getTrajectory( const FreeTrajectoryState& ftsStart,
28  const std::vector<GlobalPoint>& surfaces)
29 {
30  check_setup();
31  std::vector<GlobalPoint> trajectory;
32  TrajectoryStateOnSurface tSOSDest;
33  FreeTrajectoryState ftsCurrent = ftsStart;
34 
35  for(std::vector<GlobalPoint>::const_iterator surface_iter = surfaces.begin();
36  surface_iter != surfaces.end(); surface_iter++) {
37  // this stuff is some weird pointer, which destroy itself
38  std::unique_ptr<Cylinder> cylinder = std::make_unique<Cylinder>(surface_iter->perp(),
39  Surface::PositionType(0,0,0),
41  std::unique_ptr<Plane> forwardEndcap = std::make_unique<Plane>(Surface::PositionType(0,0,surface_iter->z()),
43  std::unique_ptr<Plane> backwardEndcap = std::make_unique<Plane>(Surface::PositionType(0,0,-surface_iter->z()),
45 
46 
47  LogTrace("StartingPoint")<< "Propagate from "<< "\n"
48  << "\tx: " << ftsStart.position().x()<< "\n"
49  << "\ty: " << ftsStart.position().y()<< "\n"
50  << "\tz: " << ftsStart.position().z()<< "\n"
51  << "\tmomentum eta: " << ftsStart.momentum().eta()<< "\n"
52  << "\tmomentum phi: " << ftsStart.momentum().phi()<< "\n"
53  << "\tmomentum: " << ftsStart.momentum().mag()<< "\n";
54 
55  float tanTheta = ftsCurrent.momentum().perp()/ftsCurrent.momentum().z();
56  float corner = surface_iter->perp()/surface_iter->z();
57 /*
58  std::cout<<"Propagate from "<< "\n"
59  << "\tx: " << ftsCurrent.position().x()<< "\n"
60  << "\ty: " << ftsCurrent.position().y()<< "\n"
61  << "\tz: " << ftsCurrent.position().z()<< "\n"
62  << "\tz: " << ftsCurrent.position().perp()<< "\n"
63  << "\tz: " << tanTheta<<" "<< corner <<"\n"
64  << "\tmomentum eta: " << ftsCurrent.momentum().eta()<< "\n"
65  << "\tmomentum phi: " << ftsCurrent.momentum().phi()<< "\n"
66  << "\tmomentum: " << ftsCurrent.momentum().mag()<<std::endl;
67 */
68  // First propage the track to the cylinder if |eta|<1, othewise to the encap
69  // and correct depending on the result
70  int ibar = 0;
71  if (fabs(tanTheta) > corner)
72  {
73  tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
74  // std::cout<<" Propagate to cylinder "<<std::endl;
75  }
76  else if(tanTheta > 0.)
77  {tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap); ibar=1; }
78  else
79  {tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap); ibar=-1; }
80 
81 // std::cout<<" Trajectory valid? "<<tSOSDest.isValid()<<" First propagation in "<<ibar<<std::endl;
82 
83  if(! tSOSDest.isValid() )
84  {
85 // barrel
86  if(ibar == 0){
87  if (tanTheta < 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*forwardEndcap);
88  if (tanTheta >= 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*backwardEndcap);
89  }
90  else
91  {
92  tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
93  }
94  } else
95  {
96 // missed target
97  if(abs(ibar) > 0)
98  {
99  if(tSOSDest.globalPosition().perp() > surface_iter->perp())
100  {
101  tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
102  }
103  }
104  else
105  {
106  if (tanTheta < 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*forwardEndcap);
107  if (tanTheta >= 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*backwardEndcap);
108  }
109  }
110 
111 
112  // If missed the target, propagate to again
113 // if ((!tSOSDest.isValid()) && point.perp() > surface_iter->perp())
114 // {tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);std::cout<<" Propagate again 1 "<<std::endl;}
115 // std::cout<<" Track is ok after repropagation to cylinder or not? "<<tSOSDest.isValid()<<std::endl;
116 // if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()>0. && fabs(ftsStart.momentum().eta())>1.)
117 // {tSOSDest = ivProp_->propagate(ftsStart, *forwardEndcap);std::cout<<" Propagate again 2 "<<std::endl;}
118 // std::cout<<" Track is ok after repropagation forward or not? "<<tSOSDest.isValid()<<std::endl;
119 // if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()<0.&&fabs(ftsStart.momentum().eta())>1.)
120 // {tSOSDest = ivProp_->propagate(ftsStart, *backwardEndcap);std::cout<<" Propagate again 3 "<<std::endl;}
121 // std::cout<<" Track is after repropagation backward ok or not? "<<tSOSDest.isValid()<<std::endl;
122 
123 
124  if (! tSOSDest.isValid()) return trajectory;
125 
126 // std::cout<<" Propagate reach something"<<std::endl;
127  LogTrace("SuccessfullPropagation") << "Great, I reached something." << "\n"
128  << "\tx: " << tSOSDest.freeState()->position().x() << "\n"
129  << "\ty: " << tSOSDest.freeState()->position().y() << "\n"
130  << "\tz: " << tSOSDest.freeState()->position().z() << "\n"
131  << "\teta: " << tSOSDest.freeState()->position().eta() << "\n"
132  << "\tphi: " << tSOSDest.freeState()->position().phi() << "\n";
133 
134 // std::cout<<" The position of trajectory "<<tSOSDest.freeState()->position().perp()<<" "<<tSOSDest.freeState()->position().z()<<std::endl;
135 
136  GlobalPoint point = tSOSDest.freeState()->position();
137  point = tSOSDest.freeState()->position();
138  ftsCurrent = *tSOSDest.freeState();
139  trajectory.push_back(point);
140  }
141  return trajectory;
142 }
143 
144 //------------------------------------------------------------------------------
145 std::set<DetId> HDetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction,
146  const int idR)
147 {
148  std::set<DetId> set;
149  check_setup();
150  int nDets=0;
151  if (! theMap_) buildMap();
152  LogTrace("MatchPoint") << "point (eta,phi): " << direction.eta() << "," << direction.phi() << "\n";
153  int ieta = iEta(direction);
154  int iphi = iPhi(direction);
155 
156  LogTrace("MatchPoint") << "(ieta,iphi): " << ieta << "," << iphi << "\n";
157 
158  if (ieta>=0 && ieta<nEta_ && iphi>=0 && iphi<nPhi_){
159  set = (*theMap_)[ieta][iphi];
160  nDets++;
161  if (idR>0){
162  LogTrace("MatchPoint") << "Add neighbors (ieta,iphi): " << ieta << "," << iphi << "\n";
163  //add neighbors
164  int maxIEta = ieta+idR;
165  int minIEta = ieta-idR;
166  if(maxIEta>=nEta_) maxIEta = nEta_-1;
167  if(minIEta<0) minIEta = 0;
168  int maxIPhi = iphi+idR;
169  int minIPhi = iphi-idR;
170  if(minIPhi<0) {
171  minIPhi+=nPhi_;
172  maxIPhi+=nPhi_;
173  }
174  LogTrace("MatchPoint") << "\tieta (min,max): " << minIEta << "," << maxIEta<< "\n";
175  LogTrace("MatchPoint") << "\tiphi (min,max): " << minIPhi << "," << maxIPhi<< "\n";
176  for (int i=minIEta;i<=maxIEta;i++)
177  for (int j=minIPhi;j<=maxIPhi;j++) {
178  if( i==ieta && j==iphi) continue; // already in the set
179  set.insert((*theMap_)[i][j%nPhi_].begin(),(*theMap_)[i][j%nPhi_].end());
180  nDets++;
181  }
182  }
183  }
184 // if(set.size() > 0) {
185 // if (ieta+idR<55 && ieta-idR>14 && set.size() != (2*idR+1)*(2*idR+1)){
186 // std::cout<<" RRRA: "<<set.size()<<" DetIds in region "<<ieta<<" "<<iphi<<std::endl;
187 // for( std::set<DetId>::const_iterator itr=set.begin(); itr!=set.end(); itr++) {
188 // GlobalPoint point = getPosition(*itr);
189 // std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<" "<<iEta(point)<<" "<<iPhi(point)<<std::endl;
190 // }
191 // }
192 // else {
193 // std::cout <<" HDetIdAssociator::getDetIdsCloseToAPoint::There are strange days "<<std::endl;
194 // }
195  return set;
196 }
197 
198 //------------------------------------------------------------------------------
200 {
201 // unequal bin sizes for endcap, following HCAL geometry
202  int iEta1 = int(point.eta()/etaBinSize_ + nEta_/2);
203  if (point.eta()>1.827 && point.eta()<=1.830) return iEta1-1;
204  else if (point.eta()>1.914 && point.eta()<=1.930) return iEta1-1;
205  else if (point.eta()>2.001 && point.eta()<=2.043) return iEta1-1;
206  else if (point.eta()>2.088 && point.eta()<=2.172) return iEta1-1;
207  else if (point.eta()>2.175 && point.eta()<=2.262) return iEta1-1;
208  else if (point.eta()>2.262 && point.eta()<=2.332) return iEta1-2;
209  else if (point.eta()>2.332 && point.eta()<=2.349) return iEta1-1;
210  else if (point.eta()>2.349 && point.eta()<=2.436) return iEta1-2;
211  else if (point.eta()>2.436 && point.eta()<=2.500) return iEta1-3;
212  else if (point.eta()>2.500 && point.eta()<=2.523) return iEta1-2;
213  else if (point.eta()>2.523 && point.eta()<=2.610) return iEta1-3;
214  else if (point.eta()>2.610 && point.eta()<=2.650) return iEta1-4;
215  else if (point.eta()>2.650 && point.eta()<=2.697) return iEta1-3;
216  else if (point.eta()>2.697 && point.eta()<=2.784) return iEta1-4;
217  else if (point.eta()>2.784 && point.eta()<=2.868) return iEta1-5;
218  else if (point.eta()>2.868 && point.eta()<=2.871) return iEta1-4;
219  else if (point.eta()>2.871 && point.eta()<=2.958) return iEta1-5;
220  else if (point.eta()>2.958) return iEta1-6;
221  else if (point.eta()<-1.827 && point.eta()>=-1.830) return iEta1+1;
222  else if (point.eta()<-1.914 && point.eta()>=-1.930) return iEta1+1;
223  else if (point.eta()<-2.001 && point.eta()>=-2.043) return iEta1+1;
224  else if (point.eta()<-2.088 && point.eta()>=-2.172) return iEta1+1;
225  else if (point.eta()<-2.175 && point.eta()>=-2.262) return iEta1+1;
226  else if (point.eta()<-2.262 && point.eta()>=-2.332) return iEta1+2;
227  else if (point.eta()<-2.332 && point.eta()>=-2.349) return iEta1+1;
228  else if (point.eta()<-2.349 && point.eta()>=-2.436) return iEta1+2;
229  else if (point.eta()<-2.436 && point.eta()>=-2.500) return iEta1+3;
230  else if (point.eta()<-2.500 && point.eta()>=-2.523) return iEta1+2;
231  else if (point.eta()<-2.523 && point.eta()>=-2.610) return iEta1+3;
232  else if (point.eta()<-2.610 && point.eta()>=-2.650) return iEta1+4;
233  else if (point.eta()<-2.650 && point.eta()>=-2.697) return iEta1+3;
234  else if (point.eta()<-2.697 && point.eta()>=-2.784) return iEta1+4;
235  else if (point.eta()<-2.784 && point.eta()>=-2.868) return iEta1+5;
236  else if (point.eta()<-2.868 && point.eta()>=-2.871) return iEta1+4;
237  else if (point.eta()<-2.871 && point.eta()>=-2.958) return iEta1+5;
238  else if (point.eta()<-2.349) return iEta1+6;
239  else return iEta1;
240 }
241 
242 //------------------------------------------------------------------------------
244 {
245  double pi=4*atan(1.);
246  int iPhi1 = int((double(point.phi())+pi)/(2*pi)*nPhi_);
247  return iPhi1;
248 }
249 
250 //------------------------------------------------------------------------------
252 {
253 // modified version: take only detector central position
254  check_setup();
255  LogTrace("HDetIdAssociator")<<"building map" << "\n";
256  if(theMap_) delete theMap_;
257  theMap_ = new std::vector<std::vector<std::set<DetId> > >(nEta_,std::vector<std::set<DetId> >(nPhi_));
258  int numberOfDetIdsOutsideEtaRange = 0;
259  int numberOfDetIdsActive = 0;
260  std::set<DetId> validIds = getASetOfValidDetIds();
261  for (std::set<DetId>::const_iterator id_itr = validIds.begin(); id_itr!=validIds.end(); id_itr++) {
262 // std::vector<GlobalPoint> points = getDetIdPoints(*id_itr);
263  GlobalPoint point = getPosition(*id_itr);
264 // reject fake DetIds (eta=0 - what are they anyway???)
265  if(point.eta()==0)continue;
266 
267  int ieta = iEta(point);
268  int iphi = iPhi(point);
269  int etaMax(-1);
270  int etaMin(nEta_);
271  int phiMax(-1);
272  int phiMin(nPhi_);
273  if ( iphi >= nPhi_ ) iphi = iphi % nPhi_;
274  assert (iphi>=0);
275  if ( etaMin > ieta) etaMin = ieta;
276  if ( etaMax < ieta) etaMax = ieta;
277  if ( phiMin > iphi) phiMin = iphi;
278  if ( phiMax < iphi) phiMax = iphi;
279 // for abs(eta)>1.8 one tower covers two phi segments
280  if ((ieta>54||ieta<15) && iphi%2==0) phiMax++;
281  if ((ieta>54||ieta<15) && iphi%2==1) phiMin--;
282 
283  if (etaMax<0||phiMax<0||etaMin>=nEta_||phiMin>=nPhi_) {
284  LogTrace("HDetIdAssociator")<<"Out of range: DetId:" << id_itr->rawId() <<
285  "\n\teta (min,max): " << etaMin << "," << etaMax <<
286  "\n\tphi (min,max): " << phiMin << "," << phiMax <<
287  "\nTower id: " << id_itr->rawId() << "\n";
288  numberOfDetIdsOutsideEtaRange++;
289  continue;
290  }
291 
292  if (phiMax-phiMin > phiMin+nPhi_-phiMax){
293  phiMin += nPhi_;
294  std::swap(phiMin,phiMax);
295  }
296  for(int ieta = etaMin; ieta <= etaMax; ieta++)
297  for(int iphi = phiMin; iphi <= phiMax; iphi++)
298  (*theMap_)[ieta][iphi%nPhi_].insert(*id_itr);
299  numberOfDetIdsActive++;
300  }
301  LogTrace("HDetIdAssociator") << "Number of elements outside the allowed range ( |eta|>"<<
302  nEta_/2*etaBinSize_ << "): " << numberOfDetIdsOutsideEtaRange << "\n";
303  LogTrace("HDetIdAssociator") << "Number of active DetId's mapped: " <<
304  numberOfDetIdsActive << "\n";
305 }
306 
307 //------------------------------------------------------------------------------
308 std::set<DetId> HDetIdAssociator::getDetIdsInACone(const std::set<DetId>& inset,
309  const std::vector<GlobalPoint>& trajectory,
310  const double dR)
311 {
312 // modified version: if dR<0, returns 3x3 towers around the input one (Michal)
313  check_setup();
314  std::set<DetId> outset;
315 
316  if(dR>=0) {
317  for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
318  for(std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end(); point_iter++)
319  if (nearElement(*point_iter,*id_iter,dR)) outset.insert(*id_iter);
320  }
321  else {
322  if (inset.size()!=1) return outset;
323  std::set<DetId>::const_iterator id_inp = inset.begin();
324  int ieta;
325  int iphi;
326  GlobalPoint point = getPosition(*id_inp);
327  ieta = iEta(point);
328  iphi = iPhi(point);
329  for (int i=ieta-1;i<=ieta+1;i++) {
330  for (int j=iphi-1;j<=iphi+1;j++) {
331 // if( i==ieta && j==iphi) continue;
332  if( i<0 || i>=nEta_) continue;
333  int j2fill = j%nPhi_;
334  if(j2fill<0) j2fill+=nPhi_;
335  if((*theMap_)[i][j2fill].empty())continue;
336  outset.insert((*theMap_)[i][j2fill].begin(),(*theMap_)[i][j2fill].end());
337  }
338  }
339  }
340 
341 // if (outset.size() > 0) {
342 // std::cout<<" RRRA: DetIds in cone:"<<std::endl;
343 // for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
344 // GlobalPoint point = getPosition(*itr);
345 // std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
346 // }
347 // }
348 
349  return outset;
350 }
351 
352 //------------------------------------------------------------------------------
353 std::set<DetId> HDetIdAssociator::getCrossedDetIds(const std::set<DetId>& inset,
354  const std::vector<GlobalPoint>& trajectory)
355 {
356  check_setup();
357  std::set<DetId> outset;
358  for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
359  for(std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end(); point_iter++)
360  if (insideElement(*point_iter, *id_iter)) outset.insert(*id_iter);
361  return outset;
362 }
363 
364 //------------------------------------------------------------------------------
365 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
367 {
368 // returns the most energetic tower in the NxN box (Michal)
369  check_setup();
370  std::set<DetId> outset;
371  std::set<DetId>::const_iterator id_max = inset.begin();
372  double Ehadmax=0;
373 
374  for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
375  DetId id(*id_iter);
376 // GlobalPoint point = getPosition(*id_iter);
377 // int ieta = iEta(point);
378 // int iphi = iPhi(point);
379  CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
380  if(tower != (*caloTowers).end() && tower->hadEnergy()>Ehadmax) {
381  id_max = id_iter;
382  Ehadmax = tower->hadEnergy();
383  }
384  }
385 
386  if (Ehadmax > 0) outset.insert(*id_max);
387 
388 // if (outset.size() > 0) {
389 // std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
390 // for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
391 // GlobalPoint point = getPosition(*itr);
392 // std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
393 // }
394 // }
395 
396  return outset;
397 }
398 
399 //------------------------------------------------------------------------------
400 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
402 {
403 // returns the most energetic tower in the NxN box - from RecHits (Michal)
404  check_setup();
405  std::set<DetId> outset;
406  std::set<DetId>::const_iterator id_max = inset.begin();
407  double Ehadmax=0;
408 
409  for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
410  DetId id(*id_iter);
411 // GlobalPoint point = getPosition(*id_iter);
412 // int ieta = iEta(point);
413 // int iphi = iPhi(point);
414  HBHERecHitCollection::const_iterator hit = (*recHits).find(id);
415  if(hit != (*recHits).end() && hit->energy()>Ehadmax) {
416  id_max = id_iter;
417  Ehadmax = hit->energy();
418  }
419  }
420 
421  if (Ehadmax > 0) outset.insert(*id_max);
422 
423 // if (outset.size() > 0) {
424 // std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
425 // for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
426 // GlobalPoint point = getPosition(*itr);
427 // std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
428 // }
429 // }
430 
431  return outset;
432 }
virtual void check_setup()
virtual void buildMap()
virtual std::vector< GlobalPoint > getTrajectory(const FreeTrajectoryState &, const std::vector< GlobalPoint > &)
virtual std::set< DetId > getDetIdsCloseToAPoint(const GlobalPoint &, const int idR=0)
const double etaBinSize_
T perp() const
Definition: PV3DBase.h:72
virtual GlobalPoint getPosition(const DetId &)=0
std::vector< std::vector< std::set< DetId > > > * theMap_
virtual bool insideElement(const GlobalPoint &, const DetId &)=0
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
std::vector< CaloTower >::const_iterator const_iterator
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
const Double_t pi
virtual bool nearElement(const GlobalPoint &point, const DetId &id, const double distance)
T mag() const
Definition: PV3DBase.h:67
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
FreeTrajectoryState const * freeState(bool withErrors=true) const
T z() const
Definition: PV3DBase.h:64
virtual int iPhi(const GlobalPoint &)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual std::set< DetId > getDetIdsInACone(const std::set< DetId > &, const std::vector< GlobalPoint > &trajectory, const double)
virtual std::set< DetId > getCrossedDetIds(const std::set< DetId > &, const std::vector< GlobalPoint > &trajectory)
#define end
Definition: vmac.h:39
virtual std::set< DetId > getASetOfValidDetIds()=0
virtual std::set< DetId > getMaxEDetId(const std::set< DetId > &, edm::Handle< CaloTowerCollection > caloTowers)
GlobalVector momentum() const
#define LogTrace(id)
Propagator * ivProp_
virtual int iEta(const GlobalPoint &)
Definition: DetId.h:18
GlobalPoint position() const
Point3DBase< float, GlobalTag > PositionType
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
T eta() const
Definition: PV3DBase.h:76
#define begin
Definition: vmac.h:32
TkRotation< float > RotationType
T x() const
Definition: PV3DBase.h:62
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5