ATLAS Offline Software
IDScanZFinderInternal.h
Go to the documentation of this file.
1 // emacs: this is -*- c++ -*-
2 
3 /*
4  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
5 */
6 
8 //
9 // filename: IDScanZFinderInternal.h
10 //
11 // author: Nikos Konstantinidis <n.konstantinidis@ucl.ac.uk>
12 //
13 //
14 // Description: NON ATHENA Internals for the ZFinder primary z vertex finding AlgTool
15 //
16 // -------------------------------
17 // ATLAS Collaboration
19 
20 
21 #ifndef IDSCANZFINDER_IDSCANZFINDERINTERNAL_H
22 #define IDSCANZFINDER_IDSCANZFINDERINTERNAL_H
23 
24 #include <cmath>
25 #include <vector>
26 #include <string>
27 
28 #include "ZFinderConstants.h"
29 #include "PhiSlice.h"
31 
32 static const std::string mZFIVER("$Id: IDScanZFinderInternal.h 700483 2015-10-14 09:37:16Z sutt $");
33 
34 
35 
36 
37 namespace Run1 {
38 
39 template<class SpacePoint> class IDScanZFinderInternal
40 {
41 public:
42 
43  struct vertex {
44  vertex( double z, double weight ) : _z(z), _weight(weight) { }
45  double _z;
46  double _weight;
47  };
48 
49 public:
50 
51  IDScanZFinderInternal( const std::string&, const std::string& );
53  // void initializeInternal(long maxLayers=20, long lastBarrel=7);
54  void initializeInternal(long maxLayers, long lastBarrel);
55 
56  std::vector<vertex>* findZInternal( const std::vector<const SpacePoint *>& spVec, const IRoiDescriptor& roi);
57 
58  const std::vector< std::vector<long> >* GetnHisto() { return m_nHisto; }
59  const std::vector< std::vector<double> >* GetzHisto() { return m_zHisto; }
60 
61  long GetNMax() { return m_NMax; }
62 
63  // void setLayers(long maxLayers=20, long lastBarrelLayer=7) {
64  void setLayers(long maxLayers, long lastBarrelLayer) {
65  m_IdScan_MaxNumLayers = maxLayers; // dphiEC depends on this value
66  m_IdScan_LastBrlLayer = lastBarrelLayer; // dphiBrl depends on this value
67  }
68 
69 protected: // member functions
70 
71  // fills phi, rho, z, layer of spacepoints to simple vectors
72 
73  long fillVectors( const std::vector<const SpacePoint *>& spVec, const IRoiDescriptor& roi,
74  std::vector<double>& phi, std::vector<double>& rho, std::vector<double>& zed,
75  std::vector<long>& lyr, std::vector<long>& filledLayers);
76 
77  const std::string& getType() const { return m_Type; }
78  const std::string& getName() const { return m_Name; }
79  const std::string& getVersion() const { return mZFIVER; }
80 
81  int GetInternalStatus() const { return m_Status; }
82  int SetInternalStatus(int s) { m_Status = s; return m_Status; }
83 
84  double computeZV(double r1, double z1, double r2, double z2) const;
85  double computeZV(double r1, double p1, double z1, double r2, double p2, double z2) const;
86 
87  void SetReturnValue(double d) { m_returnval=d; }
88  double GetReturnValue() const { return m_returnval; }
89 
90 
91 protected: // data members
92 
96 
97 
98  // To be read in from jobOptions by IdScanMain
99 
100  double m_invPhiSliceSize; // the inverse size of the phi slices
101  long m_NumPhiSlices; // the number of phi slices, given the width of the RoI
102 
103  double m_phiBinSize; // the size of the phi slices
104  bool m_forcePhiBinSize; // forces the phi bin size to be be used as configured even if below reasonable limit
105  double m_usedphiBinSize; // the size of the phi slices
106  double m_ROIphiWidth; // the phi width of the ROI
107  double m_usedROIphiWidth; // the phi width of the ROI
108  double m_minZBinSize; // z-histo bin size: m_minZBinSize+|etaRoI|*m_zBinSizeEtaCoeff (to account for worse resolution in high eta)
109  double m_zBinSizeEtaCoeff; // z-histo bin size: m_minZBinSize+|etaRoI|*m_zBinSizeEtaCoeff (to account for worse resolution in high eta)
110 
111  long m_numberOfPeaks; // how many z-positions to return in findZ
112 
113  bool m_pixOnly; // use only Pixel space points
114 
115  std::string m_Type; // type information for internal book keeping
116  std::string m_Name; // name information for the same
117 
118  int m_Status; // return status of the algorithm: 0=ok, -1=error
119 
120  bool m_chargeAware ; // maintain separate sets of z histos for + & - tracks
121  bool m_zHistoPerPhi; // maintain one set of z histos per each phi slice
122 
123  double m_dphideta; // how, as a function of eta, the number of phi neighbours decreases
124  double m_neighborMultiplier; // extra factor to manually increase the number of phi neighbors
125  // long m_extraPhi[IdScan_MaxNumLayers][IdScan_MaxNumLayers]; // number of phi neighbours to look at
126  std::vector< std::vector<long> > m_extraPhi; // ( IdScan_MaxNumLayers, std::vector<long>(IdScan_MaxNumLayers) ); // number of phi neighbours to look at
127 
128  // access the ZFinder histogram from outside the findZInternal method
129 
130  std::vector < std::vector<long> > m_nHisto[2]; // the actual z histogram count of pairs
131  std::vector < std::vector<double> > m_zHisto[2]; // summed z position histograms
132 
133  long m_NMax; // maximum number of z histogram entries
134 
135  int m_nFirstLayers; // When the pairs of SPs are made, the inner SP comes from up to this "filled layer"
136 
137  double m_vrtxDistCut; // The minimum fractional distance between two output vertices
138  double m_vrtxMixing ; // If two vertices are found to be too close, "mix" the second into first
139  int m_nvrtxSeparation; // The minimum number of zbins that any two output vertices should by separated by
140  bool m_preferCentralZ; // Among peaks of same height, should precedence go to the one with smaller |z|
141 
142  bool m_trustSPprovider; // Should we re-extract the RoI phi range from the phis of the SPs from the SPP
143 
144  double m_returnval; // return value for algorithm
145 
147 
149  double m_tripletDZ;
150  double m_tripletDK;
151  double m_halfTripletDK; // replaces m_tripletDK internally to avoid unnecessary multiplication by 2 in curvature calculation, without changing the interface
152  double m_tripletDP;
153 
155 
157 
158 };
159 
160 
161 
162 // template<typename T> long IDScanZFinderInternal<T>::IdScan_MaxNumLayers = 20; // dphiEC depends on this value
163 // template<typename T> long IDScanZFinderInternal<T>::IdScan_LastBrlLayer = 7;
164 
165 
169 
170 
171 template<class SpacePoint>
173  const std::string& name)
174  : m_Type (type),
175  m_Name (name)
176 {
177  m_phiBinSize = 0.2 ;
178  m_forcePhiBinSize = false ;
180  m_pixOnly = false ;
181  m_ROIphiWidth = 0.2 ;
183  m_minZBinSize = 0.2 ;
184  m_zBinSizeEtaCoeff = 0.1 ;
185  m_dphideta = -0.02 ;
186  m_neighborMultiplier = 1. ;
187  m_numberOfPeaks = 1 ;
188  m_chargeAware = true ;
189  m_zHistoPerPhi = true ;
190  m_NMax = 0 ;
191  m_nFirstLayers = 3 ;
192  m_vrtxDistCut = 0. ;
193  m_vrtxMixing = 0. ;
194  m_nvrtxSeparation = 0 ;
195  m_preferCentralZ = false ;
196  m_trustSPprovider = true ;
197  m_returnval = 0 ;
198  m_fullScanMode = false ;
199  m_tripletMode = 0 ;
200  m_tripletDZ = 25. ;
201  m_tripletDK = 0.005 ;
203  m_tripletDP = 0.05 ;
204 
205  // m_applyWeightThreshold = false;
206  m_weightThreshold = 0;
207 
208 
210  m_IdScan_MaxNumLayers = 20; // dphiEC depends on this value !!! 19 without IBL, 20 with IBL!!
211  m_IdScan_LastBrlLayer = 7; // dphiBrl depends on this value
212 
213  m_invPhiSliceSize = 0;
214  m_NumPhiSlices = 0;
215 
219 
220  // cppcheck-suppress missingReturn; false positive
221  m_Status = 0;
222 }
223 
224 
225 
226 template<class SpacePoint>
227 void IDScanZFinderInternal<SpacePoint>::initializeInternal(long maxLayers, long lastBarrel )
228 {
229  m_halfTripletDK = 0.5*m_tripletDK;
230 
231  m_IdScan_MaxNumLayers = maxLayers;
232  m_IdScan_LastBrlLayer = lastBarrel;
233 
234  // std::cout << "m_IdScan_MaxNumLayers " << m_IdScan_MaxNumLayers
235  // << "\tm_IdScan_LastBrlLayer " << m_IdScan_LastBrlLayer << std::endl;
236 
237  // number of phi neighbours to look at
238  // if ( m_extraPhi.size()==0 )
239  m_extraPhi = std::vector< std::vector<long> >( m_IdScan_MaxNumLayers, std::vector<long>(m_IdScan_MaxNumLayers) );
240 
241  if ( m_fullScanMode ) m_usedROIphiWidth = 2*M_PI;
242  else m_usedROIphiWidth = m_ROIphiWidth;
243 
244  // from IDScanZFinder::initialize
245  m_usedphiBinSize = m_phiBinSize;
246  if ( m_usedphiBinSize < ZFinder_MinPhiSliceSize and ! m_forcePhiBinSize) m_usedphiBinSize = ZFinder_MinPhiSliceSize;
247  if ( m_dphideta > 0 ) m_dphideta *= -m_dphideta;
248 
249  m_invPhiSliceSize = 180./(M_PI*m_usedphiBinSize);
252 
253  for (long l1=0; l1<m_IdScan_MaxNumLayers-1; ++l1) {
254  for (long l2=l1+1; l2<m_IdScan_MaxNumLayers; ++l2) {
255  m_extraPhi[l1][l2]=1; // look at a single phi neighbour
256  }
257  }
258 
260 
261  long first_layer = 0;
262  long offset_layer = 1;
263  if ( m_IdScan_MaxNumLayers<20 ) {
264  first_layer = 1;
265  offset_layer = 0;
266  }
267 
268  long lyrcnt = 0;
269  for (long l1=0; l1<m_IdScan_LastBrlLayer; ++l1) {
270  for (long l2=l1+1; l2<=m_IdScan_LastBrlLayer; ++l2) {
271  double dphi = ZF_dphiBrl[lyrcnt + 7*first_layer];
272  dphi *= m_neighborMultiplier;
273  m_extraPhi[l1][l2]=static_cast<long>( ceil( sqrt(dphi*dphi+ZF_phiRes*ZF_phiRes*2) * m_invPhiSliceSize ) );
274 
275  // std::cout << "test 1 " << l1 << " " << l2 << "\tmax : " << m_IdScan_MaxNumLayers << std::endl;
276 
277  if (m_extraPhi[l1][l2]<1) m_extraPhi[l1][l2]=1;
278  // std::cout << "Delta Phi between layers " << l1 << " and " << l2
279  // << " = "<< ZF_dphiBrl[lyrcnt]
280  // << " rads ( " << m_extraPhi[l1][l2] << " bins including phi resln.)\n";
281  lyrcnt++;
282  }
283  }
284 
285 
287 
288 
289 
290  for ( long lyrpair=12*first_layer ; lyrpair<117; ++lyrpair ) {
291 
292  double dphi = ZF_dphiEC[lyrpair*4+2];
295  long l1 = (long)ZF_dphiEC[lyrpair*4] + offset_layer;
296  long l2 = (long)ZF_dphiEC[lyrpair*4+1] + offset_layer;
297  double eta = ZF_dphiEC[lyrpair*4+3];
298  // std::cout << "Delta Phi between layers " << l1 << " and " << l2
299  // << " = " << dphi << " rads @ eta=" << eta
300  // << ". Extrapolate it to eta=0.9 to get ";
301  dphi = dphi + m_dphideta * ( 0.9 - eta );
302  dphi *= m_neighborMultiplier;
303  m_extraPhi[l1][l2]=static_cast<long>(ceil(sqrt(dphi*dphi+ZF_phiRes*ZF_phiRes*2)*m_invPhiSliceSize));
304 
305  if (m_extraPhi[l1][l2]<1) m_extraPhi[l1][l2]=1;
306 
307  // std::cout << "test 2 " << l1 << " " << l2 << "\tmax : " << m_IdScan_MaxNumLayers << std::endl;
308 
309  // std::cout << dphi << " rads ( " << m_extraPhi[l1][l2] << " bins including phi resln.)\n";
310  }
311 
312 }
313 
314 template<class SpacePoint>
315 double IDScanZFinderInternal<SpacePoint>::computeZV(double r1, double z1, double r2, double z2) const {
316  return (z2*r1-z1*r2)/(r1-r2);
317 }
318 
319 template<class SpacePoint>
320 double IDScanZFinderInternal<SpacePoint>::computeZV(double r1, double p1, double z1,
321  double r2, double p2, double z2) const {
322  double x1, y1, x2, y2;
323  //sincos( p1, &y1, &x1 ); x1 *= r1; y1 *= r1;
324  //sincos( p2, &y2, &x2 ); x2 *= r2; y2 *= r2;
325  x1 = r1 * cos(p1);
326  x2 = r2 * cos(p2);
327  y1 = r1 * sin(p1);
328  y2 = r2 * sin(p2);
329 
330 #define _COMPUTEX0_
331 #ifdef _COMPUTEX0_
332  double slope = (y2-y1)/(x2-x1);
333  double invslope = 1./slope;
334  double x0 = (slope*x1-y1)/(slope+invslope);
335  //double y0 = -1*x0*invslope;
336  double d0sqr = x0*x0*(1+invslope*invslope);
337  // s1 and s2 are the track lengths from the point of closest approach
338  double s1 = sqrt(r1*r1-d0sqr);
339  double s2 = sqrt(r2*r2-d0sqr); // or s1*(x1-x0)/(x2-x0)
340 #else
341  double inv_dels = 1./sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
342  double dotrr = x1*x2 + y1*y2;
343  double s1 = ( dotrr - r1*r1 ) * inv_dels;
344  double s2 = ( r2*r2 - dotrr ) * inv_dels;
345 #endif
346 
347  return (z2*s1-z1*s2)/(s1-s2);
348 }
349 
350 
351 template<class SpacePoint> long IDScanZFinderInternal<SpacePoint>::fillVectors (const std::vector<const SpacePoint *>& spVec,
352  const IRoiDescriptor& roi,
353  std::vector<double>& phi,
354  std::vector<double>& rho,
355  std::vector<double>& zed,
356  std::vector<long>& lyr,
357  std::vector<long>& filledLayers)
358 {
359  std::vector<bool> lcount( m_IdScan_MaxNumLayers, false );
360 
361  // full scan check
362  typename std::vector<const SpacePoint *>::const_iterator SpItr( spVec.begin() );
363 
364  long nSPs = spVec.size();
365 
366  // to shift the phi of space points as if the RoI starts at phi~0
367  // assumes that RoI->phi0() and the SPs are in range [-pi,+pi]
368  //
369  double roiPhiMin, roiPhiMax;
370 
371  if ( m_fullScanMode || roi.isFullscan() ) {
372  roiPhiMin = -M_PI;
373  roiPhiMax = M_PI;
374  }
375  else {
376  // If we trust that all the SPs are properly input, we determine the RoI phi width
377  // using the SPs themselves.
378  // If the RoI phi range is wider than pi, we keep everything as usual.
379  if ( m_trustSPprovider && m_usedROIphiWidth < M_PI )
380  {
381  double roiPhiPosMin( 9.9), roiPhiPosMax(0);
382  double roiPhiNegMin(-9.9), roiPhiNegMax(0); // least negative and most negative
383  for(long i=0; i<nSPs; ++i, ++SpItr)
384  {
385  double spphi = (*SpItr)->phi();
386  if ( spphi>0 && spphi>roiPhiPosMax ) roiPhiPosMax = spphi;
387  if ( spphi>0 && spphi<roiPhiPosMin ) roiPhiPosMin = spphi;
388 
389  if ( spphi<0 && spphi<roiPhiNegMax ) roiPhiNegMax = spphi;
390  if ( spphi<0 && spphi>roiPhiNegMin ) roiPhiNegMin = spphi;
391  }
392 
393  if ( roiPhiNegMax > roiPhiNegMin ) {
394  // if all SPs are in (0, pi):
395  roiPhiMin = roiPhiPosMin;
396  roiPhiMax = roiPhiPosMax;
397  }
398  else if ( roiPhiPosMax < roiPhiPosMin ) {
399  // if all SPs are in (-pi, 0):
400  roiPhiMin = roiPhiNegMax;
401  roiPhiMax = roiPhiNegMin;
402  }
403  else if ( roiPhiPosMin - roiPhiNegMin < M_PI ) {
404  // if we have an RoI that covers phi=0 axis
405  roiPhiMin = roiPhiNegMax;
406  roiPhiMax = roiPhiPosMax;
407  }
408  else {
409  // if we have an RoI that covers phi=pi axis
410  roiPhiMin = roiPhiPosMin;
411  roiPhiMax = roiPhiNegMin;
412  }
413 
414  roiPhiMin -= 1e-10;
415  roiPhiMax += 1e-10;
416 
417  SpItr = spVec.begin(); // rewind the iterator
418  }
419  else {
420 
422  if ( roi.phiMinus()==roi.phiPlus() ) {
423  roiPhiMin = roi.phi()-0.5*m_usedROIphiWidth;
424  roiPhiMax = roi.phi()+0.5*m_usedROIphiWidth;
425  if(roiPhiMin<-M_PI) roiPhiMin+=2*M_PI;
426  if(roiPhiMax>M_PI) roiPhiMax-=2*M_PI;
427  }
428  else {
429  // already wrapped by RoiDescriptor
430  roiPhiMin = roi.phiMinus();
431  roiPhiMax = roi.phiPlus();
432  }
433 
434  }
435 
436  }
437 
438 
439  double dphi = roiPhiMax-roiPhiMin;
440  if ( dphi<0 ) dphi+=2*M_PI;
441 
442  m_usedROIphiWidth = dphi;
443 
444  // std::cout << "m_usedROIphiWidth: " << m_usedROIphiWidth << std::endl;
445  m_NumPhiSlices = long (ceil( m_usedROIphiWidth*m_invPhiSliceSize ));
446 
447 
448  bool piBound=(roiPhiMin>roiPhiMax);
449 
450  int icount = 0;
451 
452  if(!piBound)
453  {
455  for(long i=0; i<nSPs; ++i, ++SpItr)
456  {
457  if (m_pixOnly && !(*SpItr)->isPixel()) continue;
458  double phi2 = (*SpItr)->phi() - roiPhiMin;
459 
460  if ( phi2>=0 && phi2<dphi ) {
461  phi[i] = phi2;
462  rho[i] = (*SpItr)->r();
463  zed[i] = (*SpItr)->z();
464  lyr[i] = (*SpItr)->layer();
465  lcount[lyr[i]]=true;
466  ++icount;
467  }
468  }
469  }
470  else
471  {
473  for(long i=0; i<nSPs; ++i, ++SpItr)
474  {
475  if (m_pixOnly && !(*SpItr)->isPixel()) continue;
476  double phi2 = (*SpItr)->phi() - roiPhiMin;
477  if( phi2<0.0) phi2+=2*M_PI;
478 
479  if ( phi2>=0 && phi2<dphi ) {
480  phi[i] = phi2;
481  rho[i] = (*SpItr)->r();
482  zed[i] = (*SpItr)->z();
483  lyr[i] = (*SpItr)->layer();
484  lcount[lyr[i]]=true;
485  ++icount;
486  }
487  }
488  }
489 
490  if ( icount<nSPs ) {
491 
492  // std::cout << "IDScanZFinderInternal::fillVectors() filtered some spacepoints " << nSPs
493  // << " -> " << icount << std::endl;
495  phi.resize(icount);
496  rho.resize(icount);
497  zed.resize(icount);
498  lyr.resize(icount);
499 
500  nSPs = icount;
501  }
502 
503 
504  // Store in filledLayers the layerNumber of those layers that contain hits.
505  // So, if there are hits in layers 1,3,8 filledLayers[0]=1, filledLayers[1]=3
506  // and filledLayers[2]=8
507  //
508  long filled = 0;
509  for ( long i=0; i<m_IdScan_MaxNumLayers; ++i ) {
510  if ( lcount[i] ) {
511  filledLayers[filled] = i;
512  ++filled;
513  }
514  }
515 
516  // std::cout << "SUTT NSP : " << phi.size() << " " << spVec.size() << std::endl;
517  // for ( unsigned i=0 ; i<phi.size() ; i++ ) {
518  // std::cout << "SUTT SP : " << i << "\tphi " << phi[i] << "\tz " << zed[i] << "\tr " << rho[i] << std::endl;
519  // }
520 
521  return filled;
522 }
523 
524 template<class SpacePoint>
526  const IRoiDescriptor& roi)
527 {
528  std::vector<vertex>* output = new std::vector<vertex>();
529 
530  long nsp = spVec.size();
531  if ( !nsp ) return output; //No points - return nothing
532 
533  SetInternalStatus(0);
534 
535  // Creating vectors of doubles/longs for storing phi,rho,z and layer of space points.
536  // filledLayers is used to know which of all layers contain space points
537  // and fill with relevant info...
538  //
539  std::vector<double> phi(nsp), rho(nsp), zed(nsp);
540  std::vector<long> lyr(nsp), filledLayers(m_IdScan_MaxNumLayers);
541 
542  long filled = this->fillVectors( spVec, roi, phi, rho, zed, lyr, filledLayers);
543 
544  nsp = phi.size();
545 
546  // std::cout << "SUTT roi " << roi << "nsp: " << nsp << std::endl;
547 
548  double zMin = roi.zedMinus();
549  double zMax = roi.zedPlus();
550 
551  // The bin size of the z-histo -and hence the number of bins-
552  // depends on the RoI's |eta| (to account for worsening z-resolution)
553  //
554  const double ZBinSize = m_minZBinSize + m_zBinSizeEtaCoeff*fabs(roi.eta());
555  const double invZBinSize = 1./ZBinSize;
556 
557 
558  const long HalfZhistoBins = long ( ceil( 0.5*(zMax-zMin)*invZBinSize ) );
559  const long NumZhistoBins = 2*HalfZhistoBins;
560 
561  // number of phi bins to look at will get fewer as eta increases
562  const long extraPhiNeg = static_cast<long> ( floor( (0.9 - fabs(roi.eta()))*m_dphideta*m_invPhiSliceSize*m_neighborMultiplier ) );
563 
564 
565 
566  // These are the z-Histograms
567  // Two sets are defined: {n/z}Histo[0][phi][z] will be for positively bending tracks
568  // {n/z}Histo[1][phi][z] will be for negatively bending tracks
569  //
570 
571  long numZhistos = m_zHistoPerPhi ? m_NumPhiSlices : 1 ;
572 
573  // std::vector < std::vector < std::vector<long> > > nHisto( 2, std::vector < std::vector<long> > (numZhistos, std::vector<long> () ) );
574  // std::vector < std::vector < std::vector<double> > > zHisto( 2, std::vector < std::vector<double> > (numZhistos, std::vector<double>() ) );
575 
576  for ( int i=2 ; i-- ; ) { m_nHisto[i].clear(); m_nHisto[i].resize(numZhistos); }
577  for ( int i=2 ; i-- ; ) { m_zHisto[i].clear(); m_zHisto[i].resize(numZhistos); }
578 
579  std::vector< std::vector<long> >* nHisto = m_nHisto;
580  std::vector< std::vector<double> >* zHisto = m_zHisto;
581 
582  m_NMax = 0;
583 
584 
585 
586 
587  //Make a vector of all the PhiSlice instances we need
588  std::vector< PhiSlice* > allSlices( m_NumPhiSlices );
589  for ( unsigned int sliceIndex = 0; sliceIndex < m_NumPhiSlices; sliceIndex++ )
590  {
591  allSlices[ sliceIndex ] = new PhiSlice( sliceIndex, ZBinSize, m_invPhiSliceSize,
592  m_tripletDZ, m_halfTripletDK, m_tripletDP, zMin, zMax,
593  m_IdScan_MaxNumLayers, m_IdScan_LastBrlLayer );
594  }
595 
596  int allSlicesSize = allSlices.size();
597  //Populate the slices
598  for ( long pointIndex = 0; pointIndex < nsp; pointIndex++ )
599  {
600  int phiIndex = floor( phi[ pointIndex ] * m_invPhiSliceSize );
601  if (phiIndex > allSlicesSize) {
602  continue;
603  }
604  allSlices[ phiIndex ]->AddPoint( rho[ pointIndex ], zed[ pointIndex ], phi[ pointIndex ], lyr[ pointIndex ] );
605  }
606 
607  //Read out the slices into flat structures for the whole layers
608  std::vector< std::vector< double > > allLayerRhos, allLayerZs, allLayerPhis;
609  allLayerRhos.resize( m_IdScan_MaxNumLayers );
610  allLayerZs.resize( m_IdScan_MaxNumLayers );
611  allLayerPhis.resize( m_IdScan_MaxNumLayers );
612  std::vector< std::vector< int > > allSliceWidths( m_IdScan_MaxNumLayers, std::vector< int >( m_NumPhiSlices + 1, 0 ) );
613  for ( unsigned int sliceIndex = 0; sliceIndex < m_NumPhiSlices; sliceIndex++ )
614  {
615  allSlices[ sliceIndex ]->MakeWideLayers( &allLayerRhos, &allLayerZs, &allLayerPhis, &allSliceWidths, filled, &filledLayers );
616  }
617 
618  //One histogram per phi slice?
619  if ( m_zHistoPerPhi )
620  {
621  //Allocate all the histograms
622  for ( unsigned int sliceIndex = 0; sliceIndex < m_NumPhiSlices; sliceIndex++ )
623  {
624  nHisto[0][sliceIndex].resize( NumZhistoBins + 1 );
625  zHisto[0][sliceIndex].resize( NumZhistoBins + 1 );
626  }
627  if ( m_chargeAware )
628  {
629  for ( unsigned int sliceIndex = 0; sliceIndex < m_NumPhiSlices; sliceIndex++ )
630  {
631  nHisto[1][sliceIndex].resize( NumZhistoBins + 1 );
632  zHisto[1][sliceIndex].resize( NumZhistoBins + 1 );
633  }
634  }
635 
636  //Populate the histograms
637  for ( unsigned int sliceIndex = 0; sliceIndex < m_NumPhiSlices; sliceIndex++ )
638  {
639  allSlices[ sliceIndex ]->GetHistogram( &( nHisto[0][sliceIndex] ), &( zHisto[0][sliceIndex] ),
640  &( nHisto[1][sliceIndex] ), &( zHisto[1][sliceIndex] ),
641  m_extraPhi, extraPhiNeg, m_nFirstLayers, m_tripletMode, m_chargeAware, nHisto, zHisto );
642  //Note the extra arguments here - pointers to the whole histogram collection
643  //This allows the filling of neighbouring slice histograms as required, but breaks thread safety
644 
645  delete allSlices[ sliceIndex ];
646  }
647  }
648  else
649  {
650  //Allocate the z-axis histograms
651  nHisto[0][0].resize( NumZhistoBins + 1 );
652  zHisto[0][0].resize( NumZhistoBins + 1 );
653  if ( m_chargeAware )
654  {
655  nHisto[1][0].resize( NumZhistoBins + 1 );
656  zHisto[1][0].resize( NumZhistoBins + 1 );
657  }
658 
659  //Populate the histogram - fast and memory-efficient, but not thread safe (use MakeHistogram for thread safety)
660  for ( unsigned int sliceIndex = 0; sliceIndex < m_NumPhiSlices; sliceIndex++ )
661  {
662  allSlices[ sliceIndex ]->GetHistogram( &( nHisto[0][0] ), &( zHisto[0][0] ),
663  &( nHisto[1][0] ), &( zHisto[1][0] ),
664  m_extraPhi, extraPhiNeg, m_nFirstLayers, m_tripletMode, m_chargeAware );
665  delete allSlices[ sliceIndex ];
666  }
667  }
668 
669  /*
670  printf("etaRoI: %f\n", etaRoI);
671  printf("ZBinSize: %f\n", ZBinSize);
672  printf("m_minZBinSize %f\n", m_minZBinSize);
673  printf("NumZhistoBins: %d\n",NumZhistoBins);
674  */
675 
676 
677 
680 
681  double pedestal = 0;
682 
683  if ( m_weightThreshold>0 ) {
684 
685  int count = 0;
686 
687  for ( long zpm=0; zpm<1 || ( m_chargeAware && zpm<2 ) ; ++zpm ) {
688 
689  for(std::vector< std::vector<long> >::iterator nfiter = nHisto[zpm].begin(); nfiter!=nHisto[zpm].end(); ++nfiter) {
690 
691  if((*nfiter).empty()) continue; // only check the filled zHistos
692  if((*nfiter).size() <= 2 ) continue;// this is only a protection : with proper inputs to zfinder, it is always satisfied
693 
694  for(std::vector<long>::iterator niter=nfiter->begin()+2; niter!=nfiter->end(); ++niter ) {
696  if ( *niter>=0 ) {
697  count++;
698  pedestal += *(niter) + *(niter-1) + *(niter-2);
699  }
700  }
701  }
702  }
703 
704  if ( count>0 ) pedestal /= count;
705 
706  }
707 
708 
709 
710 
711  // Now the nHisto's are filled; find the 3 consecutive bins with the highest number of entries...
712  //
713 
714  std::vector<double> zoutput;
715  std::vector<double> woutput;
716 
717 
718 
719  while((int)zoutput.size() < m_numberOfPeaks) {
720 
721  long maxh=0; // was 1 before triplets were introduced
722  long binMax=0;
723  long bending=0, bestPhi=0;
724  long ztest;
725 
726  for ( long zpm=0; zpm<1 || ( m_chargeAware && zpm<2 ) ; ++zpm ) {
727 
728  for(std::vector< std::vector<long> >::iterator nfiter = nHisto[zpm].begin(); nfiter!=nHisto[zpm].end(); ++nfiter) {
729 
730  if((*nfiter).empty()) continue; // only check the filled zHistos
731  if((*nfiter).size() <= 2 ) continue;// this is only a protection : with proper inputs to zfinder, it is always satisfied
732 
733  for(std::vector<long>::iterator niter=(*nfiter).begin()+2; niter!=(*nfiter).end(); ++niter ) {
734 
735  ztest = *(niter-2) + *(niter-1) + *(niter);
736  if ( ztest <= 0 || ztest < maxh ) continue;
739  if ( ztest >= maxh ) { // && ztest>m_weightThreshold ) {
740  long bintest = niter-(*nfiter).begin()-1;
741  if ( ztest > maxh ||
742  // for two candidates at same "height", prefer the more central one
743  (m_preferCentralZ && std::abs( bintest - HalfZhistoBins ) < std::abs( binMax - HalfZhistoBins ) ) ) {
744  maxh = ztest;
745  binMax = bintest;
746  bestPhi = nfiter-nHisto[zpm].begin();
747  bending = zpm;
748  }
749  }
750  }// end of niter loop
751  }
752  }
753  m_NMax = maxh;
756  if ( maxh==0 || ( m_tripletMode==0 && maxh<2 ) ) {
757  break;
758  }
759 
760  // ...and compute the "entries-weighted" average bin position
761 
762  double weightedMax = ( zHisto[bending][bestPhi][binMax] +
763  zHisto[bending][bestPhi][binMax+1] +
764  zHisto[bending][bestPhi][binMax-1] ) /maxh;
765 
767  if ( m_numberOfPeaks>1 ) {
768  nHisto[bending][bestPhi][binMax] = -1;
769  nHisto[bending][bestPhi][binMax-1] = -1;
770  nHisto[bending][bestPhi][binMax+1] = -1;
771  zHisto[bending][bestPhi][binMax] = 0;
772  zHisto[bending][bestPhi][binMax-1] = 0;
773  zHisto[bending][bestPhi][binMax+1] = 0;
774  }
775 
776  int closestVtx = -1; // find the closest vertex already put into the output list
777  float dist2closestVtx = 1000; // should be larger than m_ZFinder_MaxZ*2
778  for ( size_t iv = 0; iv < zoutput.size(); ++iv )
779  if ( fabs(weightedMax-zoutput[iv]) < dist2closestVtx ) {
780  closestVtx = iv;
781  dist2closestVtx = fabs(weightedMax-zoutput[iv]);
782  }
783 
784  if ( dist2closestVtx < m_nvrtxSeparation * ZBinSize ||
785  dist2closestVtx < fabs(weightedMax) * m_vrtxDistCut ) {
786  zoutput[closestVtx] = m_vrtxMixing * weightedMax + (1.0 - m_vrtxMixing) * zoutput[closestVtx] ;
787  woutput[closestVtx] = m_vrtxMixing * maxh + (1.0 - m_vrtxMixing) * woutput[closestVtx] ;
788  } else {
789  zoutput.push_back( weightedMax );
790  woutput.push_back( maxh );
791  }
792 
793  }
794 
795 
800 
808 
811 
812 
813 #if 0
814  if ( m_weightThreshold>0 ) {
815 
823 
824  for ( unsigned i=0 ; i<zoutput.size() ; i++ ) {
827  output->push_back( woutput[i] - pedestal );
828  }
829 
830  }
831  else {
832  for ( unsigned i=0 ; i<zoutput.size() ; i++ ) output->push_back( zoutput[i] );
833  }
834 
835 
836 #else
837 
841  if ( m_weightThreshold>0 ) {
842  for ( unsigned i=0 ; i<zoutput.size() ; i++ ) {
843  output->push_back( vertex( woutput[i] - pedestal, zoutput[i] ) );
844  }
845  }
846  else {
847  for ( unsigned i=0 ; i<zoutput.size() ; i++ ) {
848  output->push_back( vertex( zoutput[i], woutput[i] - pedestal ) );
849  }
850  }
851 
852 
853 #endif
854 
855  // std::cout << "SUTT zoutput size " << zoutput.size() << "\t" << roi << std::endl;
856  // for ( unsigned i=0 ; i<zoutput.size() ; i++ ) std::cout << "SUTT zoutput " << i << "\t" << zoutput[i] << std::endl;
857 
858  return output;
859 }
860 }
861 
862 #endif
863 
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Run1
USAGE: openCoraCool.exe "COOLONL_SCT/COMP200".
Definition: openCoraCool.cxx:57
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
IRoiDescriptor::phi
virtual double phi() const =0
Methods to retrieve data members.
Run1::IDScanZFinderInternal::m_preferCentralZ
bool m_preferCentralZ
Definition: IDScanZFinderInternal.h:140
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
Run1::IDScanZFinderInternal::m_zHistoPerPhi
bool m_zHistoPerPhi
Definition: IDScanZFinderInternal.h:121
Run1::IDScanZFinderInternal::m_extraPhi
std::vector< std::vector< long > > m_extraPhi
Definition: IDScanZFinderInternal.h:126
PhiSlice.h
Run1::IDScanZFinderInternal::m_nFirstLayers
int m_nFirstLayers
Definition: IDScanZFinderInternal.h:135
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
Run1::IDScanZFinderInternal::GetnHisto
const std::vector< std::vector< long > > * GetnHisto()
Definition: IDScanZFinderInternal.h:58
Merge_LAr_blobs.filled
list filled
Definition: Merge_LAr_blobs.py:43
Run1::IDScanZFinderInternal::m_forcePhiBinSize
bool m_forcePhiBinSize
Definition: IDScanZFinderInternal.h:104
Run1::IDScanZFinderInternal::vertex::vertex
vertex(double z, double weight)
Definition: IDScanZFinderInternal.h:44
Run1::IDScanZFinderInternal::m_IdScan_MaxNumLayers
long m_IdScan_MaxNumLayers
maximum number of layers and last barrel layer
Definition: IDScanZFinderInternal.h:94
Run1::IDScanZFinderInternal::m_neighborMultiplier
double m_neighborMultiplier
Definition: IDScanZFinderInternal.h:124
Run1::IDScanZFinderInternal::m_phiBinSize
double m_phiBinSize
Definition: IDScanZFinderInternal.h:103
Run1::IDScanZFinderInternal::m_fullScanMode
bool m_fullScanMode
Definition: IDScanZFinderInternal.h:146
Run1::IDScanZFinderInternal::setLayers
void setLayers(long maxLayers, long lastBarrelLayer)
Definition: IDScanZFinderInternal.h:64
hist_file_dump.d
d
Definition: hist_file_dump.py:137
Run1::IDScanZFinderInternal::getVersion
const std::string & getVersion() const
Definition: IDScanZFinderInternal.h:79
Run1::IDScanZFinderInternal::m_chargeAware
bool m_chargeAware
Definition: IDScanZFinderInternal.h:120
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
Run1::IDScanZFinderInternal::m_tripletDP
double m_tripletDP
Definition: IDScanZFinderInternal.h:152
Run1::IDScanZFinderInternal::m_Type
std::string m_Type
Definition: IDScanZFinderInternal.h:115
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Run1::IDScanZFinderInternal::getName
const std::string & getName() const
Definition: IDScanZFinderInternal.h:78
Run1::IDScanZFinderInternal::vertex::_z
double _z
Definition: IDScanZFinderInternal.h:45
Run1::IDScanZFinderInternal::m_numberOfPeaks
long m_numberOfPeaks
Definition: IDScanZFinderInternal.h:111
Run1::IDScanZFinderInternal::m_ROIphiWidth
double m_ROIphiWidth
Definition: IDScanZFinderInternal.h:106
Run1::IDScanZFinderInternal::~IDScanZFinderInternal
virtual ~IDScanZFinderInternal()
Definition: IDScanZFinderInternal.h:52
Run1::IDScanZFinderInternal::m_Status
int m_Status
Definition: IDScanZFinderInternal.h:118
Run1::IDScanZFinderInternal::vertex::_weight
double _weight
Definition: IDScanZFinderInternal.h:46
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
Run1::IDScanZFinderInternal::m_Name
std::string m_Name
Definition: IDScanZFinderInternal.h:116
Run1::IDScanZFinderInternal::m_tripletDK
double m_tripletDK
Definition: IDScanZFinderInternal.h:150
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
Run1::IDScanZFinderInternal::m_nHisto
std::vector< std::vector< long > > m_nHisto[2]
Definition: IDScanZFinderInternal.h:130
Run1::IDScanZFinderInternal::m_dphideta
double m_dphideta
Definition: IDScanZFinderInternal.h:123
Run1::IDScanZFinderInternal::m_invPhiSliceSize
double m_invPhiSliceSize
Definition: IDScanZFinderInternal.h:100
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
Run1::IDScanZFinderInternal::m_nvrtxSeparation
int m_nvrtxSeparation
Definition: IDScanZFinderInternal.h:139
Run1::IDScanZFinderInternal::m_zBinSizeEtaCoeff
double m_zBinSizeEtaCoeff
Definition: IDScanZFinderInternal.h:109
Run1::IDScanZFinderInternal::computeZV
double computeZV(double r1, double p1, double z1, double r2, double p2, double z2) const
Definition: IDScanZFinderInternal.h:320
Run1::IDScanZFinderInternal::m_usedROIphiWidth
double m_usedROIphiWidth
Definition: IDScanZFinderInternal.h:107
Run1::IDScanZFinderInternal::m_vrtxDistCut
double m_vrtxDistCut
Definition: IDScanZFinderInternal.h:137
IRoiDescriptor::eta
virtual double eta() const =0
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
Run1::IDScanZFinderInternal::IDScanZFinderInternal
IDScanZFinderInternal(const std::string &, const std::string &)
Definition: IDScanZFinderInternal.h:172
skel.l2
l2
Definition: skel.GENtoEVGEN.py:399
ZF_phiRes
const double ZF_phiRes
Definition: ZFinderConstants.h:37
lumiFormat.i
int i
Definition: lumiFormat.py:85
Run1::IDScanZFinderInternal::m_NumPhiSlices
long m_NumPhiSlices
Definition: IDScanZFinderInternal.h:101
z
#define z
Run1::IDScanZFinderInternal::GetzHisto
const std::vector< std::vector< double > > * GetzHisto()
Definition: IDScanZFinderInternal.h:59
IRoiDescriptor
Describes the API of the Region of Ineterest geometry.
Definition: IRoiDescriptor.h:23
Run1::IDScanZFinderInternal::m_usedphiBinSize
double m_usedphiBinSize
Definition: IDScanZFinderInternal.h:105
PhiSlice
Definition: PhiSlice.h:13
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
IRoiDescriptor::phiMinus
virtual double phiMinus() const =0
Run1::IDScanZFinderInternal::getType
const std::string & getType() const
Definition: IDScanZFinderInternal.h:77
Run1::IDScanZFinderInternal::m_vrtxMixing
double m_vrtxMixing
Definition: IDScanZFinderInternal.h:138
Run1::IDScanZFinderInternal::GetNMax
long GetNMax()
Definition: IDScanZFinderInternal.h:61
Run1::IDScanZFinderInternal::m_halfTripletDK
double m_halfTripletDK
Definition: IDScanZFinderInternal.h:151
Run1::IDScanZFinderInternal::GetInternalStatus
int GetInternalStatus() const
Definition: IDScanZFinderInternal.h:81
Run1::IDScanZFinderInternal::initializeInternal
void initializeInternal(long maxLayers, long lastBarrel)
Definition: IDScanZFinderInternal.h:227
Run1::IDScanZFinderInternal::m_trustSPprovider
bool m_trustSPprovider
Definition: IDScanZFinderInternal.h:142
IRoiDescriptor::phiPlus
virtual double phiPlus() const =0
extreme phi values
merge.output
output
Definition: merge.py:17
Run1::IDScanZFinderInternal::m_IdScan_LastBrlLayer
long m_IdScan_LastBrlLayer
Definition: IDScanZFinderInternal.h:95
Run1::IDScanZFinderInternal::SetInternalStatus
int SetInternalStatus(int s)
Definition: IDScanZFinderInternal.h:82
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
Run1::IDScanZFinderInternal::m_tripletDZ
double m_tripletDZ
Definition: IDScanZFinderInternal.h:149
Run1::IDScanZFinderInternal::SetReturnValue
void SetReturnValue(double d)
Definition: IDScanZFinderInternal.h:87
ZFinder_MinPhiSliceSize
const double ZFinder_MinPhiSliceSize
Definition: ZFinderConstants.h:34
eflowRec::phiIndex
unsigned int phiIndex(float phi, float binsize)
calculate phi index for a given phi
Definition: EtaPhiLUT.cxx:23
Run1::IDScanZFinderInternal::m_tripletMode
int m_tripletMode
Definition: IDScanZFinderInternal.h:148
Run1::IDScanZFinderInternal::m_zHisto
std::vector< std::vector< double > > m_zHisto[2]
Definition: IDScanZFinderInternal.h:131
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
IRoiDescriptor::zedPlus
virtual double zedPlus() const =0
the zed and eta values at the most forward and most rear ends of the RoI
Run1::IDScanZFinderInternal::m_pixOnly
bool m_pixOnly
Definition: IDScanZFinderInternal.h:113
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
ZFinderConstants.h
ReadCellNoiseFromCoolCompare.s2
s2
Definition: ReadCellNoiseFromCoolCompare.py:379
IRoiDescriptor::isFullscan
virtual bool isFullscan() const =0
is this a full detector RoI?
IRoiDescriptor::zedMinus
virtual double zedMinus() const =0
IRoiDescriptor.h
Run1::IDScanZFinderInternal::m_weightThreshold
double m_weightThreshold
to apply a hreshold to the found vertex candidates
Definition: IDScanZFinderInternal.h:156
Run1::IDScanZFinderInternal::fillVectors
long fillVectors(const std::vector< const SpacePoint * > &spVec, const IRoiDescriptor &roi, std::vector< double > &phi, std::vector< double > &rho, std::vector< double > &zed, std::vector< long > &lyr, std::vector< long > &filledLayers)
Definition: IDScanZFinderInternal.h:351
Run1::IDScanZFinderInternal::vertex
Definition: IDScanZFinderInternal.h:43
Run1::IDScanZFinderInternal::m_NMax
long m_NMax
Definition: IDScanZFinderInternal.h:133
skel.l1
l1
Definition: skel.GENtoEVGEN.py:398
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
CreatePhysValWebPage.nHisto
nHisto
Definition: CreatePhysValWebPage.py:89
Run1::IDScanZFinderInternal
Definition: IDScanZFinderInternal.h:40
Run1::IDScanZFinderInternal::m_minZBinSize
double m_minZBinSize
Definition: IDScanZFinderInternal.h:108
Run1::IDScanZFinderInternal::computeZV
double computeZV(double r1, double z1, double r2, double z2) const
Definition: IDScanZFinderInternal.h:315
Run1::IDScanZFinderInternal::m_returnval
double m_returnval
Definition: IDScanZFinderInternal.h:144
fitman.rho
rho
Definition: fitman.py:532
Run1::IDScanZFinderInternal::GetReturnValue
double GetReturnValue() const
Definition: IDScanZFinderInternal.h:88
Run1::IDScanZFinderInternal::findZInternal
std::vector< vertex > * findZInternal(const std::vector< const SpacePoint * > &spVec, const IRoiDescriptor &roi)
Definition: IDScanZFinderInternal.h:525