CMS 3D CMS Logo
List of all members | Public Types | Public Member Functions | Private Attributes
TrackProducerFP420 Class Reference

#include <TrackProducerFP420.h>

Public Types

typedef std::vector< ClusterFP420 >::const_iterator ClusterFP420Iter
 

Public Member Functions

std::vector< TrackFP420 > trackFinderSophisticated (edm::Handle< ClusterCollectionFP420 > input, int det) const
 
 TrackProducerFP420 (int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, bool, bool, bool, bool, double, double, float, float, double, int, double, double)
 

Private Attributes

float chiCutX
 
float chiCutY
 
double dXX
 
double dYY
 
double gapBlade
 
double pitchX
 
double pitchXW
 
double pitchY
 
double pitchYW
 
int pn0
 
std::vector< TrackFP420 > rhits
 
int rn0
 
int sn0
 
ClusterCollectionFP420 soutput
 
bool UseHalfPitchShiftInX
 
bool UseHalfPitchShiftInXW
 
bool UseHalfPitchShiftInY
 
bool UseHalfPitchShiftInYW
 
int verbos
 
double XsensorSize
 
int xytype
 
double YsensorSize
 
double z420
 
double zBlade
 
double zD2
 
double zD3
 
double ZGapLDet
 
double zinibeg
 
double ZSiDet
 
double ZSiPlane
 
double ZSiStep
 

Detailed Description

Definition at line 18 of file TrackProducerFP420.h.

Member Typedef Documentation

◆ ClusterFP420Iter

typedef std::vector<ClusterFP420>::const_iterator TrackProducerFP420::ClusterFP420Iter

Definition at line 20 of file TrackProducerFP420.h.

Constructor & Destructor Documentation

◆ TrackProducerFP420()

TrackProducerFP420::TrackProducerFP420 ( int  asn0,
int  apn0,
int  arn0,
int  axytype,
double  az420,
double  azD2,
double  azD3,
double  apitchX,
double  apitchY,
double  apitchXW,
double  apitchYW,
double  aZGapLDet,
double  aZSiStep,
double  aZSiPlane,
double  aZSiDet,
double  azBlade,
double  agapBlade,
bool  aUseHalfPitchShiftInX,
bool  aUseHalfPitchShiftInY,
bool  aUseHalfPitchShiftInXW,
bool  aUseHalfPitchShiftInYW,
double  adXX,
double  adYY,
float  achiCutX,
float  achiCutY,
double  azinibeg,
int  verbosity,
double  aXsensorSize,
double  aYsensorSize 
)

Definition at line 22 of file TrackProducerFP420.cc.

References gather_cfg::cout, ctppsCommonDQMSource_cfi::verbosity, FP420Track_cfi::z420, FP420Track_cfi::zD2, and FP420Track_cfi::zD3.

50  {
51  //
52  // Everything that depend on the det
53  //
54  verbos = verbosity;
55  sn0 = asn0;
56  pn0 = apn0;
57  rn0 = arn0;
58  xytype = axytype;
59  z420 = az420;
60  zD2 = azD2;
61  zD3 = azD3;
62  //zUnit= azUnit;
63  pitchX = apitchX;
64  pitchY = apitchY;
65  pitchXW = apitchXW;
66  pitchYW = apitchYW;
67  ZGapLDet = aZGapLDet;
68  ZSiStep = aZSiStep;
69  ZSiPlane = aZSiPlane;
70  ZSiDet = aZSiDet;
71  zBlade = azBlade;
72  gapBlade = agapBlade;
73 
74  UseHalfPitchShiftInX = aUseHalfPitchShiftInX;
75  UseHalfPitchShiftInY = aUseHalfPitchShiftInY;
76  UseHalfPitchShiftInXW = aUseHalfPitchShiftInXW;
77  UseHalfPitchShiftInYW = aUseHalfPitchShiftInYW;
78  dXX = adXX;
79  dYY = adYY;
80  chiCutX = achiCutX;
81  chiCutY = achiCutY;
82  zinibeg = azinibeg;
83  XsensorSize = aXsensorSize;
84  YsensorSize = aYsensorSize;
85 
86  if (verbos > 0) {
87  std::cout << "TrackProducerFP420: call constructor" << std::endl;
88  std::cout << " sn0= " << sn0 << " pn0= " << pn0 << " rn0= " << rn0 << " xytype= " << xytype << std::endl;
89  std::cout << " zD2= " << zD2 << " zD3= " << zD3 << " zinibeg= " << zinibeg << std::endl;
90  //std::cout << " zUnit= " << zUnit << std::endl;
91  std::cout << " pitchX= " << pitchX << " pitchY= " << pitchY << std::endl;
92  std::cout << " ZGapLDet= " << ZGapLDet << std::endl;
93  std::cout << " ZSiStep= " << ZSiStep << " ZSiPlane= " << ZSiPlane << std::endl;
94  std::cout << " ZSiDet= " << ZSiDet << std::endl;
95  std::cout << " UseHalfPitchShiftInX= " << UseHalfPitchShiftInX << " UseHalfPitchShiftInY= " << UseHalfPitchShiftInY
96  << std::endl;
97  std::cout << "TrackProducerFP420:----------------------" << std::endl;
98  std::cout << " dXX= " << dXX << " dYY= " << dYY << std::endl;
99  std::cout << " chiCutX= " << chiCutX << " chiCutY= " << chiCutY << std::endl;
100  }
101 }
ctppsCommonDQMSource_cfi.verbosity
verbosity
Definition: ctppsCommonDQMSource_cfi.py:14
gather_cfg.cout
cout
Definition: gather_cfg.py:144

Member Function Documentation

◆ trackFinderSophisticated()

std::vector< TrackFP420 > TrackProducerFP420::trackFinderSophisticated ( edm::Handle< ClusterCollectionFP420 >  input,
int  det 
) const

Definition at line 105 of file TrackProducerFP420.cc.

References funct::abs(), ClusterFP420::amplitudes(), ClusterFP420::barycenter(), ClusterFP420::barycenterW(), ClusterFP420::barycerror(), ClusterFP420::barycerrorW(), haddnano::cl, GetRecoTauVFromDQM_MC_cff::cl2, gather_cfg::cout, cuy::ii, input, nlayers, FP420NumberingScheme::packMYIndex(), multPhiCorr_741_25nsDY_cfi::py, FP420NumberingScheme::realzside(), corrVsCorr::ry, submitPVValidationJobs::t, RandomServiceHelper::t1, RandomServiceHelper::t2, FP420NumberingScheme::unpackCopyIndex(), FP420NumberingScheme::unpackLayerIndex(), FP420NumberingScheme::unpackOrientation(), FP420Track_cfi::z420, FP420Track_cfi::zD2, FP420Track_cfi::zD3, SiStripMonitorCluster_cfi::zmax, and ecaldqm::zside().

