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