CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Types | Public Member Functions | Private Attributes
TrackProducerFP420 Class Reference

#include <TrackProducerFP420.h>

Public Types

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

Public Member Functions

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

Private Attributes

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

Detailed Description

Definition at line 25 of file TrackProducerFP420.h.

Member Typedef Documentation

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

Definition at line 28 of file TrackProducerFP420.h.

Constructor & Destructor Documentation

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

Definition at line 23 of file TrackProducerFP420.cc.

References gather_cfg::cout, and PFRecoTauDiscriminationAgainstMuon2_cfi::verbosity.

23  {
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 }
PFRecoTauDiscriminationAgainstMuon2_cfi.verbosity
tuple verbosity
Definition: PFRecoTauDiscriminationAgainstMuon2_cfi.py:26
TrackProducerFP420::theFP420NumberingScheme
FP420NumberingScheme * theFP420NumberingScheme
Definition: TrackProducerFP420.h:48
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121

Member Function Documentation

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

Definition at line 86 of file TrackProducerFP420.cc.

References abs, ClusterFP420::amplitudes(), ClusterFP420::barycenter(), ClusterFP420::barycenterW(), ClusterFP420::barycerror(), ClusterFP420::barycerrorW(), cl, GetRecoTauVFromDQM_MC_cff::cl2, gather_cfg::cout, corrVsCorr::ry, and lumiQTWidget::t.

Referenced by FP420TrackMain::run().

86  {
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 }
ClusterFP420::barycerror
float barycerror() const
Definition: ClusterFP420.h:32
ClusterCollectionFP420::ContainerIterator
std::vector< ClusterFP420 >::const_iterator ContainerIterator
Definition: ClusterCollectionFP420.h:13
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
GetRecoTauVFromDQM_MC_cff.cl2
tuple cl2
Definition: GetRecoTauVFromDQM_MC_cff.py:43
TrackProducerFP420::theFP420NumberingScheme
FP420NumberingScheme * theFP420NumberingScheme
Definition: TrackProducerFP420.h:48
cl
float cl
Definition: Combine.cc:71
TrackProducerFP420::rhits
std::vector< TrackFP420 > rhits
Definition: TrackProducerFP420.h:46
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
ClusterCollectionFP420::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: ClusterCollectionFP420.h:14

Member Data Documentation

float TrackProducerFP420::chiCutX
private

Definition at line 84 of file TrackProducerFP420.h.

float TrackProducerFP420::chiCutY
private

Definition at line 85 of file TrackProducerFP420.h.

double TrackProducerFP420::dXX
private

Definition at line 82 of file TrackProducerFP420.h.

double TrackProducerFP420::dYY
private

Definition at line 83 of file TrackProducerFP420.h.

double TrackProducerFP420::gapBlade
private

Definition at line 80 of file TrackProducerFP420.h.

double TrackProducerFP420::pitchX
private

Definition at line 70 of file TrackProducerFP420.h.

double TrackProducerFP420::pitchXW
private

Definition at line 72 of file TrackProducerFP420.h.

double TrackProducerFP420::pitchY
private

Definition at line 71 of file TrackProducerFP420.h.

double TrackProducerFP420::pitchYW
private

Definition at line 73 of file TrackProducerFP420.h.

int TrackProducerFP420::pn0
private

Definition at line 54 of file TrackProducerFP420.h.

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

Definition at line 46 of file TrackProducerFP420.h.

int TrackProducerFP420::rn0
private

Definition at line 56 of file TrackProducerFP420.h.

int TrackProducerFP420::sn0
private

Definition at line 52 of file TrackProducerFP420.h.

ClusterCollectionFP420 TrackProducerFP420::soutput
private

Definition at line 44 of file TrackProducerFP420.h.

FP420NumberingScheme* TrackProducerFP420::theFP420NumberingScheme
private

Definition at line 48 of file TrackProducerFP420.h.

bool TrackProducerFP420::UseHalfPitchShiftInX
private

Definition at line 61 of file TrackProducerFP420.h.

bool TrackProducerFP420::UseHalfPitchShiftInXW
private

Definition at line 63 of file TrackProducerFP420.h.

bool TrackProducerFP420::UseHalfPitchShiftInY
private

Definition at line 62 of file TrackProducerFP420.h.

bool TrackProducerFP420::UseHalfPitchShiftInYW
private

Definition at line 64 of file TrackProducerFP420.h.

int TrackProducerFP420::verbos
private

Definition at line 89 of file TrackProducerFP420.h.

double TrackProducerFP420::XsensorSize
private

Definition at line 91 of file TrackProducerFP420.h.

int TrackProducerFP420::xytype
private

Definition at line 58 of file TrackProducerFP420.h.

double TrackProducerFP420::YsensorSize
private

Definition at line 92 of file TrackProducerFP420.h.

double TrackProducerFP420::z420
private

Definition at line 67 of file TrackProducerFP420.h.

double TrackProducerFP420::zBlade
private

Definition at line 79 of file TrackProducerFP420.h.

double TrackProducerFP420::zD2
private

Definition at line 68 of file TrackProducerFP420.h.

double TrackProducerFP420::zD3
private

Definition at line 69 of file TrackProducerFP420.h.

double TrackProducerFP420::ZGapLDet
private

Definition at line 74 of file TrackProducerFP420.h.

double TrackProducerFP420::zinibeg
private

Definition at line 87 of file TrackProducerFP420.h.

double TrackProducerFP420::ZSiDet
private

Definition at line 78 of file TrackProducerFP420.h.

double TrackProducerFP420::ZSiPlane
private

Definition at line 77 of file TrackProducerFP420.h.

double TrackProducerFP420::ZSiStep
private

Definition at line 76 of file TrackProducerFP420.h.

Generated for CMSSW Reference Manual by doxygen 1.8.5