106  {
108 
109  std::vector<TrackFP420> rhits;
110  int constexpr restracks = 10; // max # tracks
111  rhits.reserve(restracks);
112  rhits.clear();
113  double Ax[restracks];
114  double Bx[restracks];
115  double Cx[restracks];
116  int Mx[restracks];
117  double Ay[restracks];
118  double By[restracks];
119  double Cy[restracks];
120  int My[restracks];
121  double AxW[restracks];
122  double BxW[restracks];
123  double CxW[restracks];
124  int MxW[restracks];
125  double AyW[restracks];
126  double ByW[restracks];
127  double CyW[restracks];
128  int MyW[restracks];
129  if (verbos > 0) {
130  std::cout << "===============================================================================" << std::endl;
131  std::cout << "=================================================================" << std::endl;
132  std::cout << "==========================================================" << std::endl;
133  std::cout << "= =" << std::endl;
134  std::cout << "TrackProducerFP420: Start trackFinderSophisticated " << std::endl;
135  }
138  // xytype is the sensor grid arrangment
139  if (xytype < 1 || xytype > 2) {
140  std::cout << "TrackProducerFP420:ERROR in trackFinderSophisticated: check xytype = " << xytype << std::endl;
141  return rhits;
142  }
143  // sn0= 3 - 2St configuration, sn0= 4 - 3St configuration
144  // if( sn0 < 3 || sn0 > 4 ){
145  if (sn0 != 3) {
146  std::cout << "TrackProducerFP420:ERROR in trackFinderSophisticated: check sn0 (configuration) = " << sn0
147  << std::endl;
148  return rhits;
149  }
151  int constexpr zbeg = 1, zmax = 3; // means layer 1 and 2 in superlayer, i.e. for loop: 1,2
153  // .
154  int constexpr reshits1 = 12; // is max # cl in sensor of copyinlayer=1
155  int constexpr reshits2 = 24; // (reshits2-reshits1) is max # cl in sensors of copyinlayers= 2 or 3
156  int constexpr resplanes = 20; // Number of planes
157  int nX[resplanes] = {}, nY[resplanes] = {}; // #cl for every X and Y plane
158  int uX[resplanes] = {}, uY[resplanes] = {}; // #cl used for every X and Y plane
159  double zX[reshits2][resplanes], xX[reshits2][resplanes], wX[reshits2][resplanes];
160  double zY[reshits2][resplanes], yY[reshits2][resplanes], wY[reshits2][resplanes];
161  double yXW[reshits2][resplanes], wXW[reshits2][resplanes];
162  double xYW[reshits2][resplanes], wYW[reshits2][resplanes];
163  bool qX[reshits2][resplanes], qY[reshits2][resplanes];
164  // .
165  int txf = 0;
166  int txs1 = 0;
167  int txss = 0;
168  int tyf = 0;
169  int tys1 = 0;
170  int tyss = 0;
171  // .
172  double pitch = 0.;
173  double pitchW = 0.;
174  if (xytype == 1) {
175  pitch = pitchY;
176  pitchW = pitchYW;
177  } else if (xytype == 2) {
178  pitch = pitchX;
179  pitchW = pitchXW;
180  }
181 
182  //current change of geometry:
183  float Xshift = pitch / 2.;
184  float Yshift = pitchW / 2.;
185 
186  //
187  int nmetcurx = 0;
188  int nmetcury = 0;
189  unsigned int ii0 = 999999;
190  int allplacesforsensors = 7;
191  for (int sector = 1; sector < sn0; sector++) {
192  for (int zmodule = 1; zmodule < pn0; zmodule++) {
193  for (int zsideinorder = 1; zsideinorder < allplacesforsensors; zsideinorder++) {
194  int zside = FP420NumberingScheme::realzside(rn0, zsideinorder); // 1, 3, 5, 2, 4, 6
195  if (verbos == -49) {
196  std::cout << "TrackProducerFP420: sector= " << sector << " zmodule= " << zmodule
197  << " zsideinorder= " << zsideinorder << " zside= " << zside << " det= " << det << std::endl;
198  }
199  if (zside != 0) {
200  int justlayer = FP420NumberingScheme::unpackLayerIndex(rn0, zside); // 1, 2
201  if (justlayer < 1 || justlayer > 2) {
202  std::cout << "TrackProducerFP420:WRONG justlayer= " << justlayer << std::endl;
203  }
204  int copyinlayer = FP420NumberingScheme::unpackCopyIndex(rn0, zside); // 1, 2, 3
205  if (copyinlayer < 1 || copyinlayer > 3) {
206  std::cout << "TrackProducerFP420:WRONG copyinlayer= " << copyinlayer << std::endl;
207  }
208  int orientation = FP420NumberingScheme::unpackOrientation(rn0, zside); // Front: = 1; Back: = 2
209  if (orientation < 1 || orientation > 2) {
210  std::cout << "TrackProducerFP420:WRONG orientation= " << orientation << std::endl;
211  }
212  // ii is a continues numbering of planes(!) over two arm FP420 set up
213  // and ...[ii] massives have prepared in use of ii
214  int detfixed =
215  1; // use this treatment for each set up arm, hence no sense to repete the same for +FP420 and -FP420;
216  int nlayers = 3; // 2=3-1
217  unsigned int ii =
218  FP420NumberingScheme::packMYIndex(nlayers, pn0, sn0, detfixed, justlayer, sector, zmodule) -
219  1; // substruct 1 from 1(+1), 2(+2), 3(+3),4(+4),5...,6...,7...,8...,9...,10... (1st Station) ,11...,12...,13,...20... (2nd Station)
220  // ii = 0-19 --> 20 items
222 
223  if (verbos == -49) {
224  std::cout << "TrackProducerFP420: justlayer= " << justlayer << " copyinlayer= " << copyinlayer
225  << " ii= " << ii << std::endl;
226  }
227 
228  double zdiststat = 0.;
229  if (sn0 < 4) {
230  if (sector == 2)
231  zdiststat = zD3;
232  } else {
233  if (sector == 2)
234  zdiststat = zD2;
235  if (sector == 3)
236  zdiststat = zD3;
237  }
238  double kplane = -(pn0 - 1) / 2 - 0.5 + (zmodule - 1); //-3.5 +0...5 = -3.5,-2.5,-1.5,+2.5,+1.5
239  double zcurrent = zinibeg + z420 + (ZSiStep - ZSiPlane) / 2 + kplane * ZSiStep + zdiststat;
240  //double zcurrent = zinibeg +(ZSiStep-ZSiPlane)/2 + kplane*ZSiStep + (sector-1)*zUnit;
241 
242  if (justlayer == 1) {
243  if (orientation == 1)
244  zcurrent += (ZGapLDet + ZSiDet / 2);
245  else if (orientation == 2)
246  zcurrent += zBlade - (ZGapLDet + ZSiDet / 2);
247  } else if (justlayer == 2) {
248  if (orientation == 1)
249  zcurrent += (ZGapLDet + ZSiDet / 2) + zBlade + gapBlade;
250  else if (orientation == 2)
251  zcurrent += 2 * zBlade + gapBlade - (ZGapLDet + ZSiDet / 2);
252  }
253  // .
254  //
255  if (det == 2)
256  zcurrent = -zcurrent;
257  //
258  //
259  // .
260  // local - global systems with possible shift of every second plate:
261 
262  // for xytype=1
263  float dYYcur = dYY; // XSiDet/2.
264  float dYYWcur = dXX; //(BoxYshft+dYGap) + (YSi - YSiDet)/2. = 4.7
265  // for xytype=2
266  float dXXcur = dXX; //(BoxYshft+dYGap) + (YSi - YSiDet)/2. = 4.7
267  float dXXWcur = dYY; // XSiDet/2.
268  // .
269  if (justlayer == 2) {
270  // X-type: x-coord
271  if (UseHalfPitchShiftInX) {
272  dXXcur += Xshift;
273  }
274  // X-type: y-coord
275  if (UseHalfPitchShiftInXW) {
276  dXXWcur -= Yshift;
277  }
278  }
279  //
280  double XXXDelta = 0.0;
281  double YYYDelta = 0.0;
282  if (copyinlayer == 2) {
283  XXXDelta = XsensorSize;
284  YYYDelta = XsensorSize;
285  } else if (copyinlayer == 3) {
286  XXXDelta = 2. * XsensorSize;
287  YYYDelta = 2. * XsensorSize;
288  }
289 
290  // .
291  // GET CLUSTER collection !!!!
292  // .
293  unsigned int iu = FP420NumberingScheme::packMYIndex(rn0, pn0, sn0, det, zside, sector, zmodule);
294  if (verbos > 0) {
295  std::cout << "TrackProducerFP420: check iu = " << iu << std::endl;
296  std::cout << "TrackProducerFP420: sector= " << sector << " zmodule= " << zmodule << " zside= " << zside
297  << " det= " << det << " rn0= " << rn0 << " pn0= " << pn0 << " sn0= " << sn0
298  << " copyinlayer= " << copyinlayer << std::endl;
299  }
300  //============================================================================================================ put into currentclust
301  std::vector<ClusterFP420> currentclust;
302  currentclust.clear();
303  ClusterCollectionFP420::Range outputRange;
304  outputRange = input->get(iu);
305  // fill output in currentclust vector (for may be sorting? or other checks)
306  ClusterCollectionFP420::ContainerIterator sort_begin = outputRange.first;
307  ClusterCollectionFP420::ContainerIterator sort_end = outputRange.second;
308  for (; sort_begin != sort_end; ++sort_begin) {
309  // std::sort(currentclust.begin(),currentclust.end());
310  currentclust.push_back(*sort_begin);
311  } // for
312 
313  if (verbos > 0) {
314  std::cout << "TrackProducerFP420: currentclust.size = " << currentclust.size() << std::endl;
315  }
316  //============================================================================================================
317 
318  std::vector<ClusterFP420>::const_iterator simHitIter = currentclust.begin();
319  std::vector<ClusterFP420>::const_iterator simHitIterEnd = currentclust.end();
320 
321  if (xytype == 1) {
322  if (ii != ii0) {
323  ii0 = ii;
324  nY[ii] = 0; // # cl in every Y plane (max is reshits)
325  uY[ii] = 0; // current used # cl in every X plane
326  nmetcury = 0;
327  }
328  } else if (xytype == 2) {
329  if (ii != ii0) {
330  ii0 = ii;
331  nX[ii] = 0; // # cl in every X plane (max is reshits)
332  uX[ii] = 0; // current used # cl in every X plane
333  nmetcurx = 0;
334  }
335  }
336  // loop in #clusters of current sensor
337  for (; simHitIter != simHitIterEnd; ++simHitIter) {
338  const ClusterFP420 icluster = *simHitIter;
339 
340  // fill vectors for track reconstruction
341 
342  //disentangle complicated pattern recognition of hits?
343  // Y:
344  if (xytype == 1) {
345  nY[ii]++;
346  if (copyinlayer == 1 && nY[ii] > reshits1) {
347  nY[ii] = reshits1;
348  std::cout << "WARNING-ERROR:TrackproducerFP420: currentclust.size()= " << currentclust.size()
349  << " bigger reservated number of hits"
350  << " zcurrent=" << zY[nY[ii] - 1][ii] << " copyinlayer= " << copyinlayer << " ii= " << ii
351  << std::endl;
352  } else if (copyinlayer != 1 && nY[ii] > reshits2) {
353  nY[ii] = reshits2;
354  std::cout << "WARNING-ERROR:TrackproducerFP420: currentclust.size()= " << currentclust.size()
355  << " bigger reservated number of hits"
356  << " zcurrent=" << zY[nY[ii] - 1][ii] << " copyinlayer= " << copyinlayer << " ii= " << ii
357  << std::endl;
358  }
359  zY[nY[ii] - 1][ii] = zcurrent;
360  yY[nY[ii] - 1][ii] = icluster.barycenter() * pitch + 0.5 * pitch + YYYDelta;
361  xYW[nY[ii] - 1][ii] = icluster.barycenterW() * pitchW + 0.5 * pitchW;
362  // go to global system:
363  yY[nY[ii] - 1][ii] = yY[nY[ii] - 1][ii] - dYYcur;
364  wY[nY[ii] - 1][ii] =
365  1. / (icluster.barycerror() * pitch); //reciprocal of the variance for each datapoint in y
366  wY[nY[ii] - 1][ii] *= wY[nY[ii] - 1][ii]; //reciprocal of the variance for each datapoint in y
367  if (det == 2) {
368  xYW[nY[ii] - 1][ii] = (xYW[nY[ii] - 1][ii] + dYYWcur);
369  } else {
370  xYW[nY[ii] - 1][ii] = -(xYW[nY[ii] - 1][ii] + dYYWcur);
371  }
372  wYW[nY[ii] - 1][ii] =
373  1. / (icluster.barycerrorW() * pitchW); //reciprocal of the variance for each datapoint in y
374  wYW[nY[ii] - 1][ii] *= wYW[nY[ii] - 1][ii]; //reciprocal of the variance for each datapoint in y
375  qY[nY[ii] - 1][ii] = true;
376  if (copyinlayer == 1 && nY[ii] == reshits1)
377  break;
378  if (copyinlayer != 1 && nY[ii] == reshits2)
379  break;
380  }
381  // X:
382  else if (xytype == 2) {
383  nX[ii]++;
384  if (verbos == -49) {
385  std::cout << "TrackproducerFP420: nX[ii]= " << nX[ii] << " Ncl= " << currentclust.size()
386  << " copyinlayer= " << copyinlayer << " ii= " << ii << " zcurrent = " << zcurrent
387  << " xX= " << icluster.barycenter() * pitch + 0.5 * pitch + XXXDelta
388  << " yXW= " << icluster.barycenterW() * pitchW + 0.5 * pitchW << " det= " << det
389  << " cl.size= " << icluster.amplitudes().size() << " cl.ampl[0]= " << icluster.amplitudes()[0]
390  << std::endl;
391  }
392  if (copyinlayer == 1 && nX[ii] > reshits1) {
393  std::cout << "WARNING-ERROR:TrackproducerFP420: nX[ii]= " << nX[ii]
394  << " bigger reservated number of hits"
395  << " currentclust.size()= " << currentclust.size() << " copyinlayer= " << copyinlayer
396  << " ii= " << ii << std::endl;
397  nX[ii] = reshits1;
398  } else if (copyinlayer != 1 && nX[ii] > reshits2) {
399  std::cout << "WARNING-ERROR:TrackproducerFP420: nX[ii]= " << nX[ii]
400  << " bigger reservated number of hits"
401  << " currentclust.size()= " << currentclust.size() << " copyinlayer= " << copyinlayer
402  << " ii= " << ii << std::endl;
403  nX[ii] = reshits2;
404  }
405  zX[nX[ii] - 1][ii] = zcurrent;
406  xX[nX[ii] - 1][ii] = icluster.barycenter() * pitch + 0.5 * pitch + XXXDelta;
407  yXW[nX[ii] - 1][ii] = icluster.barycenterW() * pitchW + 0.5 * pitchW;
408  // go to global system:
409  xX[nX[ii] - 1][ii] = -(xX[nX[ii] - 1][ii] + dXXcur);
410  wX[nX[ii] - 1][ii] =
411  1. / (icluster.barycerror() * pitch); //reciprocal of the variance for each datapoint in y
412  wX[nX[ii] - 1][ii] *= wX[nX[ii] - 1][ii]; //reciprocal of the variance for each datapoint in y
413  if (det == 2) {
414  yXW[nX[ii] - 1][ii] = -(yXW[nX[ii] - 1][ii] - dXXWcur);
415  } else {
416  yXW[nX[ii] - 1][ii] = yXW[nX[ii] - 1][ii] - dXXWcur;
417  }
418  wXW[nX[ii] - 1][ii] =
419  1. / (icluster.barycerrorW() * pitchW); //reciprocal of the variance for each datapoint in y
420  wXW[nX[ii] - 1][ii] *= wXW[nX[ii] - 1][ii]; //reciprocal of the variance for each datapoint in y
421  qX[nX[ii] - 1][ii] = true;
422  if (verbos == -29) {
423  std::cout << "trackFinderSophisticated: nX[ii]= " << nX[ii] << " ii = " << ii
424  << " zcurrent = " << zcurrent << " yXW[nX[ii]-1][ii] = " << yXW[nX[ii] - 1][ii]
425  << " xX[nX[ii]-1][ii] = " << xX[nX[ii] - 1][ii] << std::endl;
426  std::cout << " XXXDelta= " << XXXDelta << " dXXcur= " << dXXcur << " -dXXWcur= " << -dXXWcur
427  << std::endl;
428  std::cout << " icluster.barycerrorW()*pitchW= " << icluster.barycerrorW() * pitchW
429  << " wXW[nX[ii]-1][ii]= " << wXW[nX[ii] - 1][ii] << std::endl;
430  std::cout << " -icluster.barycenterW()*pitchW+0.5*pitchW = "
431  << icluster.barycenterW() * pitchW + 0.5 * pitchW << std::endl;
432  std::cout << "============================================================" << std::endl;
433  }
434  if (verbos > 0) {
435  std::cout << "trackFinderSophisticated: nX[ii]= " << nX[ii] << " ii = " << ii
436  << " zcurrent = " << zcurrent << " xX[nX[ii]-1][ii] = " << xX[nX[ii] - 1][ii] << std::endl;
437  std::cout << " wX[nX[ii]-1][ii] = " << wX[nX[ii] - 1][ii]
438  << " wXW[nX[ii]-1][ii] = " << wXW[nX[ii] - 1][ii] << std::endl;
439  std::cout << " -icluster.barycenter()*pitch-0.5*pitch = "
440  << -icluster.barycenter() * pitch - 0.5 * pitch << " -dXXcur = " << -dXXcur
441  << " -XXXDelta = " << -XXXDelta << std::endl;
442  std::cout << "============================================================" << std::endl;
443  }
444 
445  if (copyinlayer == 1 && nX[ii] == reshits1)
446  break;
447  if (copyinlayer != 1 && nX[ii] == reshits2)
448  break;
449  } // if(xytype
450 
451  } // for loop in #clusters (can be breaked)
452 
453  // Y:
454  if (xytype == 1) {
455  if (nY[ii] > nmetcury) { /* # Y-planes w/ clusters */
456  nmetcury = nY[ii];
457  ++tyf;
458  if (sector == 1)
459  ++tys1;
460  if (sector == (sn0 - 1))
461  ++tyss;
462  }
463  }
464  // X:
465  else if (xytype == 2) {
466  if (nX[ii] > nmetcurx) { /* # X-planes w/ clusters */
467  nmetcurx = nX[ii];
468  ++txf;
469  if (sector == 1)
470  ++txs1;
471  if (sector == (sn0 - 1))
472  ++txss;
473  }
474  }
475  //================================== end of for loops in continuius number iu:
476  } //if(zside!=0
477  } // for zsideinorder
478  } // for zmodule
479  } // for sector
480  if (verbos > 0) {
481  std::cout << "trackFinderSophisticated: tyf= " << tyf << " tys1 = " << tys1 << " tyss = " << tyss << std::endl;
482  std::cout << "trackFinderSophisticated: txf= " << txf << " txs1 = " << txs1 << " txss = " << txss << std::endl;
483  std::cout << "============================================================" << std::endl;
484  }
485 
486  //===========================================================================================================================
487  //===========================================================================================================================
488  //===========================================================================================================================
489  //====================== start road finder =============================================================================
490  //===========================================================================================================================
491 
492  // max # iterations to find track using go over of different XZ and YZ fits to find the good chi2X and chi2Y simultaneously(!!!)
493  int constexpr nitMax = 10;
494 
495  // criteria for track selection:
496  // track is selected if for 1st station #cl >=pys1Cut
497  // int pys1Cut = 5, pyssCut = 5, pyallCut=12;
498  // int pys1Cut = 1, pyssCut = 1, pyallCut= 3;
499 
500  // int pys1Cut = 3, pyssCut = 3, pyallCut= 6; // before geom. update
501  // int pys1Cut = 2, pyssCut = 2, pyallCut= 4; // bad for 5 layers per station
502  int constexpr pys1Cut = 3, pyssCut = 1, pyallCut = 5;
503 
504  // double yyyvtx = 0.0, xxxvtx = -15; //mm
505 
507  //
508  // for configuration: 3St, 1m for 1-2 St:
509  // double sigman=0.1, ssigma = 1.0, sigmam=0.15;/* ssigma is foreseen to match 1st point of 2nd Station*/
510  //
511  // for equidistant 3 Stations:
512  //
513  // for tests:
514  // double sigman=118., ssigma = 299., sigmam=118.;
515  // RMS1=0.013, RMS2 = 1.0, RMS3 = 0.018 see plots d1XCL, d2XCL, d3XCL
516  //
517  // double sigman=0.05, ssigma = 2.5, sigmam=0.06;
518  // double sigman=0.18, ssigma = 1.8, sigmam=0.18;
519  // double sigman=0.18, ssigma = 2.9, sigmam=0.18;
520  //
522  // for 3 Stations:
523  // LAST:
524  double sigman = 0.18, ssigma = 2.5, sigmam = 0.18;
525  if (sn0 < 4) {
526  // for 2 Stations:
527  // sigman=0.24, ssigma = 4.2, sigmam=0.33;
528  // sigman=0.18, ssigma = 3.9, sigmam=0.18;
529  // sigman=0.18, ssigma = 3.6, sigmam=0.18;
530  //
531 
532  //
533 
534  // sigman=0.18, ssigma = 3.3, sigmam=0.18;// before geometry update for 4 sensors per superlayer
535  // sigman=0.30, ssigma = 7.1, sigmam=0.40;// for update
536  sigman = 0.30, ssigma = 8.0,
537  sigmam = 1.0; // for matching update to find point nearby to fit track in 1st plane of 2nd Station
538  //
539  }
540  if (verbos > 0) {
541  std::cout << "trackFinderSophisticated: ssigma= " << ssigma << std::endl;
542  }
543 
545 
546  /* ssigma = 3. * 8000.*(0.025+0.009)/sqrt(pn0-1)/100. = 2.9 mm(!!!)----
547  ssigma is reduced by factor k_reduced = (sn0-1)-sector+1 = sn0-sector
548  # Stations currentStation
549  2Stations: sector=2, sn0=3 , sn0-sector = 1 --> k_reduced = 1
550  3Stations: sector=2, sn0=4 , sn0-sector = 2 --> k_reduced = 2
551  3Stations: sector=3, sn0=4 , sn0-sector = 1 --> k_reduced = 1
552  */
553  int numberXtracks = 0, numberYtracks = 0, totpl = 2 * (pn0 - 1) * (sn0 - 1);
554  double sigma;
555 
556  for (int xytypecurrent = xytype; xytypecurrent < xytype + 1; ++xytypecurrent) {
557  if (verbos > 0) {
558  std::cout << "trackFinderSophisticated: xytypecurrent= " << xytypecurrent << std::endl;
559  }
560 
561  //
562  //
563  double tg0 = 0.;
564  int qAcl[resplanes], qAii[resplanes], fip = 0, niteration = 0;
565  int ry = 0, rys1 = 0, ryss = 0;
566  int tas1 = tys1, tass = tyss, taf = tyf;
567  bool SelectTracks = true;
568  //
570  // .
571 
572  double yA[reshits2][resplanes], zA[reshits2][resplanes], wA[reshits2][resplanes];
573  int nA[resplanes], uA[resplanes];
574  bool qA[reshits2][resplanes];
575  //
576  // Y:
577  //====================== start road finder for xytypecurrent = 1 ===========================================================
578  if (xytypecurrent == 1) {
579  //===========================================================================================================================
580  numberYtracks = 0;
581  tg0 = 3 * 1. / (800. + 20.); // for Y: 1cm/... *3 - 3sigma range
582  tas1 = tys1;
583  tass = tyss;
584  taf = tyf;
585  for (int ii = 0; ii < totpl; ++ii) {
586  if (verbos > 0) {
587  std::cout << "trackFinderSophisticated: ii= " << ii << " nY[ii]= " << nY[ii] << std::endl;
588  std::cout << "trackFinderSophisticated: ii= " << ii << " uY[ii]= " << uY[ii] << std::endl;
589  }
590  nA[ii] = nY[ii];
591  uA[ii] = uY[ii];
592  for (int cl = 0; cl < nA[ii]; ++cl) {
593  if (verbos > 0) {
594  std::cout << " cl= " << cl << " yY[cl][ii]= " << yY[cl][ii] << std::endl;
595  std::cout << " zY[cl][ii]= " << zY[cl][ii] << " wY[cl][ii]= " << wY[cl][ii] << " qY[cl][ii]= " << qY[cl][ii]
596  << std::endl;
597  }
598  yA[cl][ii] = yY[cl][ii];
599  zA[cl][ii] = zY[cl][ii];
600  wA[cl][ii] = wY[cl][ii];
601  qA[cl][ii] = qY[cl][ii];
602  }
603  }
604  //===========================================================================================================================
605  } // if xytypecurrent ==1
606  // X:
607  //====================== start road finder for superlayer = 2 ===========================================================
608  else if (xytypecurrent == 2) {
609  //===========================================================================================================================
610  numberXtracks = 0;
611  tg0 = 3 * 2. / (800. + 20.); // for X: 2cm/... *3 - 3sigma range
612  tas1 = txs1;
613  tass = txss;
614  taf = txf;
615  for (int ii = 0; ii < totpl; ++ii) {
616  if (verbos > 0) {
617  std::cout << "trackFinderSophisticated: ii= " << ii << " nX[ii]= " << nX[ii] << std::endl;
618  std::cout << "trackFinderSophisticated: ii= " << ii << " uX[ii]= " << uX[ii] << std::endl;
619  }
620  nA[ii] = nX[ii];
621  uA[ii] = uX[ii];
622  for (int cl = 0; cl < nA[ii]; ++cl) {
623  if (verbos == -29) {
624  std::cout << " cl= " << cl << " xX[cl][ii]= " << xX[cl][ii] << std::endl;
625  std::cout << " zX[cl][ii]= " << zX[cl][ii] << " wX[cl][ii]= " << wX[cl][ii] << " qX[cl][ii]= " << qX[cl][ii]
626  << std::endl;
627  }
628  yA[cl][ii] = xX[cl][ii];
629  zA[cl][ii] = zX[cl][ii];
630  wA[cl][ii] = wX[cl][ii];
631  qA[cl][ii] = qX[cl][ii];
632  }
633  }
634  //===========================================================================================================================
635  } // if xytypecurrent ==xytype
636 
637  //====================== start road finder ====================================================
638  if (verbos > 0) {
639  std::cout << " start road finder " << std::endl;
640  }
641  do {
642  double fyY[20], fzY[20], fwY[20];
643  double fyYW[20], fwYW[20];
644  int py = 0, pys1 = 0, pyss = 0;
645  bool NewStation = false, py1first = false;
646  for (int sector = 1; sector < sn0; ++sector) {
647  double tav = 0., t1 = 0., t2 = 0., t = 0., sm;
648  int stattimes = 0;
649  if (sector != 1) {
650  NewStation = true;
651  }
652  for (int zmodule = 1; zmodule < pn0; ++zmodule) {
653  for (int justlayer = zbeg; justlayer < zmax; justlayer++) {
654  // iu is a continues numbering of planes(!)
655  int detfixed =
656  1; // use this treatment for each set up arm, hence no sense to do it differently for +FP420 and -FP420;
657  int nlayers = 3; // 2=3-1
658  unsigned int ii =
659  FP420NumberingScheme::packMYIndex(nlayers, pn0, sn0, detfixed, justlayer, sector, zmodule) -
660  1; // substruct 1 from 1(+1), 2(+2), 3(+3),4(+4),5...,6...,7...,8...,9...,10... (1st Station) ,11...,12...,13,...20... (2nd Station)
661  // ii = 0-19 --> 20 items
662 
663  if (nA[ii] != 0 && uA[ii] != nA[ii]) {
664  ++py;
665  if (sector == 1)
666  ++pys1;
667  if (sector == (sn0 - 1))
668  ++pyss;
669  if (py == 2 && sector == 1) {
670  // find closest cluster in X .
671  double dymin = 9999999., df2;
672  int cl2 = -1;
673  for (int cl = 0; cl < nA[ii]; ++cl) {
674  if (qA[cl][ii]) {
675  df2 = std::abs(fyY[fip] - yA[cl][ii]);
676  if (df2 < dymin) {
677  dymin = df2;
678  cl2 = cl;
679  } //if(df2
680  } //if(qA
681  } //for(cl
682  // END of finding of closest cluster in X .
683  if (cl2 != -1) {
684  t = (yA[cl2][ii] - fyY[fip]) / (zA[cl2][ii] - fzY[fip]);
685  t1 = t * wA[cl2][ii];
686  t2 = wA[cl2][ii];
687  if (verbos > 0) {
688  std::cout << " t= " << t << " tg0= " << tg0 << std::endl;
689  }
690  if (std::abs(t) < tg0) {
691  qA[cl2][ii] = false; //point is taken, mark it for not using again
692  fyY[py - 1] = yA[cl2][ii];
693  fzY[py - 1] = zA[cl2][ii];
694  fwY[py - 1] = wA[cl2][ii];
695  qAcl[py - 1] = cl2;
696  qAii[py - 1] = ii;
697  ++uA[ii];
698  if (verbos > 0) {
699  std::cout << " point is taken, mark it for not using again uA[ii]= " << uA[ii] << std::endl;
700  }
701  if (uA[ii] == nA[ii]) { /* no points anymore for this plane */
702  ++ry;
703  if (sector == 1)
704  ++rys1;
705  if (sector == (sn0 - 1))
706  ++ryss;
707  } //if(uA
708  } //if abs
709  else {
710  py--;
711  if (sector == 1)
712  pys1--;
713  if (sector == (sn0 - 1))
714  pyss--;
715  t1 -= t * wA[cl2][ii];
716  t2 -= wA[cl2][ii];
717  } //if(abs
718  } //if(cl2!=-1
719  else {
720  py--;
721  if (sector == 1)
722  pys1--;
723  if (sector == (sn0 - 1))
724  pyss--;
725  } //if(cl2!=-1
726  } //if(py==2
727  else {
728  // .
729  bool clLoopTrue = true;
730  int clcurr = -1;
731  for (int clind = 0; clind < nA[ii]; ++clind) {
732  if (clLoopTrue) {
733  int cl = clind;
734  if (qA[cl][ii]) {
735  clcurr = cl;
736  if (py < 3) {
737  if (py == 1) {
738  py1first = true;
739  fip = py - 1;
740  qA[cl][ii] = false; //point is taken, mark it for not using again
741  fyY[py - 1] = yA[cl][ii];
742  fzY[py - 1] = zA[cl][ii];
743  fwY[py - 1] = wA[cl][ii];
744  qAcl[py - 1] = cl;
745  qAii[py - 1] = ii;
746  ++uA[ii];
747  if (verbos > 0)
748  std::cout << " point is taken, mark it uA[ii]= " << uA[ii] << std::endl;
749  } //if py=1
750  if (uA[ii] == nA[ii]) { /* no points anymore for this plane */
751  ++ry;
752  if (sector == 1)
753  ++rys1;
754  if (sector == (sn0 - 1))
755  ++ryss;
756  } //if(uA
757  } //py<3
758  else {
759  if (NewStation) {
760  if (sn0 < 4) {
761  // stattimes=0 case (point of 1st plane to be matched in new Station)
762  sigma = ssigma;
763  } else {
764  sigma = ssigma / (sn0 - 1 - sector);
765  }
766  //sigma = ssigma/(sn0-sector);
767  //if(stattimes==1 || sector==3 ) sigma = msigma * sqrt(1./wA[cl][ii]);
768  if (stattimes == 1 || sector == 3)
769  sigma = sigmam; // (1st $3rd Stations for 3St. configur. ), 1st only for 2St. conf.
770  // if(stattimes==1 || sector==(sn0-1) ) sigma = sigmam;
771 
772  double cov00, cov01, cov11, c0Y, c1Y, chisqY;
773  gsl_fit_wlinear(fzY, 1, fwY, 1, fyY, 1, py - 1, &c0Y, &c1Y, &cov00, &cov01, &cov11, &chisqY);
774 
775  // find closest cluster in X .
776  int cl2match = -1;
777  double dymin = 9999999., df2;
778  for (int clmatch = 0; clmatch < nA[ii]; ++clmatch) {
779  if (qA[clmatch][ii]) {
780  double smmatch = c0Y + c1Y * zA[clmatch][ii];
781  df2 = std::abs(smmatch - yA[clmatch][ii]);
782  if (df2 < dymin) {
783  dymin = df2;
784  cl2match = clmatch;
785  } //if(df2
786  } //if(qA
787  } //for(clmatch
788 
789  if (cl2match != -1) {
790  cl = cl2match;
791  clLoopTrue = false; // just not continue the clinid loop
792  }
793 
794  sm = c0Y + c1Y * zA[cl][ii];
795 
796  if (verbos > 0) {
797  std::cout << " sector= " << sector << " sn0= " << sn0 << " sigma= " << sigma << std::endl;
798  std::cout << " stattimes= " << stattimes << " ssigma= " << ssigma << " sigmam= " << sigmam
799  << std::endl;
800  std::cout << " sm= " << sm << " c0Y= " << c0Y << " c1Y= " << c1Y << " chisqY= " << chisqY
801  << std::endl;
802  std::cout << " zA[cl][ii]= " << zA[cl][ii] << " ii= " << ii << " cl= " << cl << std::endl;
803  for (int ct = 0; ct < py - 1; ++ct) {
804  std::cout << " py-1= " << py - 1 << " fzY[ct]= " << fzY[ct] << std::endl;
805  std::cout << " fyY[ct]= " << fyY[ct] << " fwY[ct]= " << fwY[ct] << std::endl;
806  }
807  }
808 
809  } //NewStation 1
810  else {
811  t = (yA[cl][ii] - fyY[fip]) / (zA[cl][ii] - fzY[fip]);
812  t1 += t * wA[cl][ii];
813  t2 += wA[cl][ii];
814  tav = t1 / t2;
815  sm = fyY[fip] + tav * (zA[cl][ii] - fzY[fip]);
816  //sigma = nsigma * sqrt(1./wA[cl][ii]);
817  sigma = sigman;
818  }
819 
820  double diffpo = yA[cl][ii] - sm;
821  if (verbos > 0) {
822  std::cout << " diffpo= " << diffpo << " yA[cl][ii]= " << yA[cl][ii] << " sm= " << sm
823  << " sigma= " << sigma << std::endl;
824  }
825 
826  if (std::abs(diffpo) < sigma) {
827  if (NewStation) {
828  ++stattimes;
829  if (stattimes == 1) {
830  fip = py - 1;
831  t1 = 0;
832  t2 = 0;
833  } else if (stattimes == 2) {
834  NewStation = false;
835  t = (yA[cl][ii] - fyY[fip]) / (zA[cl][ii] - fzY[fip]);
836  //t1 += t*wA[cl][ii];
837  //t2 += wA[cl][ii];
838  t1 = t * wA[cl][ii];
839  t2 = wA[cl][ii];
840  } //if(stattime
841  } //if(NewStation 2
842  fyY[py - 1] = yA[cl][ii];
843  fzY[py - 1] = zA[cl][ii];
844  fwY[py - 1] = wA[cl][ii];
845  qA[cl][ii] = false; //point is taken, mark it for not using again
846  qAcl[py - 1] = cl;
847  qAii[py - 1] = ii;
848  ++uA[ii];
849  if (verbos > 0) {
850  std::cout << " 3333 point is taken, mark it uA[ii]= " << uA[ii] << std::endl;
851  }
852  if (uA[ii] == nA[ii]) { /* no points anymore for this plane */
853  ++ry;
854  if (sector == 1)
855  ++rys1;
856  if (sector == (sn0 - 1))
857  ++ryss;
858  } //if(cl==
859  // break; // to go on next plane
860  } //if abs
861  else {
862  t1 -= t * wA[cl][ii];
863  t2 -= wA[cl][ii];
864  } //if abs
865  } // if py<3 and else py>3
866 
867  if (!qA[cl][ii])
868  break; // go on next plane if point is found among clusters of current plane;
869  } // if qA
870  } // if clLoopTrue
871  } // for cl -- can be break and return to "for zmodule"
872  // .
873  if ((py != 1 && clcurr != -1 && qA[clcurr][ii]) || (py == 1 && !py1first)) {
874  // if point is not found - continue natural loop, but reduce py
875  py--;
876  if (sector == 1)
877  pys1--;
878  if (sector == (sn0 - 1))
879  pyss--;
880  } //if(py!=1
881  } //if(py==2 else
882  } //if(nA !=0 : inside this if( - ask ++py
883  } // for justlayer
884  } // for zmodule
885  } // for sector
886  //============
888 
889  if (verbos > 0) {
890  std::cout << "END: pys1= " << pys1 << " pyss = " << pyss << " py = " << py << std::endl;
891  }
892  // apply criteria for track selection:
893  // do not take track if
894  if (pys1 < pys1Cut || pyss < pyssCut || py < pyallCut) {
895  // if( pys1<3 || pyss<2 || py<4 ){
896  }
897  // do fit:
898  else {
900  double cov00, cov01, cov11;
901  double c0Y, c1Y, chisqY;
902  gsl_fit_wlinear(fzY, 1, fwY, 1, fyY, 1, py, &c0Y, &c1Y, &cov00, &cov01, &cov11, &chisqY);
904  // collect cases where the nearby points with the same coordinate exists
905  // int pyrepete=py+2;
906  // if(py < 11 && chisqY/(py-2) < 0.5) {
907  // double fyYold=999999.;
908  // for (int ipy=0; ipy<py; ++ipy) {
909  // if( fyY[ipy]!=fyYold) --pyrepete;
910  // fyYold=fyY[ipy];
911  // }
912  // }
914  float chindfx;
915  if (py > 2) {
916  chindfx = chisqY / (py - 2);
917  } else {
918  // chindfy = chisqY;
919  chindfx = 9999;
920  } //py
921  if (verbos > 0) {
922  // std::cout << " Do FIT XZ: chindfx= " << chindfx << " chisqY= " << chisqY << " py= " << py << " pyrepete= " << pyrepete << std::endl;
923  std::cout << " Do FIT XZ: chindfx= " << chindfx << " chisqY= " << chisqY << " py= " << py << std::endl;
924  }
925 
927  if (verbos > 0) {
928  std::cout << " preparation for second order fit for Wide pixels= " << std::endl;
929  }
930  for (int ipy = 0; ipy < py; ++ipy) {
931  if (xytypecurrent == 1) {
932  fyYW[ipy] = xYW[qAcl[ipy]][qAii[ipy]];
933  fwYW[ipy] = wYW[qAcl[ipy]][qAii[ipy]];
934  if (verbos > 0) {
935  std::cout << " ipy= " << ipy << std::endl;
936  std::cout << " qAcl[ipy]= " << qAcl[ipy] << " qAii[ipy]= " << qAii[ipy] << std::endl;
937  std::cout << " fyYW[ipy]= " << fyYW[ipy] << " fwYW[ipy]= " << fwYW[ipy] << std::endl;
938  }
939  } else if (xytypecurrent == 2) {
940  fyYW[ipy] = yXW[qAcl[ipy]][qAii[ipy]];
941  fwYW[ipy] = wXW[qAcl[ipy]][qAii[ipy]];
942  if (verbos == -29) {
943  std::cout << " ipy= " << ipy << std::endl;
944  std::cout << " qAcl[ipy]= " << qAcl[ipy] << " qAii[ipy]= " << qAii[ipy] << std::endl;
945  std::cout << " fyYW[ipy]= " << fyYW[ipy] << " fwYW[ipy]= " << fwYW[ipy] << std::endl;
946  }
947  }
948  } // for
949 
950  if (verbos > 0) {
951  std::cout << " start second order fit for Wide pixels= " << std::endl;
952  }
953  double wov00, wov01, wov11;
954  double w0Y, w1Y, whisqY;
955  gsl_fit_wlinear(fzY, 1, fwYW, 1, fyYW, 1, py, &w0Y, &w1Y, &wov00, &wov01, &wov11, &whisqY);
957 
958  float chindfy;
959  if (py > 2) {
960  chindfy = whisqY / (py - 2);
961  } else {
962  // chindfy = chisqY;
963  chindfy = 9999;
964  } //py
965 
966  if (verbos > 0) {
967  std::cout << " chindfy= " << chindfy << " chiCutY= " << chiCutY << std::endl;
968  }
969 
970  if (xytypecurrent == 1) {
971  if (chindfx < chiCutX && chindfy < chiCutY) {
972  ++numberYtracks;
973  Ay[numberYtracks - 1] = c0Y;
974  By[numberYtracks - 1] = c1Y;
975  Cy[numberYtracks - 1] = chisqY;
976  // My[numberYtracks-1] = py-pyrepete;
977  My[numberYtracks - 1] = py;
978  AyW[numberYtracks - 1] = w0Y;
979  ByW[numberYtracks - 1] = w1Y;
980  CyW[numberYtracks - 1] = whisqY;
981  MyW[numberYtracks - 1] = py;
982  if (verbos > 0) {
983  if (py > 20) {
984  std::cout << " niteration = " << niteration << std::endl;
985  std::cout << " chindfy= " << chindfy << " py= " << py << std::endl;
986  std::cout << " c0Y= " << c0Y << " c1Y= " << c1Y << std::endl;
987  std::cout << " pys1= " << pys1 << " pyss = " << pyss << std::endl;
988  }
989  }
990  } //chindfy
991  } else if (xytypecurrent == 2) {
992  if (chindfx < chiCutX && chindfy < chiCutY) {
993  ++numberXtracks;
994  Ax[numberXtracks - 1] = c0Y;
995  Bx[numberXtracks - 1] = c1Y;
996  Cx[numberXtracks - 1] = chisqY;
997  // Mx[numberXtracks-1] = py-pyrepete;
998  Mx[numberXtracks - 1] = py;
999  AxW[numberXtracks - 1] = w0Y;
1000  BxW[numberXtracks - 1] = w1Y;
1001  CxW[numberXtracks - 1] = whisqY;
1002  MxW[numberXtracks - 1] = py;
1003  if (verbos > 0) {
1004  std::cout << " niteration = " << niteration << std::endl;
1005  std::cout << " chindfx= " << chindfy << " px= " << py << std::endl;
1006  std::cout << " c0X= " << c0Y << " c1X= " << c1Y << std::endl;
1007  std::cout << " pxs1= " << pys1 << " pxss = " << pyss << std::endl;
1008  }
1009  } //chindfy
1010  }
1011 
1012  } // if else
1014  // do not select tracks anymore if
1015  if (verbos > 0) {
1016  std::cout << "Current iteration, niteration >= " << niteration << std::endl;
1017  std::cout << " numberYtracks= " << numberYtracks << std::endl;
1018  std::cout << " numberXtracks= " << numberXtracks << std::endl;
1019  std::cout << " pys1= " << pys1 << " pyss = " << pyss << " py = " << py << std::endl;
1020  std::cout << " tas1= " << tas1 << " tass = " << tass << " taf = " << taf << std::endl;
1021  std::cout << " rys1= " << rys1 << " ryss = " << ryss << " ry = " << ry << std::endl;
1022  std::cout << " tas1-rys1= " << tas1 - rys1 << " tass-ryss = " << tass - ryss << " taf-ry = " << taf - ry
1023  << std::endl;
1024  std::cout << "---------------------------------------------------------- " << std::endl;
1025  }
1026  // let's decide: do we continue track finder procedure
1027  if (tas1 - rys1 < pys1Cut || tass - ryss < pyssCut || taf - ry < pyallCut) {
1028  SelectTracks = false;
1029  } else {
1030  ++niteration;
1031  }
1032 
1033  } while (SelectTracks && niteration < nitMax);
1034  //====================== finish do loop finder for xytypecurrent ====================================================
1035 
1036  //============
1037 
1038  //===========================================================================================================================
1039 
1040  //===========================================================================================================================
1041  } // for xytypecurrent
1042  //===========================================================================================================================
1043 
1044  if (verbos > 0) {
1045  std::cout << " numberXtracks= " << numberXtracks << " numberYtracks= " << numberYtracks << std::endl;
1046  }
1047  //===========================================================================================================================
1048  //===========================================================================================================================
1049  //===========================================================================================================================
1050 
1051  // case X and Y plane types are available
1052  if (xytype > 2) {
1053  //===========================================================================================================================
1054  // match selected X and Y tracks to each other: tgphi=By/Bx->phi=artg(By/Bx); tgtheta=Bx/cosphi=By/sinphi-> ================
1055  // min of |Bx/cosphi-By/sinphi| ================
1056 
1057  //
1058  if (verbos > 0) {
1059  std::cout << " numberXtracks= " << numberXtracks << " numberYtracks= " << numberYtracks << std::endl;
1060  }
1061  if (numberXtracks > 0) {
1062  int newxnum[restracks], newynum[restracks]; // max # tracks = restracks = 10
1063  int nmathed = 0;
1064  do {
1065  double dthmin = 999999.;
1066  int trminx = -1, trminy = -1;
1067  for (int trx = 0; trx < numberXtracks; ++trx) {
1068  if (verbos > 0) {
1069  std::cout << "----------- trx= " << trx << " nmathed= " << nmathed << std::endl;
1070  }
1071  for (int tr = 0; tr < numberYtracks; ++tr) {
1072  if (verbos > 0) {
1073  std::cout << "--- tr= " << tr << " nmathed= " << nmathed << std::endl;
1074  }
1075  bool YesY = false;
1076  for (int nmx = 0; nmx < nmathed; ++nmx) {
1077  if (trx == newxnum[nmx])
1078  YesY = true;
1079  if (YesY)
1080  break;
1081  for (int nm = 0; nm < nmathed; ++nm) {
1082  if (tr == newynum[nm])
1083  YesY = true;
1084  if (YesY)
1085  break;
1086  }
1087  }
1088  if (!YesY) {
1089  //-------------------------------------------------------------------- ---- ---- ---- ---- ---- ----
1090  //double yyyyyy = 999999.;
1091  //if(Bx[trx] != 0.) yyyyyy = Ay[tr]-(Ax[trx]-xxxvtx)*By[tr]/Bx[trx];
1092  //double xxxxxx = 999999.;
1093  //if(By[tr] != 0.) xxxxxx = Ax[trx]-(Ay[tr]-yyyvtx)*Bx[trx]/By[tr];
1094  //double dthdif= std::abs(yyyyyy-yyyvtx) + std::abs(xxxxxx-xxxvtx);
1095 
1096  double dthdif = std::abs(AxW[trx] - Ay[tr]) + std::abs(BxW[trx] - By[tr]);
1097 
1098  if (verbos > 0) {
1099  // std::cout << " yyyyyy= " << yyyyyy << " xxxxxx= " << xxxxxx << " dthdif= " << dthdif << std::endl;
1100  std::cout << " abs(AxW[trx]-Ay[tr]) = " << std::abs(AxW[trx] - Ay[tr])
1101  << " abs(BxW[trx]-By[tr])= " << std::abs(BxW[trx] - By[tr]) << " dthdif= " << dthdif
1102  << std::endl;
1103  }
1104  //-------------------------------------------------------------------- ---- ---- ---- ---- ---- ----
1105  if (dthdif < dthmin) {
1106  dthmin = dthdif;
1107  trminx = trx;
1108  trminy = tr;
1109  } //if dthdif
1110  //--------------------------------------------------------------------
1111  } //if !YesY
1112  } //for y
1113  } // for x
1114  ++nmathed;
1115  if (trminx != -1) {
1116  newxnum[nmathed - 1] = trminx;
1117  } else {
1118  newxnum[nmathed - 1] = nmathed - 1;
1119  }
1120  if (verbos > 0) {
1121  std::cout << " trminx= " << trminx << std::endl;
1122  }
1123  if (nmathed > numberYtracks) {
1124  newynum[nmathed - 1] = -1;
1125  if (verbos > 0) {
1126  std::cout << "!!! nmathed= " << nmathed << " > numberYtracks= " << numberYtracks << std::endl;
1127  }
1128  } else {
1129  if (verbos > 0) {
1130  std::cout << " trminy= " << trminy << std::endl;
1131  }
1132  newynum[nmathed - 1] = trminy;
1133  }
1134  } while (nmathed < numberXtracks && nmathed < restracks);
1135 
1136  //
1137  //===========================================================================================================================
1138  //
1139  for (int tr = 0; tr < nmathed; ++tr) {
1140  int tx = newxnum[tr];
1141  int ty = newynum[tr];
1142  if (ty == -1) {
1143  ty = tx;
1144  Ay[ty] = 999.;
1145  By[ty] = 999.;
1146  Cy[ty] = 999.;
1147  My[ty] = -1;
1148  } //if ty
1149  // test:
1150  // tx=tr;
1151  //ty=tr;
1152  if (verbos > 0) {
1153  if (Mx[tx] > 20) {
1154  std::cout << " for track tr= " << tr << " tx= " << tx << " ty= " << ty << std::endl;
1155  std::cout << " Ax= " << Ax[tx] << " Ay= " << Ay[ty] << std::endl;
1156  std::cout << " Bx= " << Bx[tx] << " By= " << By[ty] << std::endl;
1157  std::cout << " Cx= " << Cx[tx] << " Cy= " << Cy[ty] << std::endl;
1158  std::cout << " Mx= " << Mx[tx] << " My= " << My[ty] << std::endl;
1159  std::cout << " AxW= " << AxW[tx] << " AyW= " << AyW[ty] << std::endl;
1160  std::cout << " BxW= " << BxW[tx] << " ByW= " << ByW[ty] << std::endl;
1161  std::cout << " CxW= " << CxW[tx] << " CyW= " << CyW[ty] << std::endl;
1162  std::cout << " MxW= " << MxW[tx] << " MyW= " << MyW[ty] << std::endl;
1163  }
1164  }
1165  // rhits.push_back( TrackFP420(c0X,c1X,chisqX,nhitplanesY,c0Y,c1Y,chisqY,nhitplanesY) );
1166  rhits.push_back(TrackFP420(Ax[tx], Bx[tx], Cx[tx], Mx[tx], Ay[ty], By[ty], Cy[ty], My[ty]));
1167  } //for tr
1168  //============================================================================================================
1169  } //in numberXtracks >0
1170  //============
1171 
1172  }
1173  // case Y plane types are available only
1174  else if (xytype == 1) {
1175  for (int ty = 0; ty < numberYtracks; ++ty) {
1176  if (verbos > 0) {
1177  std::cout << " for track ty= " << ty << std::endl;
1178  std::cout << " Ay= " << Ay[ty] << std::endl;
1179  std::cout << " By= " << By[ty] << std::endl;
1180  std::cout << " Cy= " << Cy[ty] << std::endl;
1181  std::cout << " My= " << My[ty] << std::endl;
1182  std::cout << " AyW= " << AyW[ty] << std::endl;
1183  std::cout << " ByW= " << ByW[ty] << std::endl;
1184  std::cout << " CyW= " << CyW[ty] << std::endl;
1185  std::cout << " MyW= " << MyW[ty] << std::endl;
1186  }
1187  rhits.push_back(TrackFP420(AyW[ty], ByW[ty], CyW[ty], MyW[ty], Ay[ty], By[ty], Cy[ty], My[ty]));
1188  } //for ty
1189  //============
1190  }
1191  // case X plane types are available only
1192  else if (xytype == 2) {
1193  for (int tx = 0; tx < numberXtracks; ++tx) {
1194  if (verbos > 0) {
1195  std::cout << " for track tx= " << tx << std::endl;
1196  std::cout << " Ax= " << Ax[tx] << std::endl;
1197  std::cout << " Bx= " << Bx[tx] << std::endl;
1198  std::cout << " Cx= " << Cx[tx] << std::endl;
1199  std::cout << " Mx= " << Mx[tx] << std::endl;
1200  std::cout << " AxW= " << AxW[tx] << std::endl;
1201  std::cout << " BxW= " << BxW[tx] << std::endl;
1202  std::cout << " CxW= " << CxW[tx] << std::endl;
1203  std::cout << " MxW= " << MxW[tx] << std::endl;
1204  }
1205  rhits.push_back(TrackFP420(Ax[tx], Bx[tx], Cx[tx], Mx[tx], AxW[tx], BxW[tx], CxW[tx], MxW[tx]));
1206  } //for tx
1207  //============
1208  } //xytype
1209 
1211 
1212  return rhits;
1213  //============
1214 }
FP420NumberingScheme::unpackCopyIndex
static int unpackCopyIndex(int rn0, int zside)
Definition: FP420NumberingScheme.h:69
ClusterFP420::barycenter
float barycenter() const
Definition: ClusterFP420.h:30
GetRecoTauVFromDQM_MC_cff.cl2
cl2
Definition: GetRecoTauVFromDQM_MC_cff.py:44
ClusterFP420::barycenterW
float barycenterW() const
Definition: ClusterFP420.h:33
ClusterCollectionFP420::ContainerIterator
std::vector< ClusterFP420 >::const_iterator ContainerIterator
Definition: ClusterCollectionFP420.h:11
SiStripMonitorCluster_cfi.zmax
zmax
Definition: SiStripMonitorCluster_cfi.py:215
FP420NumberingScheme::realzside
static int realzside(int rn0, int zsideinorder)
Definition: FP420NumberingScheme.h:99
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
FP420NumberingScheme::unpackLayerIndex
static int unpackLayerIndex(int rn0, int zside)
Definition: FP420NumberingScheme.h:56
input
static std::string const input
Definition: EdmProvDump.cc:47
ClusterFP420::barycerror
float barycerror() const
Definition: ClusterFP420.h:31
FP420NumberingScheme::packMYIndex
static unsigned packMYIndex(int rn0, int pn0, int sn0, int det, int zside, int sector, int zmodule)
Definition: FP420NumberingScheme.h:37
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
cuy.ii
ii
Definition: cuy.py:589
TrackProducerFP420::rhits
std::vector< TrackFP420 > rhits
Definition: TrackProducerFP420.h:66
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
gather_cfg.cout
cout
Definition: gather_cfg.py:144
ClusterCollectionFP420::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: ClusterCollectionFP420.h:12
ClusterFP420::amplitudes
const std::vector< short > & amplitudes() const
Definition: ClusterFP420.h:28
FP420NumberingScheme::unpackOrientation
static int unpackOrientation(int rn0, int zside)
Definition: FP420NumberingScheme.h:88
ClusterFP420::barycerrorW
float barycerrorW() const
Definition: ClusterFP420.h:34

