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