StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
gl3LMVertexFinder.cxx
1 #include "gl3LMVertexFinder.h"
2 
3 #include "gl3Event.h"
4 #include "gl3HistoManager.h"
5 
6 #include <assert.h>
7 #include <rtsLog.h>
8 
9 //#define GL3ROOT
10 #ifdef GL3ROOT
11 static int lmvdbg=1;
12 #else
13 static int lmvdbg=0;
14 #endif
15 const float gl3LMVertexFinder::C_PI=3.1416926;
16 
17 //####################################################################
18 // Constructor
19 //####################################################################
20 gl3LMVertexFinder::gl3LMVertexFinder(int _minHitsOnTrack,
21  int _minCTBadc,
22  float _maxZdca,
23  float _zMargin,
24  float _phiMargin,
25  float _delVertMax)
26 {
27  setParameters(_minHitsOnTrack, _minCTBadc, _maxZdca,
28  _zMargin, _phiMargin, _delVertMax);
29 
30  init();
31  registerHistos();
32 }
33 
34 
35 //####################################################################
36 // Destructor
37 //####################################################################
38 gl3LMVertexFinder::~gl3LMVertexFinder()
39 {
40 }
41 
42 //####################################################################
43 //
44 //####################################################################
45 int gl3LMVertexFinder::setParameters(int _minHitsOnTrack,
46  int _minCTBadc,
47  float _maxZdca,
48  float _zMargin,
49  float _phiMargin,
50  float _delVertMax)
51 {
52  CtbRadius=213.; //(cm)
53  CtbLen2=255.; //(cm) for eta=[0,1]
54 
55  minHitsOnTrack = _minHitsOnTrack;
56  minCTBadc = _minCTBadc;
57 
58  maxZdca = _maxZdca;
59  MatchCtbMaxZ = CtbLen2/4. + _zMargin; //(cm) =slat+margin/cm
60  MatchCtbMaxPhi = C_PI/30+_phiMargin/CtbRadius;//slat half width+margin/cm
61  DelVertMax = _delVertMax; //(cm) rejection of outliers in ppLMV
62 
63  return 1;
64 }
65 
66 
67 
68 //####################################################################
69 //
70 //####################################################################
71 int gl3LMVertexFinder::setParameters(const char *filename)
72 {
73  int minHitsOnTrack;
74  int minCTBadc;
75  float maxZdca;
76  float zMargin;
77  float phiMargin;
78  float delVertMax;
79 
80  // set defaults, if reading form file fails ...
81  minHitsOnTrack = 10;
82  minCTBadc = 3; // set to -1 if no CTB data avaliable, without CTB works only for low luminosity and bXing ID is not sure any more
83 
84  maxZdca = 250.0;
85  zMargin = 10.0;
86  phiMargin = 3.0;
87  delVertMax = 3.0;
88 
89  // open jan's para txt file
90  FILE *fp = fopen(filename,"r");
91  if(fp == NULL) return -1;
92 
93  int ret=-1;
94  if(fp) {
95  ret=fscanf(fp,"%d %d %f %f %f %f",
96  &minHitsOnTrack, &minCTBadc, &maxZdca,
97  &zMargin, &phiMargin, &delVertMax);
98  fclose(fp);
99  }
100  if(ret!=6) LOG(ERR,"gl3LMVertexFinder::setParameters() ERROR while reading %s file, ret=%d, default params used\n",filename,ret,0,0,0);
101 
102  if(minCTBadc<0) LOG(ERR, "gl3LMVertexFinder::setParameters() WARN : CTB matching is off\n",0,0,0,0,0);
103 
104  setParameters(minHitsOnTrack, minCTBadc, maxZdca,
105  zMargin, phiMargin, delVertMax);
106 
107  return 0;
108 }
109 
110 //####################################################################
111 //
112 //####################################################################
113 int gl3LMVertexFinder::init () {
114 
115  totInpEve=0;
116  setCtbMap();
117 
118  if(lmvdbg) LOG(ERR, "CTB matching limits del_Phi(deg)=%.1f, delZ(cm)=%f\n",
119  MatchCtbMaxPhi/C_PI*180, MatchCtbMaxZ,0,0,0);
120 
121  //
122  // Book histos
123  //
124  hjan[0]= new gl3Histo ( "ppLMV0","ZDC vertex in Z (cm)",
125  100, -200, 200. ) ;
126  hjan[1]= new gl3Histo ( "ppLMV1","raw event track multiplicity",
127  100, 0.,500.);
128  hjan[2]= new gl3Histo ( "ppLMV2"," event multiplicity,with any CTB slat",
129  50, 0., 500. );
130  hjan[3]= new gl3Histo ( "ppLMV3"," HITS IN TRACK, all,with any CTB slat ",
131  50, -0.5, 49.5 );
132  hjan[4]= new gl3Histo ( "ppLMV4"," Zdca (cm), long tracks",
133  50, -300., 300. );
134  hjan[5]= new gl3Histo ( "ppLMV5"," Zdca-ZvertexZDC (cm)",
135  50, -200., 200. );
136  hjan[6]= new gl3Histo ( "ppLMV6","non-zero ADC's for CTB",
137  51, -0.5, 50.5 );
138  hjan[7]= new gl3Histo ( "ppLMV7"," CTB slats with ADC>th",
139  240, -0.5, 239.5 );
140  hjan[8]= new gl3Histo ( "ppLMV8"," <z> of CTB with ADC>th",
141  10, -300., 300. );
142  hjan[9]= new gl3Histo ( "ppLMV9","pasX #Delta#phi (deg) from MatchCTB",
143  100, -10., 10. );
144  hjan[10]= new gl3Histo("ppLMV10","pasX #DeltaZ (cm) from MatchCTB",
145  50, -100., 100. );
146  hjan[11]= new gl3Histo("ppLMV11","number of tracks matched to CTB/eve",
147  50, -0.5, 49.5 );
148  hjan[12]= new gl3Histo("ppLMV12","number of accepted CTB slats/eve",
149  250, -0.5, 249.5 );
150  hjan[13]= new gl3Histo( "ppLMV13","counter of accepted events",
151  7,0.5,7.5);
152  hjan[14]= new gl3Histo( "ppLMV14","counter of accepted tracks",
153  9,0.5,9.5);
154  hjan[15]= new gl3Histo ( "ppLMV15","my vertex in Z (cm)",
155  50, -300, 300. ) ;
156  hjan[16]= new gl3Histo ( "ppLMV16","my-ZDC vertex in Z (cm)",
157  50, -20, 20. ) ;
158  hjan[17]= new gl3Histo("ppLMV17","end: #DeltaZ (cm) (matched track-myVertex)",
159  50, -5., 5. );
160  hjan[18]= new gl3Histo("ppLMV18","#DeltaZ (cm) (matched track-myVertex), Nmatch=[2-5]",
161  50, -5., 5. );
162  hjan[19]= new gl3Histo("ppLMV19","#DeltaZ (cm) (matched track-myVertex), Nmatch=[6-10]",
163  50, -5., 5. );
164  hjan[20]= new gl3Histo("ppLMV20","#DeltaZ (cm) (matched track-myVertex), Nmatch=[10++]",
165  50, -5., 5. );
166  hjan[21]= new gl3Histo("ppLMV21","sigma of tracks extrapolated to vertex, a.u.",
167  50,.0,10.);
168 
169  nhjan=22;
170 
171  if(lmvdbg) {
172  for(int i=0;i<16;i++) {
173  char t1[100], t2[100];
174  sprintf(t1,"ppeve%d",i);
175  sprintf(t2,"Zdca-Zvertex (cm) for accepted eve=%2d",i);
176  heve[i]= new gl3Histo ( t1,t2, 50, -50., 50. );
177  }
178  }
179 
180  neve=0;
181 
182  return 0 ;
183 }
184 
185 
186 //####################################################################
187 //
188 //####################################################################
189 int gl3LMVertexFinder::registerHistos() {
190 
191  int j;
192  for(j=0;j<nhjan;j++) gl3HistoManager::instance()->add(hjan[j]);
193 
194  if(lmvdbg) for(j=0;j<16;j++) gl3HistoManager::instance()->add(heve[j]);
195 
196  return 0;
197 }
198 
199 
200 //####################################################################
201 //
202 //####################################################################
203 int gl3LMVertexFinder::makeVertex( gl3Event *event ) {
204  int j;
205 
206  totInpEve++;
207  myVert.z=2001;
208  myVert.x=2002;
209  myVert.y=2003;
210 
211  NtrackL=0;
212 
213  memset(trackL,0,sizeof(trackL));
214  hjan[13]->Fill(1);
215  hjan[1]->Fill(event->getNTracks(),1);
216 
217  if (event->getNTracks() > 5000)
218  return 0;
219 
220 
221  if( getCtbHits(event) <=0) return 0; // no CTB slats found
222 
223  hjan[13]->Fill(2);
224  // if( NctbH> maxCtbSlats) return 0; // too many CTB slats found
225 
226  hjan[13]->Fill(3);
227 
228  float zVertexZDC=event->getZDCVertex();
229 
230  hjan[0]->Fill(zVertexZDC);
231  hjan[2]->Fill(event->getNTracks(),1);
232 
233  //
234  // loop over gtracks
235  //
236 
237  for ( int trkcnt = 0 ; trkcnt < event->getNTracks() ; trkcnt++ ) {
238  hjan[14]->Fill(1);
239  gl3Track* gTrack = event->getTrack(trkcnt);
240 
241  gTrack->flag=0; // reset track flag to NULL
242 
243  hjan[3]->Fill(gTrack->nHits);
244  if ( gTrack->nHits < minHitsOnTrack ) continue ;
245 
246  hjan[14]->Fill(2);
247 
248  // Extrapolate track to beam-DCA
249  Ftf3DHit dcaHit1=gTrack->closestApproach ( 0., 0. ) ;
250  float zDCA = dcaHit1.z ;
251 
252  //gTrack->updateToClosestApproach ( 0., 0. ) ;
253  //float zDCA = gTrack->z0 ;
254 
255  if(fabs(zDCA)>maxZdca) continue;
256 
257  hjan[14]->Fill(3);
258 
259  if(minCTBadc>=0) { // do match to CTB
260  Ftf3DHit hit=gTrack->extraRadius(CtbRadius);
261  if(fabs(hit.x)<0.1 && fabs(hit.y)<0.1) continue; // extrapolation to CTB failed
262 
263  hjan[14]->Fill(4);
264  if( matchTrack2Ctb(&hit) <0) continue; // no match
265 
266  hjan[14]->Fill(5);
267  } else { // no matching to CTB, take all tracks
268  hjan[14]->Fill(6);
269  }
270 
271  if(lmvdbg) printf("track DCA x=%f y=%f z=%f\n",dcaHit1.x, dcaHit1.y,dcaHit1.z );
272  // gTrack->Print(30);
273 
274  if(NtrackL>= MAXTRACK) {
275  LOG(ERR, " PileupFilter WARN: %d=NtrackL>= MAXTRACK, skip this matched track\n",NtrackL,0,0,0,0);
276  continue;
277  }
278  trackL[NtrackL]=gTrack;
279  XdcaHitL[NtrackL]=dcaHit1.x;
280  YdcaHitL[NtrackL]=dcaHit1.y;
281  ZdcaHitL[NtrackL]=dcaHit1.z;
282  NtrackL++;
283 
284  hjan[4]->Fill(zDCA);
285  hjan[5]->Fill(zDCA-zVertexZDC);
286  }// end of loop over tracks
287 
288  neve++;
289  hjan[11]->Fill(NtrackL);
290 
291  if(lmvdbg) printf(" - end of eve: accepted tracks=%d, ZDCvert=%f.1 (neve=%d)\n",NtrackL,zVertexZDC,neve);
292 
293  // add beam line at X=Y=0 as as high momentum track
294  gl3Track beamTrack;
295  beamTrack.r0=0;
296  beamTrack.phi0=0;
297  beamTrack.z0=0;
298  beamTrack.tanl=888999.; // almost parallel to Z
299  beamTrack.psi=0;
300  beamTrack.id=888999;
301  beamTrack.lastHit=0;
302 
303  trackL[NtrackL]=&beamTrack;
304  XdcaHitL[NtrackL]=0.;// should not be used
305  YdcaHitL[NtrackL]=0.;// should not be used
306  ZdcaHitL[NtrackL]=0.;// should not be used
307  NtrackL++;
308 
309  if(NtrackL<2) return 0;
310 
311  hjan[13]->Fill(4);
312 
313  // vertex finder
314  int NtrackVertex=ppLMV4b(&myVert);
315  if( NtrackVertex <2) return 0; // no vertex found
316 
317 
318  // output:
319  // NtrackMatch;
320  // myVert.z;
321 
322  hjan[13]->Fill(5);
323  hjan[15]->Fill(myVert.z);
324  hjan[16]->Fill(myVert.z-zVertexZDC);
325 
326 
327  for( j=0;j<NtrackL;j++) {
328  if( trackL[j]==NULL) continue;
329  if(trackL[j]->lastHit==NULL) continue; // ignore 'beamTrack'
330  float zVer = ZdcaHitL[j] ;
331  hjan[17]->Fill(zVer-myVert.z);
332  if(lmvdbg && neve<16) heve[neve]->Fill(zVer-myVert.z);
333  }
334 
335 
336  for( j=0;j<NtrackL;j++) {
337  if( trackL[j]==NULL) continue;
338  if(trackL[j]->lastHit==NULL) continue; // ignore 'beamTrack'
339  float zDca = ZdcaHitL[j] ;
340  hjan[17]->Fill(zDca-myVert.z);
341  if( NtrackVertex <=5 ) hjan[18]->Fill(zDca-myVert.z);
342  else if( NtrackVertex <=10 ) hjan[19]->Fill(zDca-myVert.z);
343  else hjan[20]->Fill(zDca-myVert.z);
344  }
345 
346 
347  if (NtrackVertex<2)
348  return 0; // failed to find vertex
349  else
350  return 1; // Found a vertex
351 }
352 
353 //qTotal += gTrack->q ;
354 //pxTotal += gTrack->pt * cos(gTrack->psi);
355 //pyTotal += gTrack->pt * sin(gTrack->psi);
356 //pzTotal += gTrack->pt * gTrack->tanl ;
357 
358 
359 //####################################################################
360 //
361 //####################################################################
362 int gl3LMVertexFinder::getCtbHits (gl3Event *event) {
363 
364  memset(ctbH,0,sizeof(ctbH));
365  NctbH=0; // clear old response
366 
367  int nSlat=0;
368 
369  float phiZero1 = 72 ; // magic lines from Pablo & Herb
370  float phiZero2 = 108 ;
371  float deltaPhi = 6 ;
372 
373  for(int ss=0;ss<2;ss++) {
374  for(int tt=0;tt<120;tt++) {
375  int indx=ctbMap[tt][ss];
376  if(indx<0||indx>255) {
377  //printf ( "gl3LMVertexFinder::getCtbHits ( ) Ctb index %d out of bounds, skip hit\n", indx ) ;
378  continue ;
379  }
380  int ctbAdc =event->getCTB(indx);
381  //printf("ss=%d, tt=%d, indx=%d Adc=%d\n",ss,tt,indx,ctbAdc);
382 
383  if(ctbAdc)hjan[6]->Fill(ctbAdc);
384  //assert(ctbAdc==0);
385  if(ctbAdc< minCTBadc) continue;
386 
387  int iz ;
388  float phi ;
389 
390  if ( tt < 60 ) {
391  phi = phiZero1 - tt * deltaPhi ;
392  iz = 0 ;
393  } else {
394  phi = phiZero2 + (tt-60) * deltaPhi ;
395  iz = 1 ;
396  }
397  if ( phi < 0. ) phi += 360 ;
398  if ( phi > 360. ) phi -= 360 ;
399 
400  //printf ( "ctb-hit tt=%d ss=%d adc=%d phi=%f z=%f idx=%d\n",
401  // printf ( "%d %d %d %f %f %d\n",
402 // tt+1, ss+1, ctbAdc, phi/180*C_PI,
403 // zSlats[iz][ss], indx) ;
404 
405  if(nSlat>= MAXSLATS) {
406  if(lmvdbg)
407  printf(" PileupFilter WARN: %d=nSlat>= MAXSLATS,nSlat, skip this CTB hit\n",nSlat);
408  continue;
409  }
410 
411  ctbH[nSlat].z= zSlats[iz][ss];
412  ctbH[nSlat].phi=phi/180*C_PI;
413  ctbH[nSlat].adc=ctbAdc;
414  ctbH[nSlat].slatIndex=indx;
415  // printf("slat Indx: %i ADC: %i z=%.1f phi/deg=%.1f\n",indx,ctbAdc,ctbH[nSlat].z,ctbH[nSlat].phi/C_PI*180);
416  hjan[7]->Fill(tt+120*ss);
417  hjan[8]->Fill(ctbH[nSlat].z);
418  nSlat++;
419  } // end of loop over tt
420  } // end of loop over ss
421  if(lmvdbg) printf(" %d CTB slats with ADC>=%d\n",nSlat,minCTBadc); //tmp
422  NctbH=nSlat;
423  hjan[12]->Fill(NctbH);
424  return NctbH;
425 }
426 
427 //####################################################################
428 //
429 //####################################################################
430 int gl3LMVertexFinder::matchTrack2Ctb( Ftf3DHit * h) {
431  assert(h);
432 
433  float phi=atan2(h->y,h->x);
434  if(phi<0) phi+=2*C_PI;
435  // printf("track at CTB x=%f, y=%f, z=%f, phi/deg=%.1f\n",h->x,h->y,h->z,phi/3.1416*180);
436  for(int is=0;is<NctbH;is++) {
437 
438  float del_z=h->z-ctbH[is].z;
439  if(fabs(del_z) >MatchCtbMaxZ) continue;
440 
441  float del_phi=phi-ctbH[is].phi;
442  if(del_phi>C_PI) del_phi-=2*C_PI;
443  if(del_phi<-C_PI) del_phi+=2*C_PI;
444 
445  // printf("Match is=%d DZ=%f, Dphi=%f %f\n",is,del_z,del_phi,MatchCtbMaxPhi);
446 
447  if(fabs(del_phi) >MatchCtbMaxPhi) continue;
448 
449  hjan[9]->Fill(del_phi/C_PI*180.);
450  hjan[10]->Fill(del_z);
451  //printf("Match OK ctbHit_is=%d DZ=%f, Dphi=%f\n",is,del_z,del_phi);
452  return is;
453  }
454 
455  return -1 ; // no match
456 }
457 
458 //####################################################################
459 //
460 //####################################################################
461 int gl3LMVertexFinder::ppLMV4b ( Ftf3DHit *myV) {// vertex finder
462 
463  myV->x=1000;
464  myV->y=1001;
465  myV->z=1002;
466  int NtrackUsed=NtrackL;
467 
468  //Do the actual vertex fitting, continue until good
469  double A11=0.0,A12=0.0,A13=0.0,A21=0.0,A22=0.0,A23=0.0;
470  double A31=0.0,A32=0.0,A33=0.0; // Matrix Elements
471  double C11=0.0,C12=0.0,C13=0.0,C21=0.0,C22=0.0,C23=0.0;
472  double C31=0.0,C32=0.0,C33=0.0; // C = A^-1
473  int done = 0;
474  double chi2=0;
475 
476  while( done != 1 ){
477 
478  // Check that there at least are 2 tracks
479  if( NtrackUsed <= 1 ){
480  if(lmvdbg)cout<<"ppLMV4b: Fewer than 2 track remains. No vertex found."<<endl;
481  return -1;
482  }
483 
484  // Begin by doing a fit
485  A11=0.0,A12=0.0,A13=0.0,A21=0.0,A22=0.0,A23=0.0;
486  A31=0.0,A32=0.0,A33=0.0; // Matrix Elements
487  C11=0.0,C12=0.0,C13=0.0,C21=0.0,C22=0.0,C23=0.0;
488  C31=0.0,C32=0.0,C33=0.0; // C = A^-1
489  double b1=0.0,b2=0.0,b3=0.0;
490  // Compute matrix A and vector b
491  int itr;
492  for(itr=0; itr < NtrackL; itr++){
493 
494  gl3Track *trk=trackL[itr];
495  if(trk==NULL) continue;
496  double xp=XdcaHitL[itr];
497  double yp=YdcaHitL[itr];
498  double zp=ZdcaHitL[itr];
499 
500  double mag=sqrt(1+trk->tanl*trk->tanl);
501  double xhat=cos(trk->psi)/mag; // wrong, temp not at DCA
502  double yhat=sin(trk->psi)/mag; // wrong, temp not at DCA
503  double zhat=trk->tanl/mag;
504 
505  /* the transverse extrapolation error of the track DR= Dtheta*L
506  where Dtheta is angular straggling, L is path from vertex to closest point on the track
507  Below following approximations are made:
508  - all particles are pions
509  - Dtheta=const /beta/P
510  - P= PT**sqrt(1+tanl**2)
511  - L = Rxy(hit)*sqrt(1+tanl**2)
512  so DR=const * Rxy /beta/PT
513  */
514 
515  double sigma=1;
516  gl3Hit *h=(gl3Hit *)trk->lastHit;
517  if(h) {
518  float Rxy=sqrt(h->getX()*h->getX() + h->getY()*h->getY());
519  float pmom =trk->pt*mag;
520  float beta= pmom/sqrt(pmom*pmom+ 0.139*0.139);
521  sigma=0.0136 *Rxy/beta/trk->pt; // average sigma is 2.2
522  } else {
523  if(trk->id==888999 && trk->tanl>88899.) sigma=0.5; // 'beamTrack'
524  }
525 
526  hjan[21]->Fill(sigma);
527 
528  A11=A11+(yhat*yhat+zhat*zhat)/sigma;
529  A12=A12-(xhat*yhat)/sigma;
530  A13=A13-(xhat*zhat)/sigma;
531  A22=A22+(xhat*xhat+zhat*zhat)/sigma;
532  A23=A23-(yhat*zhat)/sigma;
533  A33=A33+(xhat*xhat+yhat*yhat)/sigma;
534  b1=b1 + ( (yhat*yhat+zhat*zhat)*xp - xhat*yhat*yp - xhat*zhat*zp )/sigma;
535  b2=b2 + ( (xhat*xhat+zhat*zhat)*yp - xhat*yhat*xp - yhat*zhat*zp )/sigma;
536  b3=b3 + ( (xhat*xhat+yhat*yhat)*zp - xhat*zhat*xp - yhat*zhat*yp )/sigma;
537  }// end of 1-st loop over tracks
538 
539  A21 = A12; A31=A13; A32=A23;
540 
541  // Invert A
542  double detA = A11*A22*A33 + A12*A23*A31 + A13*A21*A32;
543  detA = detA - A31*A22*A13 - A32*A23*A11 - A33*A21*A12;
544  // cout<<"Determinant= "<<detA<<endl;
545  // cout<<"A11,A12,A13: "<<A11<<" "<<A12<<" "<<A13<<endl;
546  // cout<<"A21,A22,A23: "<<A21<<" "<<A22<<" "<<A23<<endl;
547  // cout<<"A31,A32,A33: "<<A31<<" "<<A32<<" "<<A33<<endl;
548  // cout<<"b1,b2,b3 "<<b1<<" "<<b2<<" "<<b3<<endl;
549  C11=(A22*A33-A23*A32)/detA; C12=(A13*A32-A12*A33)/detA; C13=(A12*A23-A13*A22)/detA;
550  C21=C12; C22=(A11*A33-A13*A31)/detA; C23=(A13*A21-A11*A23)/detA;
551  C31=C13; C32=C23; C33=(A11*A22-A12*A21)/detA;
552 
553  // Find Vertex Position
554  double Xv = C11*b1 + C12*b2 + C13*b3;
555  double Yv = C21*b1 + C22*b2 + C23*b3;
556  double Zv = C31*b1 + C32*b2 + C33*b3;
557 
558  if(lmvdbg) {
559  cout<<"current Vertex Position : "<<Xv<<" "<<Yv<<" "<<Zv<<" nTR used="<<NtrackUsed<<endl;
560  cout<<" Vertex Error : "<<sqrt(C11)<<" "<<sqrt(C22)<<" "<<sqrt(C33)<<endl;
561  }
562 
563  // Check if the fit is any good
564  // Loop over tracks again to get Chi2 and check each track's deviation
565  double dmax=0.0;
566  int iworse=-1;
567 
568  for(itr=0; itr < NtrackL; itr++){
569  gl3Track *trk=trackL[itr];
570  if(trk==NULL) continue;
571  if(trk->lastHit==NULL) continue; // ignore 'beamTrack'
572 
573  double xp=XdcaHitL[itr];
574  double yp=YdcaHitL[itr];
575  double zp=ZdcaHitL[itr];
576 
577  if(lmvdbg) printf("itr=%d, Zdca=%f Zvert=%f\n",itr,zp,Zv);
578  double d2=(xp-Xv)*(xp-Xv) +(yp-Yv)*(yp-Yv) +(zp-Zv)*(zp-Zv);
579  double d=sqrt(d2);
580  double sig=1; // temp2
581  chi2 = chi2 + d2/(sig*sig);
582  double drel = d/sig;
583  if(lmvdbg) printf(" itr=%d, drel/sig=%f d=%f/sig sig=%f\n",itr,drel,d, sig);
584  if( drel > dmax ){
585  // Save the track that deviates the most from vertex
586  dmax = drel;
587  iworse = itr;
588  }
589  }
590 
591  if( dmax > DelVertMax ){
592  if(lmvdbg) cout<<"Removing a track "<<iworse<< " with dmax= "<<dmax<<endl;
593  trackL[iworse]=NULL;
594  NtrackUsed--;
595  done=0;
596  }
597  else{
598  done=1;
599  myV->x=Xv;
600  myV->y=Yv;
601  myV->z=Zv;
602  }
603  } // End While Loop
604 
605  //double chi2pdof = chi2/(NtrackUsed-1);
606  if(lmvdbg) {
607  cout << "ppLMV4b: Primary Vertex found: used tracks="
608  << NtrackUsed << " out of " << NtrackL << " matched and "
609  << /*event->getNTracks() << */" L3-tracks" << endl
610  <<" at Position : " << myV->x << " " << myV->y << " "
611  << myV->z << endl;
612  }
613  //tag used tracks
614  for(int itr=0; itr < NtrackL; itr++){
615  gl3Track *trk=trackL[itr];
616  if(trk==NULL) continue;
617  trk->flag+=0x2; // mark track as used by lmv-l3
618  }
619 
620  return NtrackUsed;
621 }
622 
623 
624 
625 //####################################################################
626 //
627 //####################################################################
628 void gl3LMVertexFinder::setCtbMap() {
629  ctbMap[1-1][1-1]=109;
630  ctbMap[2-1][1-1]=108;
631  ctbMap[3-1][1-1]=107;
632  ctbMap[4-1][1-1]=106;
633  ctbMap[5-1][1-1]=105;
634  ctbMap[6-1][1-1]=7;
635  ctbMap[7-1][1-1]=6;
636  ctbMap[8-1][1-1]=5;
637  ctbMap[9-1][1-1]=4;
638  ctbMap[10-1][1-1]=3;
639  ctbMap[11-1][1-1]=2;
640  ctbMap[12-1][1-1]=1;
641  ctbMap[13-1][1-1]=0;
642  ctbMap[14-1][1-1]=15;
643  ctbMap[15-1][1-1]=14;
644  ctbMap[16-1][1-1]=13;
645  ctbMap[17-1][1-1]=12;
646  ctbMap[18-1][1-1]=11;
647  ctbMap[19-1][1-1]=10;
648  ctbMap[20-1][1-1]=9;
649  ctbMap[21-1][1-1]=39;
650  ctbMap[22-1][1-1]=38;
651  ctbMap[23-1][1-1]=37;
652  ctbMap[24-1][1-1]=36;
653  ctbMap[25-1][1-1]=35;
654  ctbMap[26-1][1-1]=34;
655  ctbMap[27-1][1-1]=33;
656  ctbMap[28-1][1-1]=32;
657  ctbMap[29-1][1-1]=47;
658  ctbMap[30-1][1-1]=46;
659  ctbMap[31-1][1-1]=45;
660  ctbMap[32-1][1-1]=44;
661  ctbMap[33-1][1-1]=43;
662  ctbMap[34-1][1-1]=42;
663  ctbMap[35-1][1-1]=41;
664  ctbMap[36-1][1-1]=71;
665  ctbMap[37-1][1-1]=70;
666  ctbMap[38-1][1-1]=69;
667  ctbMap[39-1][1-1]=68;
668  ctbMap[40-1][1-1]=67;
669  ctbMap[41-1][1-1]=66;
670  ctbMap[42-1][1-1]=65;
671  ctbMap[43-1][1-1]=64;
672  ctbMap[44-1][1-1]=79;
673  ctbMap[45-1][1-1]=78;
674  ctbMap[46-1][1-1]=77;
675  ctbMap[47-1][1-1]=76;
676  ctbMap[48-1][1-1]=75;
677  ctbMap[49-1][1-1]=74;
678  ctbMap[50-1][1-1]=73;
679  ctbMap[51-1][1-1]=103;
680  ctbMap[52-1][1-1]=102;
681  ctbMap[53-1][1-1]=101;
682  ctbMap[54-1][1-1]=100;
683  ctbMap[55-1][1-1]=99;
684  ctbMap[56-1][1-1]=98;
685  ctbMap[57-1][1-1]=97;
686  ctbMap[58-1][1-1]=96;
687  ctbMap[59-1][1-1]=111;
688  ctbMap[60-1][1-1]=110;
689  ctbMap[1-1][2-1]=125;
690  ctbMap[2-1][2-1]=124;
691  ctbMap[3-1][2-1]=123;
692  ctbMap[4-1][2-1]=122;
693  ctbMap[5-1][2-1]=121;
694  ctbMap[6-1][2-1]=23;
695  ctbMap[7-1][2-1]=22;
696  ctbMap[8-1][2-1]=21;
697  ctbMap[9-1][2-1]=20;
698  ctbMap[10-1][2-1]=19;
699  ctbMap[11-1][2-1]=18;
700  ctbMap[12-1][2-1]=17;
701  ctbMap[13-1][2-1]=16;
702  ctbMap[14-1][2-1]=31;
703  ctbMap[15-1][2-1]=30;
704  ctbMap[16-1][2-1]=29;
705  ctbMap[17-1][2-1]=28;
706  ctbMap[18-1][2-1]=27;
707  ctbMap[19-1][2-1]=26;
708  ctbMap[20-1][2-1]=25;
709  ctbMap[21-1][2-1]=55;
710  ctbMap[22-1][2-1]=54;
711  ctbMap[23-1][2-1]=53;
712  ctbMap[24-1][2-1]=52;
713  ctbMap[25-1][2-1]=51;
714  ctbMap[26-1][2-1]=50;
715  ctbMap[27-1][2-1]=49;
716  ctbMap[28-1][2-1]=48;
717  ctbMap[29-1][2-1]=63;
718  ctbMap[30-1][2-1]=62;
719  ctbMap[31-1][2-1]=61;
720  ctbMap[32-1][2-1]=60;
721  ctbMap[33-1][2-1]=59;
722  ctbMap[34-1][2-1]=58;
723  ctbMap[35-1][2-1]=57;
724  ctbMap[36-1][2-1]=87;
725  ctbMap[37-1][2-1]=86;
726  ctbMap[38-1][2-1]=85;
727  ctbMap[39-1][2-1]=84;
728  ctbMap[40-1][2-1]=83;
729  ctbMap[41-1][2-1]=82;
730  ctbMap[42-1][2-1]=81;
731  ctbMap[43-1][2-1]=80;
732  ctbMap[44-1][2-1]=95;
733  ctbMap[45-1][2-1]=94;
734  ctbMap[46-1][2-1]=93;
735  ctbMap[47-1][2-1]=92;
736  ctbMap[48-1][2-1]=91;
737  ctbMap[49-1][2-1]=90;
738  ctbMap[50-1][2-1]=89;
739  ctbMap[51-1][2-1]=119;
740  ctbMap[52-1][2-1]=118;
741  ctbMap[53-1][2-1]=117;
742  ctbMap[54-1][2-1]=116;
743  ctbMap[55-1][2-1]=115;
744  ctbMap[56-1][2-1]=114;
745  ctbMap[57-1][2-1]=113;
746  ctbMap[58-1][2-1]=112;
747  ctbMap[59-1][2-1]=127;
748  ctbMap[60-1][2-1]=126;
749  ctbMap[61-1][1-1]=141;
750  ctbMap[62-1][1-1]=140;
751  ctbMap[63-1][1-1]=139;
752  ctbMap[64-1][1-1]=138;
753  ctbMap[65-1][1-1]=137;
754  ctbMap[66-1][1-1]=167;
755  ctbMap[67-1][1-1]=166;
756  ctbMap[68-1][1-1]=165;
757  ctbMap[69-1][1-1]=164;
758  ctbMap[70-1][1-1]=163;
759  ctbMap[71-1][1-1]=162;
760  ctbMap[72-1][1-1]=161;
761  ctbMap[73-1][1-1]=160;
762  ctbMap[74-1][1-1]=175;
763  ctbMap[75-1][1-1]=174;
764  ctbMap[76-1][1-1]=173;
765  ctbMap[77-1][1-1]=172;
766  ctbMap[78-1][1-1]=171;
767  ctbMap[79-1][1-1]=170;
768  ctbMap[80-1][1-1]=169;
769  ctbMap[81-1][1-1]=199;
770  ctbMap[82-1][1-1]=198;
771  ctbMap[83-1][1-1]=197;
772  ctbMap[84-1][1-1]=196;
773  ctbMap[85-1][1-1]=195;
774  ctbMap[86-1][1-1]=194;
775  ctbMap[87-1][1-1]=193;
776  ctbMap[88-1][1-1]=192;
777  ctbMap[89-1][1-1]=207;
778  ctbMap[90-1][1-1]=206;
779  ctbMap[91-1][1-1]=205;
780  ctbMap[92-1][1-1]=204;
781  ctbMap[93-1][1-1]=203;
782  ctbMap[94-1][1-1]=202;
783  ctbMap[95-1][1-1]=201;
784  ctbMap[96-1][1-1]=231;
785  ctbMap[97-1][1-1]=230;
786  ctbMap[98-1][1-1]=229;
787  ctbMap[99-1][1-1]=228;
788  ctbMap[100-1][1-1]=227;
789  ctbMap[101-1][1-1]=226;
790  ctbMap[102-1][1-1]=225;
791  ctbMap[103-1][1-1]=224;
792  ctbMap[104-1][1-1]=239;
793  ctbMap[105-1][1-1]=239;
794  ctbMap[106-1][1-1]=237;
795  ctbMap[107-1][1-1]=236;
796  ctbMap[108-1][1-1]=235;
797  ctbMap[109-1][1-1]=234;
798  ctbMap[110-1][1-1]=233;
799  ctbMap[111-1][1-1]=135;
800  ctbMap[112-1][1-1]=134;
801  ctbMap[113-1][1-1]=133;
802  ctbMap[114-1][1-1]=132;
803  ctbMap[115-1][1-1]=131;
804  ctbMap[116-1][1-1]=130;
805  ctbMap[117-1][1-1]=129;
806  ctbMap[118-1][1-1]=128;
807  ctbMap[119-1][1-1]=143;
808  ctbMap[120-1][1-1]=142;
809  ctbMap[61-1][2-1]=157;
810  ctbMap[62-1][2-1]=156;
811  ctbMap[63-1][2-1]=155;
812  ctbMap[64-1][2-1]=154;
813  ctbMap[65-1][2-1]=153;
814  ctbMap[66-1][2-1]=183;
815  ctbMap[67-1][2-1]=182;
816  ctbMap[68-1][2-1]=181;
817  ctbMap[69-1][2-1]=180;
818  ctbMap[70-1][2-1]=179;
819  ctbMap[71-1][2-1]=178;
820  ctbMap[72-1][2-1]=177;
821  ctbMap[73-1][2-1]=176;
822  ctbMap[74-1][2-1]=191;
823  ctbMap[75-1][2-1]=190;
824  ctbMap[76-1][2-1]=189;
825  ctbMap[77-1][2-1]=188;
826  ctbMap[78-1][2-1]=187;
827  ctbMap[79-1][2-1]=186;
828  ctbMap[80-1][2-1]=185;
829  ctbMap[81-1][2-1]=215;
830  ctbMap[82-1][2-1]=214;
831  ctbMap[83-1][2-1]=213;
832  ctbMap[84-1][2-1]=212;
833  ctbMap[85-1][2-1]=211;
834  ctbMap[86-1][2-1]=210;
835  ctbMap[87-1][2-1]=209;
836  ctbMap[88-1][2-1]=208;
837  ctbMap[89-1][2-1]=223;
838  ctbMap[90-1][2-1]=222;
839  ctbMap[91-1][2-1]=221;
840  ctbMap[92-1][2-1]=220;
841  ctbMap[93-1][2-1]=219;
842  ctbMap[94-1][2-1]=218;
843  ctbMap[95-1][2-1]=217;
844  ctbMap[96-1][2-1]=247;
845  ctbMap[97-1][2-1]=246;
846  ctbMap[98-1][2-1]=245;
847  ctbMap[99-1][2-1]=244;
848  ctbMap[100-1][2-1]=243;
849  ctbMap[101-1][2-1]=242;
850  ctbMap[102-1][2-1]=241;
851  ctbMap[103-1][2-1]=240;
852  ctbMap[104-1][2-1]=255;
853  ctbMap[105-1][2-1]=254;
854  ctbMap[106-1][2-1]=253;
855  ctbMap[107-1][2-1]=252;
856  ctbMap[108-1][2-1]=251;
857  ctbMap[109-1][2-1]=250;
858  ctbMap[110-1][2-1]=249;
859  ctbMap[111-1][2-1]=151;
860  ctbMap[112-1][2-1]=150;
861  ctbMap[113-1][2-1]=149;
862  ctbMap[114-1][2-1]=148;
863  ctbMap[115-1][2-1]=147;
864  ctbMap[116-1][2-1]=146;
865  ctbMap[117-1][2-1]=145;
866  ctbMap[118-1][2-1]=144;
867  ctbMap[119-1][2-1]=159;
868  ctbMap[120-1][2-1]=158;
869 
870  zSlats[0][0] = 64;;
871  zSlats[0][1] = 191. ;
872  zSlats[1][0] = -64.;
873  zSlats[1][1] = -191.;
874 
875 
876 #ifdef GL3ROOT
877  printf("the Ctb Map from Pablo from herb is set\n");
878 #endif
879  return ;
880 }
881 
882 
883 
gl3Hit
Definition: gl3Hit.h:24
gl3Histo::Fill
int Fill(double x, double weight=1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ...
Definition: gl3Histo.cxx:53
Generated by   doxygen 1.8.5