StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFgtStraightTrackMaker.cxx
1 #include "StFgtStraightTrackMaker.h"
3 #include "StRoot/StEvent/StFgtCollection.h"
4 #include "StRoot/StEvent/StFgtHitCollection.h"
5 #include "StRoot/StEvent/StFgtHit.h"
6 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
7 #include "StRoot/StEvent/StEvent.h"
8 #include "StRoot/StEvent/StEventInfo.h"
9 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
10 
11 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
12 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
13 #include "StMuDSTMaker/COMMON/StMuDst.h"
14 #include "StMuDSTMaker/COMMON/StMuEvent.h"
15 #include "StarClassLibrary/StThreeVectorF.hh"
16 
17 
18 
19 #include <TH2D.h>
20 #include <TROOT.h>
21 #include <TStyle.h>
22 #include <TCanvas.h>
23 #include <utility>
24 #include <TArc.h>
25 #include <TLine.h>
26 #include <set>
27 
28 //#define PRINT_1D
29 //max num clusters any disk is allowed to have
30 
31 //#define COSMIC
32 #include "StFgtCosmicAlignment.h"
33 #define CHARGE_MEASURE clusterCharge
34 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
35 //#define LEN_CONDITION
36 
37 
38 #include "StRoot/StEvent/StEvent.h"
39 #include "StRoot/StEvent/StFgtCollection.h"
40 //disk for which I want to calculate efficieny
41 
42 
43 #define MAX_CHARGE_RATIO
44 #define MIN_CHARGE_RATIO
45 
46 //#define DISK_EFF 2
47 #define MY_PI 3.14159265359
48 
49 #define DISK_DIM 40
50 #define NUM_EFF_BIN 100
51 
52 
53 pair<double,double> StFgtStraightTrackMaker::getDca( vector<AVTrack>::iterator it)
54 {
55 
56  float oldDist=100000;
57  float zStep=1.0;
58  float optZ=-200;
59  float dist=0;
60  for(float z=-100;z<100;z=z+zStep)
61  {
62  float x=it->mx*z+it->ax;
63  float y=it->my*z+it->ay;
64  dist=x*x+y*y;
65  if(dist>oldDist)
66  {
68  dist=oldDist;
69  optZ=z-zStep;//last z was better
70  break;
71  }
72  else
73  oldDist=dist;
74  }
75  return pair<double,double>(optZ,sqrt(dist));
76 }
77 
78 
79 
80 //if required only return clusters with matches...
81 pair<Double_t,Double_t> StFgtStraightTrackMaker::findCluChargeSize(Int_t iD,Char_t layer, Double_t ordinate)
82 {
83  if(iD<0 || iD >5)
84  {
85  return pair<Double_t,Double_t>(-99999,-99999);
86  }
87  vector<generalCluster> &hitVec=*(pClusters[iD]);
88  Double_t charge=-99999;
89  Double_t cluSize=-9999;
90  Double_t minDist=99999;
91  for(vector<generalCluster>::iterator it=hitVec.begin();it!=hitVec.end();it++)
92  {
93  if(it->layer!=layer)
94  continue;
95  Float_t ord=-9999;
96  if(layer=='R')
97  {
98  ord=it->posR;
99  }
100  else
101  {
102  ord=it->posPhi;
103  if(fabs(ordinate-(ord+MY_PI))<fabs(ord-ordinate))
104  ord=ord+MY_PI;
105  if(fabs(ordinate-(ord-MY_PI))<fabs(ord-ordinate))
106  ord=ord-MY_PI;
107  }
108 
109  if((fabs(ordinate-ord)<minDist) && (!useChargeMatch || it->hasMatch))
110  {
111  charge=it->clusterCharge;
112  cluSize=it->clusterSize;
113  minDist=fabs(ordinate-ord);
114  }
115  }
116  return pair<Double_t,Double_t>(charge,cluSize);
117 }
118 
119 Double_t StFgtStraightTrackMaker::findClosestPoint(float mx, float bx, float my, float by, double xE, double yE, Int_t iD)
120 {
121  // cout <<"expecting point at " << xE <<", " <<yE <<endl;
122  if(iD<0 || iD >5)
123  {
124  return 99999;
125  }
126  vector<generalCluster> &hitVec=*(pClusters[iD]);
127  Double_t dist2=99999;
128  for(vector<generalCluster>::iterator it=hitVec.begin();it!=hitVec.end();it++)
129  {
130  for(vector<generalCluster>::iterator it2=hitVec.begin();it2!=hitVec.end();it2++)
131  {
132  //
133  if(useChargeMatch && !StFgtGeneralBase::arePointsMatched(it,it2))
134  continue;
135 
136  if(it->layer==it2->layer)
137  continue;
138  Float_t r=it->posR;
139  Float_t phi=it2->posPhi;
140  if(it->layer!='R')
141  {
142  phi=it->posPhi;
143  r=it2->posR;
144  }
145 
146  Float_t x=r*cos(phi);
147  Float_t y=r*sin(phi);
148  // cout <<"we have " << x <<", " << y <<endl;
149  Double_t mDist=(x-xE)*(x-xE)+(y-yE)*(y-yE);
150  if(mDist<dist2)
151  {
152  dist2=mDist;
153  //recalculate distance with proper alignment
154  float tmpX, tmpY,tmpZ,tmpP,tmpR;
155  if(isCosmic)
156  {
157 
158  getAlign(iD,phi,r,tmpX,tmpY,tmpZ,tmpP,tmpR);
159  Double_t xExpUpdate=mx*tmpZ+bx;
160  Double_t yExpUpdate=my*tmpZ+by;
161  // cout<<"tmpx: " << tmpX <<" old: " << x <<" xE old: " << xE << " updated: " << xExpUpdate;
162  // cout<<"tmpy: " << tmpY <<" old: " << y <<" yE old: " << yE << " updated: " << yExpUpdate<<endl;
163  mDist=(tmpX-xExpUpdate)*(tmpX-xExpUpdate)+(tmpY-yExpUpdate)*(tmpY-yExpUpdate);
164  dist2=mDist;
165 
166  }
168  // Double_t yExp=my*StFgtGeom::getDiscZ(i)+by;
169 
170 
171  // (*outTxtFile) <<"point found, x: " << x <<" y: " << y << " dist: " << dist2 <<endl;
172  }
173  }
174  }
175  // (*outTxtFile) <<"returning : " << dist2<<endl;
176  return dist2;
177 }
178 
179 
181 Short_t StFgtStraightTrackMaker::getQuadFromCoo(Double_t x, Double_t y)
182 {
183  cout <<"do not use this function!!!" <<endl;
184  if(x>0 && y>0)
185  return 0;
186  if(x>0 && y<0)
187  return 1;
188  if(x<0 && y<0)
189  return 2;
190  if(x<0 && y>0)
191  return 3;
192 
193  return -9999;
194 }
195 
196 
197 //print the strips around the place where we expect hit
198 pair<Double_t,Double_t> StFgtStraightTrackMaker::getChargeRatio(Float_t r, Float_t phi, Int_t iD, Int_t iq)
199 {
200  //first r:
201  Double_t maxRCharge=-9999;
202  Double_t maxPhiCharge=-9999;
203  Int_t maxRInd=-1;
204  Int_t maxPInd=-1;
205  for(unsigned int i=0;i< pStrips[iD*4+iq].size();i++)
206  {
207  Int_t geoId=pStrips[iD*4+iq][i].geoId;
208  generalStrip& pStrip=pStrips[iD*4+iq][i];
209  Short_t disc, quadrant,strip;
210  Char_t layer;
211  Double_t ordinate, lowerSpan, upperSpan;//, prvOrdinate;
212  StFgtGeom::getPhysicalCoordinate(geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
213  StFgtGeom::decodeGeoId(geoId,disc, quadrant, layer, strip);
214  // if(layer=='P' && disc==iD && iq==quadrant)
215  // cout <<"looking for " << phi << " have: " << ordinate <<" diff: " << fabs(ordinate-phi) <<endl;
216  if(disc==iD && iq==quadrant && ((layer =='R' && fabs(ordinate-r)<0.7) || (layer=='P' && fabs(ordinate-phi)<0.03) || (layer=='P' && fabs(ordinate-phi+2*MY_PI)<0.03 ) || (layer=='P' && fabs(ordinate-phi-2*MY_PI)<0.03)|| (layer=='P' && fabs(ordinate-phi+MY_PI)<0.03 ) || (layer=='P' && fabs(ordinate-phi-MY_PI)<0.03)))
217  {
218  if(layer=='P')
219  {
220  if(pStrip.charge>maxPhiCharge)
221  {
222  maxPhiCharge=pStrip.charge;
223  maxPInd=i;
224  }
225  }
226  else
227  {
228  if(pStrip.charge>maxRCharge)
229  {
230  maxRCharge=pStrip.charge;
231  maxRInd=i;
232  }
233  }
234  }
235  }
236  if(maxRInd>=0 && maxPInd>=0)
237  {
238  if(maxRCharge>1000 && maxPhiCharge>1000)
239  {
241  if(maxRInd>0)
242  maxRCharge+=pStrips[iD*4+iq][maxRInd-1].charge;
243  if(maxRInd< (int)(pStrips[iD*4+iq].size()-1))
244  maxRCharge+=pStrips[iD*4+iq][maxRInd+1].charge;
245  if(maxPInd>0)
246  maxPhiCharge+=pStrips[iD*4+iq][maxPInd-1].charge;
247  if(maxPInd< (int)(pStrips[iD*4+iq].size()-1))
248  maxPhiCharge+=pStrips[iD*4+iq][maxPInd+1].charge;
249 
250  return pair<Double_t,Double_t>(maxRCharge,maxPhiCharge);
251  }
252 
253  }
254  return pair<Double_t,Double_t>(-9999,-9999);
255 }
256 
257 
258 
259 Bool_t StFgtStraightTrackMaker::getTrack(vector<AVPoint>& points, Double_t ipZ)
260 {
261  // (*outTxtFile) <<"getTrack" <<endl;
262  // cout <<"get track" <<endl;
263  ipZ=-9999; //get ourselves
264  StFgtGeneralBase *fgtGenMkr = static_cast<StFgtGeneralBase * >( GetMaker("fgtGenBase"));
265 
266  ipZ=fgtGenMkr->vtxZ;
267  Int_t vtxRank=fgtGenMkr->vtxRank;
268  vector<AVPoint>::iterator iter=points.begin();
269  Double_t A = 0, B = 0, Cx = 0, Cy = 0, D = 0, Ex = 0, Ey = 0;
270  Double_t dist = 0;
271 
272  //LOG_INFO << " --------------------------> calling makeLine with " << points.size() << " 0x" << std::hex << discBitArray << std::dec << endm;
273  // cout <<"get track2" <<endl;
274  //do the alignment
275  if(isCosmic)
276  {
277 
278  float tmpX, tmpY,tmpZ,tmpP,tmpR;
279  for( ; iter != points.end(); ++iter ){
280  getAlign(iter->dID,iter->phi,iter->r,tmpX,tmpY,tmpZ,tmpP,tmpR);
281  // cout <<"before: " << iter->phi << ", " << iter->r <<" " << iter->x <<" " << iter->y << ", " << iter->z <<endl;
282  iter->phi=tmpP;
283  iter->r=tmpR;
284  iter->x=tmpX;
285  iter->y=tmpY;
286  iter->z=tmpZ;
287  // cout <<"after: " << iter->phi << ", " << iter->r <<" " << iter->x <<" " << iter->y << ", " << iter->z <<endl;
288  }
289  }
290  else
291  {
292  // cout <<"no cosmic" <<endl;
293  }
294 
295  // cout <<" we have " << points.size() << " points " << endl;
296  for( iter=points.begin(); iter != points.end(); ++iter ){
297 
298  Double_t x = iter->x;
299  Double_t y = iter->y;
300  Double_t z = iter->z;
301  // cout <<"x: " << x << " y: " << y << " z: " << z <<endl;
302  A += z*z;
303  B += z;
304  Cx += x*z;
305  Cy += y*z;
306  Ex += x;
307  Ey += y;
308 
309  // cout << "*** Point located at " << x << ' ' << y << ' ' << z << " Disk: " << iter->dID <<endl;
310  }
311  // cout <<"ipZ: " << ipZ <<endl;
312  // cout <<"get track3" <<endl;
313  D = points.size();
314  //invalid is -999
315 
316  if(doFitWithVertex)
317  {
318 
319  if(isMuDst && vtxRank>0&&ipZ>-100 && ipZ<100)
320  {
321  cout <<"fit with vertex!!" <<endl;
322  A+=ipZ*ipZ;
323  B+=ipZ;
324  D++;
325  }
326  }
327  //cout << "*** Consts " << A << ' ' << B << ' ' << Cx << ' ' << Cy << ' ' << D << ' ' << Ex << ' ' << Ey << endl;
328  Double_t denom = D*A - B*B;
329  if( denom )
330  {
331  Double_t bx = (-B*Cx + A*Ex)/denom;
332  Double_t by = (-B*Cy + A*Ey)/denom;
333  Double_t mx = ( D*Cx - B*Ex)/denom;
334  Double_t my = ( D*Cy - B*Ey)/denom;
335  // cout <<"bx: " << bx <<" by: " << by <<" mx: " << mx << " my: " << my <<endl;
336  dist = 0;
337  // cout <<"get track4" <<endl;
338  vector<AVPoint> redPoints;
339  // find closest points...
340  for(int iDx=0;iDx<6;iDx++)
341  {
342  Double_t minDistance=9999;
343  Int_t pointIdx=-1;
344  Int_t cnt=0;
345  for( iter = points.begin(); iter != points.end(); ++iter ){
346  Int_t dId=iter->dID;
347  Double_t distX=fabs((iter->z*mx+bx)-(iter->x));
348  Double_t distY=fabs((iter->z*my+by)-(iter->y));
349  // cout <<"distX: " << distX <<" distY: " << distY <<endl;
350  Double_t distance=distX*distX+distY*distY;
351  // cout << " got distance " << distance << endl;
352  if((iDx==dId)&&(distance<minDistance))
353  {
354  minDistance=distance;
355  pointIdx=cnt;
356  // cout <<"min dist now:" << minDistance<<" for this dik " << iDx <<" pointIdx: " << pointIdx <<endl;
357  }
358  cnt++;
359  }
360  if(pointIdx>=0)
361  {
362  // cout <<"pushing back " << pointIdx <<endl;
363  redPoints.push_back(points[pointIdx]);
364  }
365  }//end of looping over discs
366 
369  // cout <<"doing refit... " <<endl;
370  // cout <<"get track5" <<endl;
371  A=0.0;
372  B=0.0;
373  Cx=0.0;
374  Cy=0.0;
375  Ex=0.0;
376  Ey=0.0;
377 
378  for(vector<AVPoint>::iterator iterR=redPoints.begin() ; iterR != redPoints.end(); ++iterR )
379  {
380  Double_t x = iterR->x;
381  Double_t y = iterR->y;
382  Double_t z = iterR->z;
383 
384  A += z*z;
385  B += z;
386  Cx += x*z;
387  Cy += y*z;
388  Ex += x;
389  Ey += y;
390  }
391  D = redPoints.size();
392  if(doRefitWithVertex)
393  {
394  if(isMuDst && vtxRank>0&&ipZ>-100 && ipZ<100)
395  {
396  cout <<"refit with vertex! " <<endl;
397  A+=ipZ*ipZ;
398  B+=ipZ;
399  D++;
400  }
401  }
402 
403  Double_t denom = D*A - B*B;
404  // cout <<"get track6" <<endl;
405  if( denom )
406  {
407  bx = (-B*Cx + A*Ex)/denom;
408  by = (-B*Cy + A*Ey)/denom;
409  mx = ( D*Cx - B*Ex)/denom;
410  my = ( D*Cy - B*Ey)/denom;
411  }
412  // cout <<"after refit: bx: " << bx <<" by: " << by <<" mx: " << mx << " my: " << my <<endl;
413  // cout <<"we have refit line: "<< bx << " by: " << by <<" mx: " << mx << " my: " << my <<endl;
415  for(vector<AVPoint>::iterator iterR = redPoints.begin(); iterR != redPoints.end(); ++iterR )
416  {
417  Double_t distX, distY;
418  distX=fabs((iterR->z*mx+bx)-(iterR->x));
419  distY=fabs((iterR->z*my+by)-(iterR->y));
420  // cout <<"distX: " << distX <<" distY: " << distY <<endl;
421  dist += (distX*distX+distY*distY);
422  // cout <<"adding " << (distX*distX+distY*distY) <<" to chi2: " << endl;
423  // cout << "*** DistSq " << dist << endl;
424  }
425  // cout <<"get track8" <<endl;
426  dist/=D;
427  // cout <<" end chi2: " <<dist << ", ip: " << ipZ <<endl;
428  //this ipZ is the TPC vertex! Is reset below...
429  m_tracks.push_back(AVTrack(mx,my,bx,by,ipZ,dist));
430  // cout <<" we have " <<m_tracks.size() <<" track now " <<endl;
431  (m_tracks.back()).vtxRank=vtxRank;
432  points.clear();
433  for(vector<AVPoint>::iterator it=redPoints.begin();it!=redPoints.end();it++)
434  {
435  points.push_back(*it);
436  }
437 
438  for(vector<AVPoint>::iterator iter = points.begin(); iter != points.end(); ++iter ){
439  // cout << "--- Location at each disc at z: " << iter->z << " "
440  // << "X: " << mx*iter->z+bx << " vs " << iter->x << ' '
441  // << "Y: " << my*iter->z+by << " vs " << iter->y << " "
442  // << " charge phi: " << iter->phiCharge <<" rcharge: "<< iter->rCharge <<endl;
443 
444  }
445  // cout <<endl<<endl;
446  // cout <<"dist again: " <<dist <<endl;
447  vector<AVTrack>::iterator it_lastTrack=m_tracks.end();
448  it_lastTrack--;
449  pair<double,double> dca=getDca(it_lastTrack);
450  Double_t vertZ = ( -( it_lastTrack->mx*it_lastTrack->ax + it_lastTrack->my*it_lastTrack->ay )/(it_lastTrack->mx*it_lastTrack->mx+it_lastTrack->my*it_lastTrack->my));
451  (it_lastTrack)->trkZ=vertZ;
452  // cout <<"reset vertex to: " << vertZ <<endl;
453  it_lastTrack->dca=dca.second;
454  it_lastTrack->ipZ=dca.first;
455 
456  // cout <<"get track10" <<endl;
457 
458  // cout <<" returning true " <<endl;
459  return true;
460  }
461  // cout <<" false... " <<endl;
462  return false;
463 };
464 
465 Double_t StFgtStraightTrackMaker::getRPhiRatio(vector<generalCluster>::iterator hitIterBegin, vector<generalCluster>::iterator hitIterEnd)
466 {
467  Short_t quad;
468  Char_t layer;
469  Int_t numR=0;
470  Int_t numPhi=0;
471  vector<generalCluster>::iterator hitIter=hitIterBegin;
472  for(;hitIter!=hitIterEnd;hitIter++)
473  {
474  quad=hitIter->quad;
475  layer=hitIter->layer;
476 
477  if(layer=='R')
478  numR++;
479  else
480  numPhi++;
481  }
482 
483  if(numR+numPhi>0)
484  return (numR-numPhi)/((Double_t)(numR+numPhi));
485  else
486  return -1;
487 };
488 
492 Int_t StFgtStraightTrackMaker::Make()
493 {
494  // cout <<"tracker! had " << m_tracks.size() << " old tracks" <<endl;
495  for(vector<AVTrack>::iterator it=m_tracks.begin();it!=m_tracks.end();it++)
496  {
497  if(it->points)
498  delete it->points;
499  }
500  m_tracks.clear();
501 
502  // cout <<"ave make " <<endl;
503  Int_t ierr = kStOk;
504  StFgtGeneralBase *fgtGenMkr = static_cast<StFgtGeneralBase * >( GetMaker("fgtGenBase"));
505 
506 
507  pClusters=fgtGenMkr->getClusters();
508  for(int i=0;i<6;i++)
509  {
510  // cout <<"there are " << pClusters[i]->size() << " clusters in disk " << i <<endl;
511  for(int j=0;j<pClusters[i]->size();j++)
512  {
513 
514  /* if((*(pClusters[i]))[j].layer=='R')
515  cout <<"R layer, ";
516  else
517  cout <<"Phi layer, ";*/
518 
519  Double_t posPhi=(*(pClusters[i]))[j].posPhi;
520  Double_t posR=(*(pClusters[i]))[j].posR;
521  Int_t clusSize=(*(pClusters[i]))[j].clusterSize;
522  Double_t charge=(*(pClusters[i]))[j].clusterCharge;
523  Int_t cntGeoId=(*(pClusters[i]))[j].centralStripGeoId;
524  Int_t seedType=(*(pClusters[i]))[j].seedType;
525  // cout <<"cluster pos phi: " << posPhi <<" posR: " << posR <<" clusSize: " << clusSize << " charge: "<< charge <<" geoId: "<< cntGeoId <<" seedT : " << seedType<<endl;
526  }
527  }
528  pStrips=fgtGenMkr->getStrips();
529 
530  // cout <<" mtracks size: " << m_tracks.size() << " oldnum: "<< oldNumTracks <<endl;
531 
532  Float_t x;
533  Float_t y;
534  // if(vtxRank<=1)
535  // return kStOk;
536 
537  //look at the r phi ratio for each disk
538  // for(int i=0;i<6;i++)
539  // {
540  // Double_t ratio=getRPhiRatio(pClusters[i]->begin(),pClusters[i]->end());
541  // rPhiRatioPlots[i]->Fill(ratio);
542  // cout << "ratio for disk: " << i << " is " << ratio <<" disk has: " << tmpClusterCol->getHitVec().size() << "hits" <<endl;
543  // }
544 
545  //vector<generalCluster> &hitVecD1=*(pClusters[0]);
546  //vector<generalCluster> &hitVecD6=*(pClusters[5]);
547  set<Int_t> usedPoints;//saves the points that have been used for tracks (and shouldn't be reused)
548  Double_t D1Pos=StFgtGeneralBase::getLocDiscZ(0);
549  Double_t D6Pos=StFgtGeneralBase::getLocDiscZ(5);
550  Double_t zArm=D6Pos-D1Pos;
551  vector<generalCluster>::reverse_iterator hitIterD1,hitIterD6;
552  vector<generalCluster>::iterator hitIterD1R, hitIterD6R, hitIter, hitIter2;
553  StFgtHit* fgtHitD1Phi;
554  StFgtHit* fgtHitD1R;
555  StFgtHit* fgtHitD6Phi;
556  StFgtHit* fgtHitD6R;
557 
558 
559  for(int locMinNumFitPoints=6;locMinNumFitPoints>=minNumFitPoints;locMinNumFitPoints--)
560  {
561  // cout <<"eff disk : " << m_effDisk <<endl;
562  for(int iSeed1=0;iSeed1<5;iSeed1++)
563  {
564  for(int iSeed2=iSeed1+1;iSeed2<6;iSeed2++)
565  {
566  if((iSeed2-iSeed1)<1)//to have large enough lever arm. Also, since three points are required shouldn't matter?
567  continue;
568  if(iSeed1==m_effDisk || iSeed2==m_effDisk)
569  continue;
570  if(pClusters[iSeed1]->size() > maxClusters || pClusters[iSeed2]->size() > maxClusters)
571  {
572  // cout <<"too many clusters in the disk!!!"<<endl<<endl;
573  continue;
574  }
575  //track where we have hits in disks
576  // cout <<"seed 1: " << iSeed1 <<" seed2: "<< iSeed2 <<endl;
577  vector<generalCluster> &hitVecSeed1=*(pClusters[iSeed1]);
578  vector<generalCluster> &hitVecSeed2=*(pClusters[iSeed2]);
579 
580  D1Pos=StFgtGeneralBase::getLocDiscZ(iSeed1);
581  D6Pos=StFgtGeneralBase::getLocDiscZ(iSeed2);
582  zArm=D6Pos-D1Pos;
583 
584  // for(hitIterD1=hitVecSeed1.begin();hitIterD1 != hitVecSeed1.end();hitIterD1++)
585  for(hitIterD1=hitVecSeed1.rbegin();hitIterD1 != hitVecSeed1.rend();hitIterD1++)
586  {
587  //this is from the loose clustering and the cluster doesn't have energy match
588  if(useChargeMatch && !hitIterD1->hasMatch)
589  continue;
590  Short_t quadP=hitIterD1->quad;
591  Char_t layer=hitIterD1->layer;
592 
593  Double_t seed1ChargePhi=hitIterD1->CHARGE_MEASURE;
594  Double_t seed1SizePhi=hitIterD1->clusterSize;
595  // if(quadP>2)
596  // continue;
597 
598  //do 1D 'fit' with r strips and the (x,y) thing
599  Int_t geoIdSeed1=hitIterD1->centralStripGeoId;
600  if(usedPoints.find(geoIdSeed1)!=usedPoints.end())
601  continue;
602  if(layer!='P')
603  continue;
604  // cout <<"ave make1 " <<endl;
605  Float_t phiD1=hitIterD1->posPhi;
606  fgtHitD1Phi=hitIterD1->fgtHit;
607  for(hitIterD6=hitVecSeed2.rbegin();hitIterD6 != hitVecSeed2.rend();hitIterD6++)
608  {
609  if(useChargeMatch && !hitIterD6->hasMatch)
610  continue;
611 
612  Int_t geoIdSeed2=hitIterD6->centralStripGeoId;
613  Short_t quadP_2=hitIterD6->quad;
614  //Short_t disc=hitIterD6->disc;
615  //Short_t strip=hitIterD6->strip;
616  Char_t layer=hitIterD6->layer;
617  Double_t seed2ChargePhi=hitIterD6->CHARGE_MEASURE;
618  Double_t seed2SizePhi=hitIterD6->clusterSize;
619  if(layer!='P')
620  continue;
621  if(usedPoints.find(geoIdSeed2)!=usedPoints.end())
622  continue;
623  Float_t phiD6=hitIterD6->posPhi;
624  fgtHitD6Phi=hitIterD6->fgtHit;
625  if(fabs(phiD6-phiD1)>maxPhiDiff)
626  continue;
627 
628  for(hitIterD1R=hitVecSeed1.begin();hitIterD1R != hitVecSeed1.end();hitIterD1R++)
629  {
630  if(useChargeMatch && !hitIterD1R->hasMatch)
631  continue;
632  Int_t geoIdSeed1R=hitIterD1R->centralStripGeoId;
633  Short_t quadR=hitIterD1R->quad;
634  //Short_t disc=hitIterD1R->disc;
635  //Short_t strip=hitIterD1R->strip;
636  Char_t layer=hitIterD1R->layer;
637  Double_t seed1ChargeR=hitIterD1R->CHARGE_MEASURE;
638  Double_t seed1SizeR=hitIterD1R->clusterSize;
639  Int_t quadSeed1=-1;
640  fgtHitD1R=hitIterD1R->fgtHit;
641  if(layer!='R')
642  continue;
643  if(quadR!=quadP)
644  continue;
645  quadSeed1=quadR;
646  if(usedPoints.find(geoIdSeed1R)!=usedPoints.end())
647  continue;
648 
649  Float_t rD1=hitIterD1R->posR;
650  Float_t xD1=rD1*cos(phiD1);
651  Float_t yD1=rD1*sin(phiD1);
652 
653  // cout <<"disk: " << iSeed1<<", phiD1: " << phiD1 <<" xD1: " << xD1 <<" yD1: " << yD1 <<" rD1: " << rD1 <<endl;
654  // cout <<"ave make3 " <<endl;
655 
656  for(hitIterD6R=hitVecSeed2.begin();hitIterD6R != hitVecSeed2.end();hitIterD6R++)
657  {
658  if(useChargeMatch && !hitIterD6R->hasMatch)
659  continue;
660  Int_t geoIdSeed2R=hitIterD6R->centralStripGeoId;
661  Short_t quadR_2=hitIterD6R->quad;
662  //Short_t disc=hitIterD6R->disc;
663  //Short_t strip=hitIterD6R->strip;
664  Char_t layer=hitIterD6R->layer;
665  Double_t seed2ChargeR=hitIterD6R->CHARGE_MEASURE;
666 
667  Double_t seed2SizeR=hitIterD6R->clusterSize;
668  Int_t quadSeed2=-1;
669  fgtHitD6R=hitIterD6R->fgtHit;
670  if(quadP_2!=quadR_2)
671  continue;
672  quadSeed2=quadP_2;
673  if(layer!='R')
674  continue;
675 
676 
677  if(usedPoints.find(geoIdSeed2R)!=usedPoints.end())
678  continue;
679  Float_t rD6=hitIterD6R->posR;
680  //track goes towards smaller radii
681  if(!isCosmic)
682  {
683  if(rD1>rD6)
684  continue;
685  }
686  vector<AVPoint>* v_points=new vector<AVPoint>;
687  vvPoints.push_back(v_points);
688 
689  //add the seed points to the points
690  Double_t xD6=rD6*cos(phiD6);
691  Double_t yD6=rD6*sin(phiD6);
692  // cout <<"Disk " << iSeed2 <<", phiD6: " << phiD6 <<" xD6: " << xD6 <<" yD6: " << yD6 <<" rD6: " << rD6 <<endl;
693  AVPoint avp1(xD1,yD1,D1Pos,rD1,phiD1,iSeed1,quadSeed1,seed1ChargeR, seed1ChargePhi, seed1SizeR, seed1SizePhi);
694  avp1.fgtHitR=fgtHitD1R;
695  avp1.fgtHitPhi=fgtHitD1Phi;
696  v_points->push_back(avp1);
697  AVPoint avp2(xD6,yD6,D6Pos,rD6,phiD6,iSeed2,quadSeed2,seed2ChargeR, seed2ChargePhi, seed2SizeR, seed2SizePhi);
698  avp2.fgtHitR=fgtHitD6R;
699  avp2.fgtHitPhi=fgtHitD6Phi;
700  v_points->push_back(avp2);
702 
703  int iFound=0;
704  int iFoundR=0;
705  // cout <<"ave make4 " <<endl;
706  //zarm is d6 position - D1
707  //Double_t xIpExp=xD1+(xD6-xD1)*(-D1Pos)/zArm;
708  //Double_t yIpExp=yD1+(yD6-yD1)*(-D1Pos)/zArm;
710  //Double_t zIpExpX0=(D6Pos-(xD6/xD1)*D1Pos)*1/(1-xD6/xD1);
712  //Double_t zIpExpY0=(D6Pos-yD6/yD1*D1Pos)*1/(1-yD6/yD1);
713  // cout <<"
714  //now find other points
715  vector<generalCluster>::iterator iterClosestPhi;
716  vector<generalCluster>::iterator iterClosestR;
717 
718  Double_t closestDist=999999;
719  Int_t closestQuad;
720  Int_t quadTestR=-1;
721  // cout <<"looking for more hits..." <<endl;
722  for(int iD=0;iD<6;iD++)
723  {
724  // cout <<"testting disk: " << iD <<endl;
725  if(iD==iSeed1 || iD==iSeed2 || iD==m_effDisk)
726  continue;
727 
728  // cout <<"looking at disk: " << iD <<" seed1: " << iSeed1 << " seed2: " << iSeed2 <<endl;
729  // cout <<"not seed " << endl;
730  Bool_t found=false;
731  Bool_t foundR=false;
732  //check for hit
733  Double_t diskZ=StFgtGeneralBase::getLocDiscZ(iD);
734  //expected
735 
736  Double_t xPosExp=xD1+(xD6-xD1)*(diskZ-D1Pos)/zArm;
737  Double_t yPosExp=yD1+(yD6-yD1)*(diskZ-D1Pos)/zArm;
738  // cout <<"hope to see something x: " << xPosExp <<" y: " << yPosExp <<endl;
739  Double_t rPosExp=rD1+(rD6-rD1)*(diskZ-D1Pos)/zArm;
740  vector<generalCluster> &hitVec=*(pClusters[iD]);
741  for(hitIter=hitVec.begin();hitIter!=hitVec.end();hitIter++)
742  {
743  if(useChargeMatch && !hitIter->hasMatch)
744  continue;
745  //do 1D 'fit' with r strips and the (x,y) thing
746  Int_t geoIdPhi=hitIter->centralStripGeoId;
747  Short_t quad=hitIter->quad;
748  Int_t quadTestPhi=quad;
749  Short_t disc=hitIter->disc;
750  Short_t strip=hitIter->strip;
751  Char_t layer=hitIter->layer;
752 
753  if(usedPoints.find(geoIdPhi)!=usedPoints.end())
754  continue;
755  if(layer!='P')
756  continue;
757  Float_t phi=hitIter->posPhi;
758  // cout <<"phi : " << phi << " phiD1: " << phiD1 <<" phiD6: " << phiD6 <<endl;
759  StFgtHit* fgtHitPhi=hitIter->fgtHit;
760  if(fabs(phi-phiD1)>maxPhiDiff)
761  continue;
762  if(fabs(phi-phiD6)>maxPhiDiff)
763  continue;
764 
765  // cout <<" survived max_phi_diff cuts " <<endl;
766  //Double_t phiCharge=hitIter->CHARGE_MEASURE;
767 
768  // Double_t phiCharge=hitIter->maxAdc;
769  // Int_t clusterSizePhi=hitIter->clusterSize;
770  // if(clusterSizePhi<=1)
771  // continue;
772  // cout <<"ave make5 " <<endl;
773  for(hitIter2=hitVec.begin();hitIter2!=hitVec.end();hitIter2++)
774  {
775  if(useChargeMatch && !hitIter2->hasMatch)
776  continue;
777  Int_t geoIdR=hitIter2->centralStripGeoId;
778  StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);//ok
779  // cout <<" r? " <<endl;
780  if(usedPoints.find(geoIdR)!=usedPoints.end())
781  continue;
782  // cout <<"not used yet " <<endl;
783  if(layer!='R')
784  continue;
785  quadTestR=quad;
786  if(quadTestR!=quadTestPhi)
787  continue;
788  Float_t r=hitIter2->posR;
789  StFgtHit* fgtHitR=hitIter->fgtHit;
790  //make sure that the radius makes sense for a track that goes from inner to outer radious
791  if(!isCosmic)
792  {
793  if(r>rD6 || r<rD1)
794  continue;
795  }
796  x=r*cos(phi);
797  y=r*sin(phi);
798  // cout <<"checking with x: " << x << " y: " << y << " phi: " << phi <<" r: " << r <<endl;
799  // cout <<" x, y: " << x <<", " << y << " exp: " << xPosExp << " , " << yPosExp <<endl;
800  Double_t dist2=(x-xPosExp)*(x-xPosExp)+(y-yPosExp)*(y-yPosExp);
801  // cout <<" dist2: " << dist2 <<endl;
802 
803  if(doAddMultiplePoints)
804  {
805  if(dist2<maxDist2)
806  {
807  Double_t rCharge=hitIter2->CHARGE_MEASURE;
808  Double_t phiCharge=hitIter->CHARGE_MEASURE;
809  Int_t clusterSizeR=hitIter2->clusterSize;
810  Int_t clusterSizePhi=hitIter->clusterSize;
811  AVPoint avp(x,y,diskZ,r,phi,iD,quadTestR, rCharge,phiCharge, clusterSizeR,clusterSizePhi);
812  avp.fgtHitR=fgtHitR;
813  avp.fgtHitPhi=fgtHitPhi;
814  v_points->push_back(avp);
815  }
816  }
817 
818  if(dist2<closestDist)
819  {
820  closestDist=dist2;
821  closestQuad=quadTestR;
822  iterClosestPhi=hitIter;
823  iterClosestR=hitIter2;
824  }
825  }
826  }
827  // cout <<"closest dist for disk: "<< iD <<" is : " << closestDist <<endl;
828  // cout <<"ave make6 " <<endl;
829  if(closestDist<1000 && closestDist<maxDist2)
830  {
831  found=true;
832  // cout <<"accepted " <<endl;
833  double r=iterClosestR->posR;
834  double phi=iterClosestPhi->posPhi;
835  StFgtHit* fgtHitR=iterClosestR->fgtHit;
836  StFgtHit* fgtHitPhi=iterClosestPhi->fgtHit;
837 
838  Int_t geoIdR=iterClosestR->centralStripGeoId;
839  Int_t geoIdPhi=iterClosestPhi->centralStripGeoId;
840 
841  double x=r*cos(phi);
842  double y=r*sin(phi);
843  // cout<<" adding point with r: "<< r <<" phi: " << phi <<" x: " << x <<" y: " << y <<endl;
844  Double_t rCharge=iterClosestR->CHARGE_MEASURE;
845  Double_t phiCharge=iterClosestPhi->CHARGE_MEASURE;
846  Int_t clusterSizeR=iterClosestR->clusterSize;
847  Int_t clusterSizePhi=iterClosestPhi->clusterSize;
848  // cout <<"charge R of middle disk: "<< iD<<": "<< rCharge <<" phicharge: " << phiCharge<<endl;
849 
850 // already added before if this flag is set
851  if(!doAddMultiplePoints)
852  {
853  AVPoint avp(x,y,diskZ,r,phi,iD,closestQuad, rCharge,phiCharge, clusterSizeR,clusterSizePhi);
854  // cout<<" adding point with r: "<< r <<" phi: " << phi <<" x: " << x <<" y: " << y <<endl;
855  avp.fgtHitR=fgtHitR;
856  avp.fgtHitPhi=fgtHitPhi;
857  v_points->push_back(avp);
858  }
859 
860  }
861  // cout <<"ave make8 " <<endl;
862  //only one per disk
863  if(found)
864  iFound++;
865  else
866  {
867 
868  }
869  if(foundR)
870  iFoundR++;
871  else
872  {
873  // cout <<"failed to find r " << rPosExp<<endl;
874  // v_rFail.push_back(pair<Int_t, Double_t>(iD,rPosExp));
875  }
876  closestDist=999999;
877 
878  }
879 
880  // cout <<"ave make9 " <<endl;
881  // if(iFound>=2 && fabs(xIpExp)<20 && fabs(yIpExp)<20 && fabs(zIpExpX0)<40 && fabs(zIpExpY0)<40) //at least one hit plus seed and pointing roughly to the interaction region
882  //with hard cuts on the track we can get the whole vertex distribution
883  // cout << " we found " << iFound <<" points " << endl;
884  if(iFound>=(locMinNumFitPoints-2)) //at least one hit plus seed and pointing roughly to the interaction region
885  {
886  // cout <<"found " <<endl;
887  Bool_t validTrack=false;
888  Bool_t passQCuts=true;
889  // if(v_x.size()>iFound)
890  {
891  Double_t ipZ;
892  // cout <<"about to get track" <<endl;
893  // cout <<"check for valid track " << endl;
894  //this also manipulates the points vector to only put points in there that fit
895  //get Track gets the ipZ, so it is ok to put in a random value here, but it is not a reference! So you don't get ipZ back
896  validTrack=getTrack(*v_points, ipZ);
897 
898 
899  if(m_tracks.size()>0)
900  {
901  passQCuts=trackQCuts(m_tracks.back());
902  }
903 
904  if(validTrack && passQCuts)
905  {
907  (m_tracks.back()).points=v_points;
908  //at least don't duplicate seeds..., the other ones might belong to multiple tracks
909  usedPoints.insert(geoIdSeed1);
910  usedPoints.insert(geoIdSeed1R);
911  usedPoints.insert(geoIdSeed2);
912  usedPoints.insert(geoIdSeed2R);
913  //put the rest of it in as well, doesn't matter that we duplicate the seeds, this is a set<int>
914  for(vector<AVPoint>::iterator it=v_points->begin();it!=v_points->end();it++)
915  {
916  usedPoints.insert(it->geoIdR);
917  usedPoints.insert(it->geoIdPhi);
918  }
919 
920 
921  }
922  else
923  {
924  v_points->clear();
925  m_tracks.pop_back();
926  delete v_points;
927  }
928  }
929 
930 
931  hitCounter++;
932  }
933  //start over
934  iFound=0;
935  iFoundR=0;
936 
937  }
938 
939  }
940 
941  }
942 
943  }
944 
945  }
946 
947  }
948  }
949  // cout <<"we have " << m_tracks.size()-oldNumTracks<< " tracks in this event " <<endl;
950  // cout <<"oldNumtracsk: " << oldNumTracks << " new size: " << m_tracks.size() <<endl;
951  // numTracks->Fill(m_tracks.size()-oldNumTracks);
952 
953  return ierr;
954 
955  };
956 
957 bool StFgtStraightTrackMaker::trackQCuts(AVTrack& trk)
958 {
959 
960  if(trk.chi2>maxChi2 || trk.trkZ> vertexCutPos || trk.trkZ< vertexCutNeg|| trk.dca> dcaCut )
961  {
962  return false;
963  }
964  return true;
965 };
966 
967 StFgtStraightTrackMaker::StFgtStraightTrackMaker( const Char_t* name): StMaker( name ),useChargeMatch(false),runningEvtNr(0),hitCounter(0),hitCounterR(0),maxChi2(2),dcaCut(1),vertexCutPos(70),vertexCutNeg(-120)
968 {
969  // cout <<"AVE constructor!!" <<endl;
970  int numTb=7;
971  m_effDisk=10;//
972  isMuDst=true; //might want to change this...
973  maxPhiDiff=0.1;
974  maxClusters=10;
975 
976  doFitWithVertex=false;
977  doRefitWithVertex=false;
978  doAddMultiplePoints=false;
979  maxDist2=1.0;
980  minNumFitPoints=3;
981 
982 };
983 
984 StFgtStraightTrackMaker::~StFgtStraightTrackMaker()
985 {
986 
987  //delete histogram arrays
988 };
989 
990 Int_t StFgtStraightTrackMaker::Finish(){
991  // cout<<" straight tracker finish" <<endl;
992  // cout <<" closing txt file " << endl;
993  // gStyle->SetPalette(1);
994  // cout <<"AVE finish function " <<endl;
995  Int_t ierr = kStOk;
996 
998  // vector<AVTrack>::iterator it=m_tracks.begin();
999  // cout <<"we found " << m_tracks.size() <<" tracks" <<endl;
1000  int counter=0;
1001 
1004  // cout <<"st done " <<endl;
1005  return ierr;
1006 };
1007 
1008 // construct histograms
1009 Int_t StFgtStraightTrackMaker::Init(){
1010  Int_t ierr = kStOk;
1011  return ierr;
1012 };
1013 ClassImp(StFgtStraightTrackMaker);
StFgtStraightTrackMaker::findClosestPoint
Double_t findClosestPoint(float mx, float bx, float my, float by, double xE, double yE, Int_t iD)
Definition: StFgtStraightTrackMaker.cxx:119
StFgtStraightTrackMaker::getDca
pair< double, double > getDca(vector< AVTrack >::iterator it)
Definition: StFgtStraightTrackMaker.cxx:53
StFgtStraightTrackMaker::getQuadFromCoo
Short_t getQuadFromCoo(Double_t x, Double_t y)
this is too naive..., assumes non-rotated quads
Definition: StFgtStraightTrackMaker.cxx:181
StFgtStraightTrackMaker::getTrack
Bool_t getTrack(vector< AVPoint > &points, Double_t ipZ)
Definition: StFgtStraightTrackMaker.cxx:259
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5