Member Data Documentation

◆ chiCutX

float TrackProducerFP420::chiCutX
private

Definition at line 101 of file TrackProducerFP420.h.

◆ chiCutY

float TrackProducerFP420::chiCutY
private

Definition at line 102 of file TrackProducerFP420.h.

◆ dXX

double TrackProducerFP420::dXX
private

Definition at line 99 of file TrackProducerFP420.h.

◆ dYY

double TrackProducerFP420::dYY
private

Definition at line 100 of file TrackProducerFP420.h.

◆ gapBlade

double TrackProducerFP420::gapBlade
private

Definition at line 97 of file TrackProducerFP420.h.

◆ pitchX

double TrackProducerFP420::pitchX
private

Definition at line 87 of file TrackProducerFP420.h.

◆ pitchXW

double TrackProducerFP420::pitchXW
private

Definition at line 89 of file TrackProducerFP420.h.

◆ pitchY

double TrackProducerFP420::pitchY
private

Definition at line 88 of file TrackProducerFP420.h.

◆ pitchYW

double TrackProducerFP420::pitchYW
private

Definition at line 90 of file TrackProducerFP420.h.

◆ pn0

int TrackProducerFP420::pn0
private

Definition at line 71 of file TrackProducerFP420.h.

◆ rhits

std::vector<TrackFP420> TrackProducerFP420::rhits
private

