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