Definition at line 66 of file TrackProducerFP420.h.

◆ rn0

int TrackProducerFP420::rn0
private

Definition at line 73 of file TrackProducerFP420.h.

◆ sn0

int TrackProducerFP420::sn0
private

Definition at line 69 of file TrackProducerFP420.h.

◆ soutput

ClusterCollectionFP420 TrackProducerFP420::soutput
private

Definition at line 64 of file TrackProducerFP420.h.

◆ UseHalfPitchShiftInX

bool TrackProducerFP420::UseHalfPitchShiftInX
private

Definition at line 78 of file TrackProducerFP420.h.

◆ UseHalfPitchShiftInXW

bool TrackProducerFP420::UseHalfPitchShiftInXW
private

Definition at line 80 of file TrackProducerFP420.h.

◆ UseHalfPitchShiftInY

bool TrackProducerFP420::UseHalfPitchShiftInY
private

Definition at line 79 of file TrackProducerFP420.h.

◆ UseHalfPitchShiftInYW

bool TrackProducerFP420::UseHalfPitchShiftInYW
private

Definition at line 81 of file TrackProducerFP420.h.

◆ verbos

int TrackProducerFP420::verbos
private

Definition at line 106 of file TrackProducerFP420.h.

◆ XsensorSize

double TrackProducerFP420::XsensorSize
private

Definition at line 108 of file TrackProducerFP420.h.

◆ xytype

int TrackProducerFP420::xytype
private

Definition at line 75 of file TrackProducerFP420.h.

◆ YsensorSize

double TrackProducerFP420::YsensorSize
private

Definition at line 109 of file TrackProducerFP420.h.

◆ z420

double TrackProducerFP420::z420
private

Definition at line 84 of file TrackProducerFP420.h.

◆ zBlade

double TrackProducerFP420::zBlade
private

Definition at line 96 of file TrackProducerFP420.h.

◆ zD2

double TrackProducerFP420::zD2
private

Definition at line 85 of file TrackProducerFP420.h.

◆ zD3

double TrackProducerFP420::zD3
private

Definition at line 86 of file TrackProducerFP420.h.

◆ ZGapLDet

double TrackProducerFP420::ZGapLDet
private

Definition at line 91 of file TrackProducerFP420.h.

◆ zinibeg

double TrackProducerFP420::zinibeg
private

Definition at line 104 of file TrackProducerFP420.h.

◆ ZSiDet

double TrackProducerFP420::ZSiDet
private

Definition at line 95 of file TrackProducerFP420.h.

◆ ZSiPlane

double TrackProducerFP420::ZSiPlane
private

Definition at line 94 of file TrackProducerFP420.h.

◆ ZSiStep

double TrackProducerFP420::ZSiStep
private

Definition at line 93 of file TrackProducerFP420.h.

Generated for CMSSW Reference Manual by doxygen 1.8.14