StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFgtAlignmentMaker.cxx
1 
11 #include <stdio.h>
12 #include <math.h>
13 
14 #include "StFgtAlignmentMaker.h"
15 #include "StEventTypes.h"
16 #include "StMessMgr.h"
17 #include "StDcaGeometry.h"
18 #include "TNtuple.h"
19 #include "TMinuit.h"
20 #include "TRandom.h"
21 #include "TCanvas.h"
22 #include "TText.h"
23 #include "TGraph.h"
24 #include "TH1F.h"
25 #include "TH2F.h"
26 #include "TArrayL64.h"
27 #include "StThreeVectorF.hh"
28 #include "StPhysicalHelix.hh"
29 #include "SystemOfUnits.h"
30 #include "THelixTrack.h"
31 #include "StDetectorName.h"
32 #include "fgtAlignment.h"
33 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
34 #include "StTpcDb/StTpcDb.h"
35 #include "StDetectorDbMaker/St_vertexSeedC.h"
36 #include "StRoot/StFgtDbMaker/StFgtDbMaker.h"
37 #include "StRoot/StFgtDbMaker/StFgtDb.h"
38 
39 #include "StMuDSTMaker/COMMON/StMuDst.h"
40 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
41 #include "StMuDSTMaker/COMMON/StMuDst.h"
42 #include "StMuDSTMaker/COMMON/StMuEvent.h"
43 #include "StarClassLibrary/StThreeVectorF.hh"
44 
45 #include "StRoot/StFgtPool/StFgtClusterTools/StFgtGeneralBase.h"
46 #include "StRoot/StFgtPool/StFgtClusterTools/StFgtStraightTrackMaker.h"
47 
48 #include "StRoot/StEEmcPool/StEEmcIUPi0/StEEmcIUPointMaker.h"
49 #include "StRoot/StEmcUtil/geometry/StEmcGeom.h"
50 
51 static const int mDebug=0; // debug mesasge level
52 static const int mEventDisp=0; // event display
53 
54 static const int MAXTRK=500000; //max number of tracks
55 static const int MAXHIT=8; //max number of hits
56 static const int NPAR=24*6; //number of parameters=(6 discs)*(4 quads)*(3 angles + 3 offset)
57 
58 static const int NAXIS=4; // number of axis for hist
59 static TH1F *hist1[kFgtNumDiscs+6][kFgtNumQuads][NAXIS]; // 1d residual histos
60 static TH2F *hist2[kFgtNumDiscs+6][kFgtNumQuads][NAXIS*2]; // 2d residual histos
61 static TH1F *histdz;
62 
63 static const double PI = 4.0*atan(1);
64 
65 const float ISOCUT=0.5;
66 const float maxdfgt=0.15, maxdvtx=0.5, maxdtpc=10.0, maxdemc=10.0;
67 const float maxpfgt=0.005, maxpvtx=0.5, maxptpc= 0.3, maxpemc=0.3;
68 
69 struct TRKHIT{
70  int run,evt,nhit,nhitUse,nhitFgt,nhitVtx,nhitTpc,nhitPrompt,nhitEemc,used;
71  float dca, chi2, trkz, dz, eta, phi, opt; //track parameters (last one is 1/pT)
72  int det[MAXHIT]; //0-23=FGT(disc*4+quad), 24-27=vertex, 28-31=TPC, 32-35=Prompt, 36-39=EEMC separated by quad
73  float x[MAXHIT],y[MAXHIT],z[MAXHIT]; //xyz
74  float ex[MAXHIT],ey[MAXHIT],ez[MAXHIT]; //xyz error
75  float lr[MAXHIT],lp[MAXHIT],r[MAXHIT],p[MAXHIT]; //local and globak r and phi
76  float dx[MAXHIT],dy[MAXHIT],dr[MAXHIT],dp[MAXHIT]; //residuals
77  float tw[MAXHIT],p1[MAXHIT],p2[MAXHIT],po[MAXHIT]; //EEMC energies (FGT: tw=ChgPhi, p1=ChgR, p2=PhiRAsy, po=EvenOddAsy)
78  float su[MAXHIT],sv[MAXHIT]; //ESMD energies
79  int nrl[MAXHIT],nsu[MAXHIT],nsv[MAXHIT]; //eemc number of relatives, strips in smd u & v
80  float ele[4],trk[4],bemc[4],eemc[4],tot[4],rec[4],iso[4]; //4 momentums of ele/tracks/bemc/eemc/total(trk+bemc+eemc)/recoil(=tot-ele)/isolation cone
81  float ptbal,riso; //pt balance & energy ratio ele/total in isolation cone
82  bool use[MAXHIT]; //false=do not use in fit true=use in fit
83 };
84 
85 static int mNtrk[kFgtNumQuads]; // number of tracks per quad
86 static int mNtrkUse[kFgtNumQuads]; // number of usable tracks per quad after hitmask
87 static TRKHIT mHit[kFgtNumQuads][MAXTRK]; // storage to keep hit for track
88 static TRKHIT mHitAlg; // current working hits to be modified with aliggnment parameters
89 static int mQuad; // quad currently working on
90 static int mStep; // step number
91 static int mHitMaskDisc; // disc mask to define hits to be used hits in track fits. bit0-5=fgt, bit6=vtx, bit7=tpc, bit8=prompt, bit9=eemc
92 static int mResidMaskDisc; // disc mask to include residual sum of all tracks for alignments
93 static int mFgtInputRPhi; // 0=x/y/z or 1=r/phi/z for fgt hits in mHit
94 static int mFillHist; // 0=not filling histo 1=filling histo before alignments 2=filling histo after alignment
95 static int mReqHit=0; // required # of hits
96 static int mReqFgtHit=0; // required # of FGT hits
97 static int mReqVtx=0; // required vertex
98 static int mReqTpcHit=0; // required # of TPC hits
99 static int mReqPromptHit=0; // required # of TPC prompt hits
100 static int mReqEemcHit=0; // required # of EEMC hits
101 static float mErrFGT=0;
102 static float mErrVTX=0;
103 static float mErrVTXZ=0;
104 static float mErrTPCI=0;
105 static float mErrTPCO=0;
106 static float mErrTPCZ=0;
107 static float mErrPPT=0;
108 static float mErrEMC=0;
109 static int mBeamLine=1;
110 
111 static StFgtDb* mDb;
112 static fgtAlignment_st* orig_algpar;
113 
114 static StThreeVectorF ele,tracks,bemcs,eemcs,tot,rec,iso;
115 static float ptbal,riso;
116 
117 #define SMALL_NUM 1e-4
118 
119 class Vector{
120 public:
121  Vector(double vx=0, double vy=0, double vz=0):x(vx),y(vy),z(vz){}
122  Vector(const Vector& v):x(v.x),y(v.y),z(v.z){}
123  Vector(const Vector& v1, const Vector& v2):x(v1.x-v2.x),y(v1.y-v2.y),z(v1.z-v2.z){}
124  // Vector operator=(const Vector& v) {return Vector(v.x,v.y,v.z);}
125  Vector operator=(const Vector& v) {x=v.x; y=v.y; z=v.z;}
126  Vector operator-() {return Vector(-x,-y,-z);}
127  friend Vector operator+(const Vector& v1, const Vector& v2) {return Vector(v1.x+v2.x,v1.y+v2.y,v1.z+v2.z);}
128  friend Vector operator-(const Vector& v1, const Vector& v2) {return Vector(v1.x-v2.x,v1.y-v2.y,v1.z-v2.z);}
129  friend Vector operator*(const double d, const Vector& v) {return Vector(d*v.x,d*v.y,d*v.z);}
130  friend Vector operator*(const Vector& v, const double d) {return Vector(d*v.x,d*v.y,d*v.z);}
131  friend double operator*(const Vector& v1, const Vector& v2) {return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z;}
132  double length() {return sqrt((*this)*(*this));}
133  void print(){printf("x=%8.3f y=%8.3f z=%8.3f\n",x,y,z);}
134  double x,y,z;
135 };
136 
137 class Line{
138 public:
139  Line(Vector& u0, Vector& u1):v0(u0),v1(u1) {};
140  Vector v0,v1;
141 };
142 
143 double distance(Line l1, Line l2, Vector& p1, Vector& p2){
144  Vector u(l1.v1,l1.v0);
145  Vector v(l2.v1,l2.v0);
146  Vector w(l1.v0,l2.v0);
147  double a=u*u;
148  double b=u*v;
149  double c=v*v;
150  double d=u*w;
151  double e=v*w;
152  double D=a*c - b*b;
153  double sc,tc;
154  if(D<SMALL_NUM) { // the lines are almost parallel
155  sc = 0.0;
156  tc = (b>c ? d/b : e/c); // use the largest denominator
157  printf("SMALL NUM!!!\n");
158  }else{
159  sc = (b*e - c*d) / D;
160  tc = (a*e - b*d) / D;
161  }
162  // 2 closest point
163  p1=l1.v0 + (sc * u);
164  p2=l2.v0 + (tc * v);
165  // get the difference of the two closest points
166  Vector dP = p1 - p2;
167  //Vector dP = w + (sc * u) - (tc * v); // = L1(sc) - L2(tc)
168  return dP.length(); // return the closest distance
169 }
170 
171 void getZvtxAndDca(double* p, double& dca, double& z, double zmin=-900, double zmax=900){
172  double x0=0,y0=0,dxdz=0,dydz=0;
173  if(mFgtInputRPhi>0){
174  St_vertexSeedC* vSeed = St_vertexSeedC::instance();
175  x0=vSeed->x0(); dxdz=vSeed->dxdz();
176  y0=vSeed->y0(); dydz=vSeed->dydz();
177  }
178  Vector p1(p[0]+p[1]*zmin, p[2]+p[3]*zmin, zmin);
179  Vector p2(p[0]+p[1]*zmax, p[2]+p[3]*zmax, zmax);
180  Vector p3(x0+dxdz*zmin, y0+dydz*zmin, zmin);
181  Vector p4(x0+dxdz*zmax, y0+dydz*zmax, zmax);
182  Vector p5,p6;
183  Line l1(p1,p2), l2(p3,p4);
184  dca=distance(l1,l2,p5,p6); //3d DCA length
185  z=p6.z; //z of DCA point on beamline constrain
186  //printf("BEAM=%8.2f %8.2f %8.2f FGT=%8.2f %8.2f %8.2f\n",p6.x,p6.y,p6.z,p5.x,p5.y,p5.z);
187 }
188 
189 inline StPhysicalHelix getStPhysicalHelix(double x0, double y0, double curve, double dip, double slp){
190  //printf("x0=%8.4f y0=%8.4f curv=%8.4f dip=%8.4f slp=%8.4f -> ",x0,y0,curve,dip,slp);
191  StThreeVectorD o(x0,y0,0);
192  int sign=int(curve/fabs(curve));
193  double phase=slp - sign*PI/2.0;
194  StPhysicalHelix h(fabs(curve),dip,phase,o,sign);
195  //printf("GetHelix c=%10.7f sign=%2d h=%2d\n",curve,sign,h.h());
196  //StThreeVectorD t1=h.at(0.0);
197  //printf(" x1=%8.4f y1=%8.4f z1=%8.4f \n",t1.x(),t1.y(),t1.z());
198  //StThreeVectorD t2=h.at(100.0);
199  //printf(" x2=%8.4f y2=%8.4f z2=%8.4f \n",t2.x(),t2.y(),t2.z());
200  //StThreeVectorD t3=h.at(200.0);
201  //printf(" x3=%8.4f y3=%8.4f z3=%8.4f \n",t3.x(),t3.y(),t3.z());
202  return h;
203 };
204 
205 inline StPhysicalHelix getStPhysicalHelix(double *p) {return getStPhysicalHelix(p[0],p[1],p[2],p[3],p[4]);}
206 
207 inline StPhysicalHelix getStPhysicalHelix(double x0, double y0, double dxdz, double dydz, double curve, int sign) {
208  double dr=sqrt(dxdz*dxdz+dydz*dydz);
209  double dip = atan2(1.0,dr);
210  double slp = atan2(dydz,dxdz);
211  return getStPhysicalHelix(x0,y0,curve*sign,dip,slp);
212 }
213 
214 void getZvtxAndDcaFromHelix(double* par, double& dca, double& z, double &eta, double &phi, double& opt){
215  St_vertexSeedC* vSeed = St_vertexSeedC::instance();
216 
217  //beamline
218  //double xyz[3]={vSeed->x0(),vSeed->y0(), 0.0};
219  //double dir[3]={vSeed->dxdz(),vSeed->dydz(),1.0};
220  //if(mFgtInputRPhi==0){ xyz[0]=0.0; xyz[1]=0.0; xyz[2]=0.0; dir[0]=0.00001; dir[1]=0.00001; dir[2]=1;}
221  //THelixTrack bb(xyz,dir,0.0000000001);
222  //bb.Move(0.0); printf(" BLT x1=%8.4f y1=%8.4f z1=%8.4f \n",bb.Pos()[0],bb.Pos()[1],bb.Pos()[2]);
223  //bb.Move(100.0); printf(" BLT x2=%8.4f y2=%8.4f z2=%8.4f \n",bb.Pos()[0],bb.Pos()[1],bb.Pos()[2]); bb.Move(-100);
224  //bb.Move(z); printf(" BLT x3=%8.4f y3=%8.4f z3=%8.4f \n",bb.Pos()[0],bb.Pos()[1],bb.Pos()[2]); bb.Move(-z);
225 
226  //track
227  double xyz2[3]={par[0], par[1], 0.0};
228  double dir2[3]={cos(par[4])/tan(par[3]),sin(par[4])/tan(par[3]),1.0};
229  THelixTrack hh(xyz2,dir2,par[2]);
230  //hh.Move(0.0); printf(" TRT x1=%8.4f y1=%8.4f z1=%8.4f \n",hh.Pos()[0],hh.Pos()[1],hh.Pos()[2]);
231  //hh.Move(100.0); printf(" TRT x2=%8.4f y2=%8.4f z2=%8.4f \n",hh.Pos()[0],hh.Pos()[1],hh.Pos()[2]); hh.Move(-100);
232  //hh.Move(z); printf(" TRT x3=%8.4f y3=%8.4f z3=%8.4f \n",hh.Pos()[0],hh.Pos()[1],hh.Pos()[2]); hh.Move(-z);
233 
234  //get DCA
235  double x0=0, y0=0, z0=0;
236  if(mFgtInputRPhi>0){ x0=vSeed->x0(); y0=vSeed->y0();}
237  double ss3=0;
238  for(int i=0; i<5; i++){
239  ss3=hh.Path(x0,y0); //closest point to x0,y0 line
240  hh.Move(ss3);
241  z0=hh.Pos()[2];
242  if(mFgtInputRPhi>0){
243  x0=vSeed->x0()+vSeed->dxdz()*z0;
244  y0=vSeed->y0()+vSeed->dydz()*z0;
245  }
246  double x=hh.Pos()[0], y=hh.Pos()[1], z=hh.Pos()[2];
247  double dx=x-x0, dy=y-y0, dz=z-z0;
248  dca=sqrt(dx*dx+dy*dy+dz*dz);
249  //printf("length=%8.4f x0=%8.4f y0=%8.4f z0=%8.4f : x=%8.4f y=%8.4f z=%8.4f dca=%8.4f\n",ss3,x0,y0,z0,hh.Pos()[0],hh.Pos()[1],hh.Pos()[2],dca);
250  hh.Move(-ss3);
251  }
252  z=z0; //z of DCA point on beamline constrain
253  double bfield=0.5*tesla;
254  StPhysicalHelix h = getStPhysicalHelix(par);
255  StThreeVectorD p = h.momentumAt(ss3,bfield);
256  //StThreeVectorD x = h.at(ss3);
257  //printf("ST x=%8.4f y=%8.4f z=%8.4f\n",x.x(),x.y(),x.z());
258  int sign=h.charge(bfield);
259  eta=p.pseudoRapidity();
260  phi=p.phi();
261  opt=sign/p.mag();
262 
263  /*doesn't seems to work
264  //make StPhysicalHelix from fitted parameters
265  StPhysicalHelix h = getStPhysicalHelix(par);
266  //make StPhysicalHelix from beamLine with very small curvature
267  StPhysicalHelix b;
268  if(mFgtInputRPhi==0){
269  b = getStPhysicalHelix(0,0,0.000001,0.000001,0.0000000001,1);
270  }else{
271  b = getStPhysicalHelix(vSeed->x0(),vSeed->y0(),vSeed->dxdz(),vSeed->dydz(),0.0000000001,1);
272  }
273 
274  StThreeVectorD t1=b.at(0.0);
275  printf(" BL x1=%8.4f y1=%8.4f z1=%8.4f \n",t1.x(),t1.y(),t1.z());
276  StThreeVectorD t2=b.at(100.0);
277  printf(" BL x2=%8.4f y2=%8.4f z2=%8.4f \n",t2.x(),t2.y(),t2.z());
278  StThreeVectorD t3=b.at(z);
279  printf(" BL x3=%8.4f y3=%8.4f z3=%8.4f \n",t3.x(),t3.y(),t3.z());
280 
281  StThreeVectorD s1=h.at(0.0);
282  printf(" TR x1=%8.4f y1=%8.4f z1=%8.4f \n",s1.x(),s1.y(),s1.z());
283  StThreeVectorD s2=h.at(100.0);
284  printf(" TR x2=%8.4f y2=%8.4f z2=%8.4f \n",s2.x(),s2.y(),s2.z());
285  StThreeVectorD s3=h.at(z);
286  printf(" TR x3=%8.4f y3=%8.4f z3=%8.4f d=%8.4f\n",s3.x(),s3.y(),s3.z(),abs(s3-t3));
287 
288  //get DCA
289  pair<double, double> d=b.pathLengths(h);
290  printf("path %8.4f %8.4f\n",d.first,d.second);
291  StThreeVectorD d1 = b.at(d.first);
292  StThreeVectorD d2 = h.at(d.second);
293  StThreeVectorD d3 = d1-d2;
294  dca=d3.mag(); //3d dca
295  z=d2.z(); //z of DCA point on beamline constrain
296  double bfield=0.5*tesla;
297  StThreeVectorD p = h.momentumAt(d.second,bfield);
298  int sign=h.charge(bfield);
299  eta=p.pseudoRapidity();
300  phi=p.phi();
301  opt=sign/p.mag();
302  //printf("getZvtxAndDcaFromHelix c=%10.6f sign=%2d h=%2dcurv=%10.6f mag=%10.6f opt=%10.6f\n",par[2],sign,h.h(),h.curvature(),p.mag(), opt);
303  */
304 }
305 
306 //Move fgt hits with alignment parameter
307 void getAlign(int itrk, fgtAlignment_st* algpar){
308  memcpy(&mHitAlg,&mHit[mQuad][itrk],sizeof(mHitAlg));
309  if(mFgtInputRPhi==0) return; //skip alignment for fake data
310  if(mDebug>3) printf("getAlign quad=%1d trk=%d mHitAlg.nhit=%d %d\n",mQuad,itrk,mHitAlg.nhit,mHit[mQuad][itrk].nhit);
311  for(int i=0; i<mHitAlg.nhit; i++){
312  if(mHitAlg.det[i]<24){
313  double r=mHitAlg.lr[i];
314  double phi=mHitAlg.lp[i];
315  int disc=mHitAlg.det[i]/4;
316  int quad=mHitAlg.det[i]%4;
317  TVector3 xyz,gxyz;
318  mDb->getStarXYZ(disc,quad,r,phi,xyz,2,algpar); //fgt internal alignment
319  //printf("y local=%8.2f\n",xyz.Y());
320  if(gStTpcDb){
321  double local[3]={0,0,0}, global[3]={0,0,0};
322  xyz.GetXYZ(local);
323  TGeoHMatrix globalMatrix = gStTpcDb->Tpc2GlobalMatrix(); //TPC global offsets
324  globalMatrix.LocalToMaster(local,global);
325  gxyz.SetXYZ(global[0],global[1],global[2]);
326  //printf("FGT local %9.6f %9.6f %9.6f\n",local[0],local[1],local[2]);
327  //printf("FGT globl %9.6f %9.6f %9.6f\n",global[0],global[1],global[2]);
328  //globalMatrix.Print();
329  }else{
330  static int nmess=0;
331  if(nmess<100){
332  printf("StFgtAlignmentMaker::Make could not get gStTpcDb... global xyz is same as fgt local xyz\n");
333  nmess++;
334  }
335  gxyz=xyz;
336  }
337  mHitAlg.x[i]=gxyz.X();
338  mHitAlg.y[i]=gxyz.Y();
339  mHitAlg.z[i]=gxyz.Z();
340  mHitAlg.r[i]=sqrt(mHitAlg.x[i]*mHitAlg.x[i]+mHitAlg.y[i]*mHitAlg.y[i]);
341  mHitAlg.p[i]=atan2(mHitAlg.y[i],mHitAlg.x[i]);
342  //printf("after y=%8.2f\n",mHitAlg.y[i]);
343  }
344  if(mDebug>3) printf("mHit[%3d] x=%8.3f y=%8.3f z=%8.3f det=%2d\n",i,
345  mHitAlg.x[i],mHitAlg.y[i],mHitAlg.z[i],mHitAlg.det[i]);
346  }
347 }
348 
349 //fill residual histos
350 void fillHist(int disc, int quad, float dx, float dy, float dr, float dp, float x, float y, float r, float p){
351  hist1[disc][quad][0]->Fill(dx);
352  hist1[disc][quad][1]->Fill(dy);
353  hist1[disc][quad][2]->Fill(dr);
354  hist1[disc][quad][3]->Fill(dp);
355  hist2[disc][quad][0]->Fill(x,dx);
356  hist2[disc][quad][1]->Fill(x,dy);
357  hist2[disc][quad][2]->Fill(r,dr);
358  hist2[disc][quad][3]->Fill(r,dp);
359  hist2[disc][quad][4]->Fill(y,dx);
360  hist2[disc][quad][5]->Fill(y,dy);
361  hist2[disc][quad][6]->Fill(p,dr);
362  hist2[disc][quad][7]->Fill(p,dp);
363  if(mDebug>3) cout << Form("FillHist d=%1d q=%1d dx=%8.3f dy=%8.3f dr=%8.3f dp=%9.5f",
364  disc,quad,dx,dy,dr,dp) <<endl;
365  //if(disc>6) cout << Form("FillHist d=%1d q=%1d dx=%8.3f dy=%8.3f dr=%8.3f dp=%9.5f",
366  // disc,quad,dx,dy,dr,dp) <<endl;
367 }
368 
369 //event display
370 void eventDisplay(double *par){
371  //setup helix track
372  StPhysicalHelix h = getStPhysicalHelix(par);
373  static StThreeVectorD n(0,0,1);
374  //set up TGraphs
375  TGraph *hitxz=new TGraph(mHitAlg.nhit); hitxz->SetMarkerStyle(21); hitxz->SetMarkerSize(0.4); hitxz->SetMarkerColor(kBlue);
376  TGraph *hityz=new TGraph(mHitAlg.nhit); hityz->SetMarkerStyle(21); hityz->SetMarkerSize(0.4); hityz->SetMarkerColor(kBlue);
377  TGraph *hitpr=new TGraph(mHitAlg.nhit); hitpr->SetMarkerStyle(21); hitpr->SetMarkerSize(0.4); hitpr->SetMarkerColor(kBlue);
378  //Plot hits
379  for(int i=0; i<mHitAlg.nhit; i++){
380  StThreeVectorD z(0,0,mHitAlg.z[i]);
381  double s=h.pathLength(z,n);
382  double x=h.x(s);
383  double y=h.y(s);
384  double dx = x - mHitAlg.x[i];
385  double dy = y - mHitAlg.y[i];
386  //double dr = sqrt(x*x+y*y) - sqrt(mHitAlg.x[i]*mHitAlg.x[i]+mHitAlg.y[i]*mHitAlg.y[i]);
387  double dp = atan2(y,x) - atan2(mHitAlg.y[i],mHitAlg.x[i]);
388  if(dp>PI) dp-=2*PI;
389  if(dp<-PI) dp+=2*PI;
390  double hx=mHitAlg.x[i];
391  double hy=mHitAlg.y[i];
392  double hz=mHitAlg.z[i];
393  double hr=sqrt(hx*hx+hy*hy);
394  //double hp=atan2(hy,hx);
395  hitxz->SetPoint(i,hz,dx);
396  hityz->SetPoint(i,hz,dy);
397  if(fabs(dp)<0.1) hitpr->SetPoint(i,hr,dp);
398  //if(mDebug>3) printf("residHelix %3d dx=%12.8f f=%12.8\n",i,dx,dy,f);
399  }
400  TCanvas *c1=new TCanvas("EventDisp","Event Disp",0,0,800,800);
401  c1->Divide(2,2);
402  c1->cd(1); hitxz->Draw("AP");
403  c1->cd(2); hityz->Draw("AP");
404  c1->cd(3); hitpr->Draw("AP");
405  c1->Update();
406  cin.get();
407 }
408 
409 //Helix fit, called from TMinuit
410 void fitHelix(Int_t &npar, Double_t* gin, Double_t &f, Double_t *par, Int_t iflag){
411  if(mDebug>4) printf("fitHelix nhit=%d\n",mHitAlg.nhit);
412  f=0;
413  StPhysicalHelix h = getStPhysicalHelix(par);
414  static StThreeVectorD n(0,0,1);
415  int ifvtx=0;
416  int ndf=0;
417  for(int i=0; i<mHitAlg.nhit; i++){
418  if(mHitAlg.use[i]==false) continue;
419  if(mHitAlg.det[i]/4==6) ifvtx=1;
420  StThreeVectorD z(0,0,mHitAlg.z[i]);
421  double s=h.pathLength(z,n);
422  double dx = h.x(s) - mHitAlg.x[i];
423  double dy = h.y(s) - mHitAlg.y[i];
424  f += (dx*dx+dy*dy)/(mHitAlg.ex[i]*mHitAlg.ex[i]+mHitAlg.ey[i]*mHitAlg.ey[i]);
425  ndf+=2;
426  if(mDebug>5) printf("Helix1 Curve=%12.9f x=%8.3f y=%8.3f dx=%8.4f dy=%8.4f f=%8.4fe\n",par[2],h.x(s),h.y(s),dx,dy,f);
427  }
428  if(ifvtx==0 && mBeamLine==1){
429  double z,dca,eta,phi,opt;
430  getZvtxAndDcaFromHelix(par,dca,z,eta,phi,opt);
431  f += dca*dca/mErrVTX/mErrVTX;
432  ndf++;
433  }
434  ndf-=5;
435  if(ndf<=0) printf("ERROR fitHelix ndf=%d\n",ndf);
436  f/=double(ndf);
437  if(mDebug>4) printf("fitHelix x0=%12.8f y0=%12.8f c=%12.8f dip=%12.8f phase=%12.8f f=%12.8f\n",
438  par[0],par[1],par[2],par[3],par[4],f);
439 }
440 
441 //Get residuals with helix fit
442 double residHelix(double *par){
443  double f=0, chi2=0;
444  StPhysicalHelix h = getStPhysicalHelix(par);
445  static StThreeVectorD n(0,0,1);
446  int ifvtx=0, ifvtx2=0;
447  int ndf=0,ndf2=0;
448  for(int i=0; i<mHitAlg.nhit; i++){
449  StThreeVectorD z(0,0,mHitAlg.z[i]);
450  double s=h.pathLength(z,n);
451  double x=h.x(s);
452  double y=h.y(s);
453  double dx = mHitAlg.x[i] - x;
454  double dy = mHitAlg.y[i] - y;
455  int disc=mHitAlg.det[i]/4;
456  double ff=(dx*dx+dy*dy)/(mHitAlg.ex[i]*mHitAlg.ex[i]+mHitAlg.ey[i]*mHitAlg.ey[i]);
457  if(mResidMaskDisc & (1<<disc)) {
458  if(disc==6) ifvtx=1;
459  f += ff;
460  ndf+=2;
461  }
462  if(mHitMaskDisc & (1<<disc)) {
463  if(disc==6) ifvtx2=1;
464  chi2 += ff;
465  ndf2+=2;
466  }
467  if(mDebug>3) printf("residHelix %3d dx=%12.8f dy=%12.8f f=%12.8f\n",i,dx,dy,f);
468  if(mFillHist>0) {
469  double r=sqrt(x*x+y*y);
470  double phi=atan2(y,x);
471  double hx=mHitAlg.x[i];
472  double hy=mHitAlg.y[i];
473  double hr=mHitAlg.r[i];
474  double hp=mHitAlg.p[i];
475  double dr = hr-r;
476  double dp = hp-phi;
477  //if(disc>8) printf("EEMCDR n=%2d used=%d %8.4f ->%8.4f\n",mHitAlg.nhitEemc,mHitAlg.used,mHitAlg.dr[i],dr);
478  while(dp>PI) dp-=2*PI;
479  while(dp<-PI) dp+=2*PI;
480  mHitAlg.dx[i]=dx;
481  mHitAlg.dy[i]=dy;
482  mHitAlg.dr[i]=dr;
483  mHitAlg.dp[i]=dp;
484  int quad=mHitAlg.det[i]%4;
485  if(disc<kFgtNumDiscs){
486  while(phi>StFgtGeom::phiQuadXaxis(2)) phi-=2*PI;
487  while(phi<StFgtGeom::phiQuadXaxis(1)) phi+=2*PI;
488  }
489  fillHist(disc,quad,dx,dy,dr,dp,x,y,r,phi);
490  }
491  }
492  if(ifvtx+ifvtx2<2 && mBeamLine==1){
493  double z,dca,eta,phi,opt;
494  getZvtxAndDcaFromHelix(par,dca,z,eta,phi,opt);
495  if(ifvtx==0) {f += dca*dca/mErrVTX/mErrVTX; ndf++; }
496  if(ifvtx2==0) {chi2 += dca*dca/mErrVTX/mErrVTX; ndf2++;}
497  }
498  ndf-=5;
499  if(ndf<=0) printf("ERROR resudHelix ndf=%d\n",ndf);
500  f/=double(ndf);
501  if(mFillHist>0) {
502  ndf2-=5;
503  if(ndf2<=0) {printf("ERROR residHelix ndf2=%d\n",ndf2);}
504  else {chi2/=double(ndf2);}
505  double z,dca,eta,phi,opt;
506  getZvtxAndDcaFromHelix(par,dca,z,eta,phi,opt);
507  if(mDebug>0) printf("UPDATE Dca=%8.4f->%8.4f trkz=%8.4f->%8.4f chi2=%8.4f->%8.4f ndf=%d\n",mHitAlg.dca,dca,mHitAlg.trkz,z,mHitAlg.chi2,chi2,ndf2);
508  mHitAlg.dca=dca;
509  mHitAlg.dz=mHitAlg.dz-mHitAlg.trkz+z;
510  mHitAlg.trkz=z;
511  mHitAlg.chi2=chi2;
512  mHitAlg.eta=eta;
513  mHitAlg.phi=phi;
514  mHitAlg.opt=opt;
515  }
516  if(mFillHist>1 && mEventDisp) eventDisplay(par);
517  if(mDebug>2)
518  printf("residHelix x0=%12.8f y0=%12.8f c=%12.8f dip=%12.8f phase=%12.8f f=%12.8f\n",
519  par[0],par[1],par[2],par[3],par[4],f);
520  return f;
521 }
522 
523 //a track fit to helix
524 void fitTrackHelix(fgtAlignment_st* algpar, int itrk, double *par){
525  int flag;
526  static double arg[10];
527  static TMinuit *m = 0;
528  if(m==0){
529  m=new TMinuit(5);
530  m->SetFCN(fitHelix);
531  arg[0] =-1; m->mnexcm("SET PRI", arg, 1,flag);
532  arg[0] =-1; m->mnexcm("SET NOWarning", arg, 0,flag);
533  arg[0] = 1; m->mnexcm("SET ERR", arg, 1,flag);
534  arg[0] = 500; arg[1] = 1.;
535  }
536  //make "aligned" hits in mHitAlg
537  getAlign(itrk,algpar);
538  //first guess of helix fit parameters and set up
539  double dx = mHitAlg.x[1]-mHitAlg.x[0];
540  double dy = mHitAlg.y[1]-mHitAlg.y[0];
541  double dz = mHitAlg.z[1]-mHitAlg.z[0];
542  double dr = sqrt(dx*dx+dy*dy);
543  double dip = atan2(dz,dr);
544  double slp = atan2(dy,dx);
545  double x0= mHitAlg.x[0]-dx/dz*mHitAlg.z[0];
546  double y0= mHitAlg.y[0]-dy/dz*mHitAlg.z[0];
547  m->mnparm(0,"x0" ,x0 ,0.1, 0, 0,flag);
548  m->mnparm(1,"y0" ,y0 ,0.1, 0, 0,flag);
549  m->mnparm(2,"curv",0.0001 ,1.0, 0, 0,flag);
550  m->mnparm(3,"dip" ,dip ,0.1, 0, 0,flag);
551  m->mnparm(4,"slp" ,slp ,0.1, 0, 0,flag);
552  //do helix fit
553  m->mnexcm("MIGRAD", arg ,2, flag);
554  double r,e;
555  for(int j=0; j<5; j++){ m->GetParameter(j,r,e); par[j]=r;}
556 }
557 
558 //Minimize residuals with helix fit, called from TMinuit
559 void funcHelix(Int_t &npar, Double_t* gin, Double_t &f, Double_t *par, Int_t iflag){
560  if(mDebug>3) printf("funcHelix mQuad=%d mNtrk=%d\n",mQuad,mNtrk[mQuad]);
561  f=0;
562  fgtAlignment_st* algpar= (fgtAlignment_st*)par;
563  double p[5];
564  for(int itrk=0; itrk<mNtrk[mQuad]; itrk++){
565  if(mHit[mQuad][itrk].used==0) continue;
566  fitTrackHelix(algpar,itrk,p);
567  f+=residHelix(p); //add up residuals
568  if(mFillHist>0) {
569  memcpy(&mHit[mQuad][itrk],&mHitAlg,sizeof(mHitAlg));
570  //printf("mHitAlg.z[0]=%8.2f\n",mHitAlg.z[0]);
571  //printf("mHit[mQuad][itrk].z[0]=%8.2f\n",mHit[mQuad][itrk].z[0]);
572  }
573  }
574  if(mDebug>3) printf("funcHelix f=%12.6f xoff=%12.8f yoff=%12.8f\n",f,algpar->xoff[8],algpar->xoff[8]);
575  static int ncall[4]={0,0,0,0};
576  ncall[mQuad]++;
577  if(ncall[mQuad]%100==0) printf("funcHelix quad=%d ncall=%d f=%12.6f\n",mQuad,ncall[mQuad],f);
578 }
579 
580 //Straight line fit, called from TMinuit
581 void fitLine(Int_t &npar, Double_t* gin, Double_t &f, Double_t *par, Int_t iflag){
582  if(mDebug>4) printf("fitLine nhit=%d\n",mHitAlg.nhit);
583  f=0;
584  int ifvtx=0;
585  int ndf=0;
586  for(int i=0; i<mHitAlg.nhit; i++){
587  if(mHitAlg.use[i]==false) continue;
588  if(mHitAlg.det[i]/4==6) ifvtx=1;
589  double dx = par[0] + par[1] * mHitAlg.z[i] - mHitAlg.x[i];
590  double dy = par[2] + par[3] * mHitAlg.z[i] - mHitAlg.y[i];
591  f += (dx*dx+dy*dy)/(mHitAlg.ex[i]*mHitAlg.ex[i]+mHitAlg.ey[i]*mHitAlg.ey[i]);
592  ndf+=2;
593  if(mDebug>5) printf("dx=%8.4f dy=%8.4f dr=%8.4f\n",dx,dy,f);
594  }
595  if(ifvtx==0 && mBeamLine==1){
596  double z,dca;
597  getZvtxAndDca(par,dca,z);
598  f += dca*dca/mErrVTX/mErrVTX;
599  ndf++;
600  }
601  ndf-=4;
602  if(ndf<=0) {printf("ERROR fitLine ndf=%d\n",ndf);}
603  else {f/=double(ndf);}
604  if(mDebug>4) printf("fitLine x0=%12.8f x1=%12.8f y0=%12.8f y1=%12.8f f=%12.8f\n",par[0],par[1],par[2],par[3],f);
605 }
606 
607 //Get residuals with straight line fit
608 double residLine(double *par){
609  //printf("residLine mQuad=%d mNtrk=%d mFillHist=%d\n",mQuad,mNtrk[mQuad],mFillHist);
610  double f=0, chi2=0;
611  int ifvtx=0, ifvtx2=0;
612  int ndf=0, ndf2=0;
613  for(int i=0; i<mHitAlg.nhit; i++){
614  double x=par[0] + par[1] * mHitAlg.z[i];
615  double y=par[2] + par[3] * mHitAlg.z[i];
616  double dx = mHitAlg.x[i]-x;
617  double dy = mHitAlg.y[i]-y;
618  int disc=mHitAlg.det[i]/4;
619  if(mResidMaskDisc & (1<<disc)) {
620  f += (dx*dx+dy*dy)/(mHitAlg.ex[i]*mHitAlg.ex[i]+mHitAlg.ey[i]*mHitAlg.ey[i]);
621  ndf+=2;
622  }
623  if(mHitMaskDisc & (1<<disc)) {
624  chi2 += (dx*dx+dy*dy)/(mHitAlg.ex[i]*mHitAlg.ex[i]+mHitAlg.ey[i]*mHitAlg.ey[i]);
625  ndf2+=2;
626  }
627  if(mDebug>5) printf("residLine %3d dx=%12.8f dy=%12.8f dr=%8.2f\n",i,dx,dy,dx*dx+dy*dy);
628  if(mFillHist>0) {
629  double r=sqrt(x*x+y*y);
630  double phi=atan2(y,x);
631  double hx=mHitAlg.x[i];
632  double hy=mHitAlg.y[i];
633  double hr=mHitAlg.r[i];
634  double hp=mHitAlg.p[i];
635  double dr = hr-r;
636  double dp = hp-phi;
637  while(dp>PI) dp-=2*PI;
638  while(dp<-PI) dp+=2*PI;
639  mHitAlg.dx[i]=dx;
640  mHitAlg.dy[i]=dy;
641  mHitAlg.dr[i]=dr;
642  mHitAlg.dp[i]=dp;
643  int quad=mHitAlg.det[i]%4;
644  if(disc<kFgtNumDiscs){
645  while(phi>StFgtGeom::phiQuadXaxis(2)) phi-=2*PI;
646  while(phi<StFgtGeom::phiQuadXaxis(1)) phi+=2*PI;
647  }
648  fillHist(disc,quad,dx,dy,dr,dp,x,y,r,phi);
649  //if(disc==6 && fabs(dx)<0.01 && fabs(dy)<0.01){
650  // printf("VTX %3d q=%1d d=%1d trk=%6.4f %6.4f %6.4f %6.4f fgt=%6.4f %6.4f %6.4f %6.4f xyrp=%6.4f %6.4f %6.4f %6.4f dxyrp=%6.4f %6.4f %6.4f %6.4f\n",
651  // i,quad,disc,par[0],par[1],par[2],par[3],x,y,r,phi,hx,hy,hr,hp,dx,dy,dr,dp);
652  //}
653  //if(disc>6)
654  //printf("EEMC2 %3d q=%1d d=%1d trk=%6.2f %6.2f %6.2f %6.2f fgt=%6.1f %6.1f %6.1f %6.2f xyrp=%6.1f %6.1f %6.1f %6.2f dxyrp=%6.1f %6.1f %6.1f %6.2f\n",
655  // i,quad,disc,par[0],par[1],par[2],par[3],x,y,r,phi,hx,hy,hr,hp,dx,dy,dr,dp);
656  }
657  }
658  if(ifvtx+ifvtx2<2 && mBeamLine==1){
659  double z,dca;
660  getZvtxAndDca(par,dca,z);
661  if(ifvtx==0) {f += dca*dca/mErrVTX/mErrVTX; ndf++; }
662  if(ifvtx2==0) {chi2 += dca*dca/mErrVTX/mErrVTX; ndf2++;}
663  }
664  ndf-=4;
665  if(ndf<=0) printf("ERROR residLine ndf=%d\n",ndf);
666  f/=double(ndf);
667  if(mFillHist>0) {
668  ndf2-=4;
669  if(ndf2<=0) printf("ERROR residLine ndf2=%d\n",ndf2);
670  chi2/=double(ndf2);
671  //update DCA and Trkz
672  double dca,z;
673  getZvtxAndDca(par,dca,z);
674  if(mDebug>0) printf("UPDATE Dca=%8.4f->%8.4f trkz=%8.4f->%8.4f chi2=%8.4f->%8.4f ndf=%d\n",mHitAlg.dca,dca,mHitAlg.trkz,z,mHitAlg.chi2,chi2,ndf2);
675  mHitAlg.dca=dca;
676  mHitAlg.dz=mHitAlg.dz-mHitAlg.trkz+z;
677  mHitAlg.trkz=z;
678  mHitAlg.chi2=chi2;
679  double theta=atan2(sqrt(par[1]*par[1]+par[3]*par[3]),1.0);
680  mHitAlg.eta=-log(tan(theta/2.0)) ;
681  mHitAlg.phi=atan2(par[3],par[1]);
682  mHitAlg.opt=0.0;
683  }
684  if(mFillHist>1 && mEventDisp) eventDisplay(par);
685  if(mDebug>4) printf("residLine x0=%12.8f x1=%12.8f y0=%12.8f y1=%12.8f f=%12.8f\n",par[0],par[1],par[2],par[3],f);
686  return f;
687 }
688 
689 //a track fit to line
690 void fitTrackLine(fgtAlignment_st* algpar, int itrk, double *par){
691  int flag;
692  static double arg[10];
693  static TMinuit *m = 0;
694  if(m==0){
695  m=new TMinuit(4);
696  m->SetFCN(fitLine);
697  arg[0] =-1; m->mnexcm("SET PRI", arg, 1,flag);
698  arg[0] =-1; m->mnexcm("SET NOWarning", arg, 0,flag);
699  arg[0] = 1; m->mnexcm("SET ERR", arg, 1,flag);
700  arg[0] = 500; arg[1] = 1.;
701  }
702  //make "aligned" hits in mHitAlg
703  getAlign(itrk,algpar);
704  //first guess of helix fit parameters and set up
705  double x1 = (mHitAlg.x[1]-mHitAlg.x[0])/(mHitAlg.z[1]-mHitAlg.z[0]);
706  double y1 = (mHitAlg.y[1]-mHitAlg.y[0])/(mHitAlg.z[1]-mHitAlg.z[0]);
707  double x0 = mHitAlg.x[0] - x1*mHitAlg.z[0];
708  double y0 = mHitAlg.y[0] - y1*mHitAlg.z[0];
709  m->mnparm(0,"x0",x0,0.1,0,0,flag);
710  m->mnparm(1,"x1",x1,0.1,0,0,flag);
711  m->mnparm(2,"y0",y0,0.1,0,0,flag);
712  m->mnparm(3,"y1",y1,0.1,0,0,flag);
713  //do Line fit
714  m->mnexcm("MIGRAD", arg ,2, flag);
715  double r,e;
716  for(int j=0; j<4; j++){ m->GetParameter(j,r,e); par[j]=r;}
717 }
718 
719 //Minimize residuals with line fit, called from TMinuit
720 void funcLine(Int_t &npar, Double_t* gin, Double_t &f, Double_t *par, Int_t iflag){
721  //printf("funcLine mQuad=%d mNtrk=%d mFillHist=%d\n",mQuad,mNtrk[mQuad],mFillHist);
722  if(mDebug>3) printf("funcLine mQuad=%d mNtrk=%d\n",mQuad,mNtrk[mQuad]);
723  f=0;
724  fgtAlignment_st* algpar= (fgtAlignment_st*)par;
725  double p[4];
726  for(int itrk=0; itrk<mNtrk[mQuad]; itrk++){
727  if(mHit[mQuad][itrk].used==0) continue;
728  fitTrackLine(algpar,itrk,p); //track fit
729  f+=residLine(p); //add up residuals
730  if(mFillHist>0) {
731  memcpy(&mHit[mQuad][itrk],&mHitAlg,sizeof(mHitAlg));
732  //printf("mHitAlg.z[0]=%8.2f\n",mHitAlg.z[0]);
733  //printf("mHit[mQuad][itrk].z[0]=%8.2f\n",mHit[mQuad][itrk].z[0]);
734  }
735  }
736  static int ncall[4]={0,0,0,0};
737  ncall[mQuad]++;
738  if(ncall[mQuad]%100==0) printf("funcLine quad=%d ncall=%d f=%12.6f\n",mQuad,ncall[mQuad],f);
739 }
740 
741 ClassImp(StFgtAlignmentMaker);
742 
743 StFgtAlignmentMaker::StFgtAlignmentMaker(const Char_t *name) : StMaker(name),mEventCounter(0),
744  mErrFgt(0.02), mErrVtx(0.02), mErrVtxZ(1.0),
745  mErrTpcI(0.6), mErrTpcO(0.12),mErrTpcZ(0.12),mErrPpt(0.1),mErrEmc(0.3),
746  mFakeNtrk(2000),mFakeEmin(40),mFakeEmax(40),mFakeEtamin(1.6),mFakeEtamax(1.6),
747  mFakePhimin(0),mFakePhimax(0),mFakeVtxSig(0),
748  mDataSource(0),
749  mOutTreeFile(0),mInTreeFile(0),mReadParFile(0),
750  mRunNumber(0),mSeqNumber(0),mDay(0),mNStep(0){
751  memset(mDzCut ,0,sizeof(mDzCut ));
752  memset(mDcaCut ,0,sizeof(mDcaCut ));
753  memset(mFgtRCut,0,sizeof(mFgtRCut));
754  memset(mFgtPCut,0,sizeof(mFgtPCut));
755  memset(mTpcRCut,0,sizeof(mTpcRCut));
756  memset(mTpcPCut,0,sizeof(mTpcPCut));
757  memset(mEmcRCut,0,sizeof(mEmcRCut));
758  memset(mEmcPCut,0,sizeof(mEmcPCut));
759 }
760 
761 Int_t StFgtAlignmentMaker::Init(){
762  memset(mNtrk,0,sizeof(mNtrk));
763  memset(mHit,0,sizeof(mHit));
764  bookHist();
765  return kStOK;
766 }
767 
768 Int_t StFgtAlignmentMaker::InitRun(Int_t runnum){
769  LOG_INFO << "StFgtAlignmentMaker::InitRun for " << runnum << endm;
770  StFgtDbMaker *fgtDbMkr = static_cast< StFgtDbMaker* >( GetMakerInheritsFrom( "StFgtDbMaker" ) );
771  //StFgtDbMaker *fgtDbMkr = static_cast<StFgtDbMaker * >( GetMaker("fgtDb"));
772  if( !fgtDbMkr ){
773  LOG_FATAL << "StFgtDb not provided and error finding StFgtDbMaker" << endm;
774  return kStFatal;
775  } else {
776  mDb = fgtDbMkr->getDbTables();
777  if( !mDb ){
778  LOG_FATAL << "StFgtDb not provided and error retrieving pointer from StFgtDbMaker '"
779  << fgtDbMkr->GetName() << endm;
780  return kStFatal;
781  }
782  }
783  if(mReadParFile==0){
784  cout << "mDb="<<mDb<<endl;
785  orig_algpar=mDb->getAlignment();
786  }else{
787  readPar(orig_algpar);
788  }
789  return kStOK;
790 }
791 
792 static const int mMaxStep=100;
793 
794 void StFgtAlignmentMaker::setStep(int discmask,int quadmask, int parmask, int hitmask_disc, int residmask,
795  int trackType, int minHit, int minFgtHit, int minVtx, int minTpcHit, int minPromptHit, int minEemcHit,
796  float dzcut, float dcacut, float fgtrcut, float fgtpcut, float tpcrcut, float tpcpcut, float emcrcut, float emcpcut){
797  if(mNStep>=mMaxStep) {printf("Reached MaxStep\n"); return; }
798  if(mNStep==0) { printf("Step0 is making before histo, and masks are set to 0\n"); discmask=0; quadmask=0; parmask=0; }
799  mDiscMask[mNStep]=discmask;
800  mQuadMask[mNStep]=quadmask;
801  mParMask[mNStep]=parmask;
802  mHitMask[mNStep]=hitmask_disc;
803  mResidMask[mNStep]=residmask;
804  mTrackType[mNStep]=trackType;
805  mMinHit[mNStep]=minHit;
806  mMinFgtHit[mNStep]=minFgtHit;
807  mMinVtx[mNStep]=minVtx;
808  mMinTpcHit[mNStep]=minTpcHit;
809  mMinPromptHit[mNStep]=minPromptHit;
810  mMinEemcHit[mNStep]=minEemcHit;
811  mDzCut[mNStep]=dzcut;
812  mDcaCut[mNStep]=dcacut;
813  mFgtRCut[mNStep]=fgtrcut;
814  mFgtPCut[mNStep]=fgtpcut;
815  mTpcRCut[mNStep]=tpcrcut;
816  mTpcPCut[mNStep]=tpcpcut;
817  mEmcRCut[mNStep]=emcrcut;
818  mEmcPCut[mNStep]=emcpcut;
819  printf("Adding step=%d with discMask=%x quadMask=%x parMask=%x hitMask=%x trkType=%d minHit=%d minFgt=%d minTpc=%d minPrompt=%d minEemc=%d\n",
820  mNStep,discmask,quadmask,parmask,hitmask_disc,trackType,minHit,minFgtHit,minTpcHit,minPromptHit,minEemcHit);
821  printf(" Cuts dz=%8.3f dca=%8.3f fgtdr=%8.3f fgtdp=%8.3f tpcdr=%8.3f tpcdp=%8.3f emcdr=%8.3f emcdp=%8.3f\n",
822  mDzCut[mNStep],mDcaCut[mNStep],mFgtRCut[mNStep],mFgtPCut[mNStep],mTpcRCut[mNStep],mTpcPCut[mNStep],mEmcRCut[mNStep],mEmcPCut[mNStep]);
823  mNStep++;
824 }
825 
826 
827 Int_t StFgtAlignmentMaker::Make() {
828  orig_algpar=mDb->getAlignment();
829  mErrFGT=mErrFgt;
830  mErrVTX=mErrVtx;
831  mErrVTXZ=mErrVtxZ;
832  mErrTPCI=mErrTpcI;
833  mErrTPCO=mErrTpcO;
834  mErrTPCZ=mErrTpcZ;
835  mErrPPT=mErrPpt;
836  mErrEMC=mErrEmc;
837  if(mDataSource==0){
838  readFromStEvent();
839  }else if(mDataSource==1){
840  readFromStEventGlobal();
841  }else if(mDataSource==2){
842  readFromStraightTrackMaker();
843  DispFromStraightTrackMaker();
844  }else if(mDataSource==5){
845  calcETBalanceFromStEvent();
846  readFromStraightTrackAndStEvent();
847  }else if(mDataSource==6){
848  readFromStraightTrackAndMudst();
849  }
850  return kStOK;
851 }
852 
853 Int_t StFgtAlignmentMaker::Finish() {
854  gMessMgr->Info() << "StFgtAlignmentMaker::Finish()" << endm;
855 
856  if(mDataSource==4) {fakeData();}
857  else if(mDataSource==3) {readFromTree();}
858  overWriteError();
859  fgtAlignment_st result;
860  memcpy(&result,orig_algpar,sizeof(fgtAlignment_st));
861 
862  cout << "Doing Alignment with Number of steps = "<<mNStep<<endl;
863  for(int s=0; s<mNStep; s++){
864  mStep=s;
865  if(mDiscMask[s]+mParMask[s]>0) {mFillHist=0;}
866  else {mFillHist=1; resetHist(); printf("Saving hist for step=%d\n",mStep);}
867  for(int quad=0; quad<kFgtNumQuads; quad++){
868  mQuad=quad;
869  int quadmask= 1<<quad;
870  if( (mFillHist>0) || (quadmask & mQuadMask[s]) ){
871  cout << Form("Doing alignment for quad=%1d with Ntrk=%4d",quad,mNtrk[quad])<<endl;
872  if(mNtrk[mQuad]>0){
873  printf("Quad=%d Step=%d with discMask=%x quadMask=%x parMask=%x hitMask=%x residMask=%x trkType=%d minHit=%d minFgt=%d vtx=%d minTpc=%d minPrompt=%d minEEmc=%d\n",
874  quad,s,mDiscMask[s],quadmask,mParMask[s],mHitMask[s],mResidMask[s],
875  mTrackType[s],mMinHit[s],mMinFgtHit[s],mMinVtx[s],mMinTpcHit[s],mMinPromptHit[s],mMinEemcHit[s]);
876  printf(" Cuts dz=%8.3f dca=%8.3f fgtdr=%8.3f fgtdp=%8.3f tpcdr=%8.3f tpcdp=%8.3f emcdr=%8.3f emcdp=%8.3f\n",
877  mDzCut[s],mDcaCut[s],mFgtRCut[s],mFgtPCut[s],mTpcRCut[s],mTpcPCut[s],mEmcRCut[s],mEmcPCut[s]);
878  doAlignment(&result,mDiscMask[s],quadmask,mParMask[s],mHitMask[s],mResidMask[s],
879  mTrackType[s],mMinHit[s],mMinFgtHit[s],mMinVtx[s],mMinTpcHit[s],mMinPromptHit[s],mMinEemcHit[s],&result);
880  }
881  }
882  }
883  if(mFillHist>0){
884  saveHist();
885  if(mOutTreeFile) writeTree();
886  }
887  }
888  writePar(&result);
889  return kStOK;
890 }
891 
892 void StFgtAlignmentMaker::fakeData() {
893  mFgtInputRPhi=0;
894  orig_algpar=new fgtAlignment_st;
895  memset(orig_algpar,0,sizeof(fgtAlignment_st));
896  int quad=0;
897  mNtrk[quad]=0;
898  memset(mHit,0,sizeof(mHit));
899  int opt=1;
900  if(opt==0){
901  for(int itrk=0; itrk<1; itrk++){
902  mHit[quad][itrk].nhit=0;
903  for(int d=0; d<10; d++){
904  double z;
905  mHit[quad][itrk].nhit++;
906  if (d==0) {mHit[quad][itrk].det[d]=24+quad; z=0;}
907  else if(d>6) {mHit[quad][itrk].det[d]=28+quad; z=200.0+d;}
908  else{
909  mHit[quad][itrk].det[d]=(d-1)*4+quad;
910  z=StFgtGeom::getDiscZ(d-1);
911  }
912  mHit[quad][itrk].x[d]=20.0/67.399*(z-10); mHit[quad][itrk].ex[d]=0.1;
913  mHit[quad][itrk].y[d]=0.0; mHit[quad][itrk].ey[d]=0.1;
914  mHit[quad][itrk].z[d]=z; mHit[quad][itrk].ez[d]=0.1;
915  if(d==3) mHit[quad][itrk].x[d]+=0.2; // added mis-alignment
916  if(mDebug>0){
917  cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
918  itrk,d,quad,mHit[quad][itrk].det[d],
919  mHit[quad][itrk].x[d],mHit[quad][itrk].y[d],mHit[quad][itrk].z[d],
920  mHit[quad][itrk].ex[d],mHit[quad][itrk].ey[d],mHit[quad][itrk].ez[d])
921  <<endl;
922  }
923  }
924  mNtrk[quad]++;
925  }
926  }else{
927  //fake track inputs
928  mQuad=0;
929  TRandom rand;
930  int mFgt=1,mVtx=1,mTpc=0,mPpt=1,mEmc=1;
931  const double B = 0.5*tesla;
932  StThreeVectorD zvec(0,0,1);
933  //TPC and EEMC z positions
934  static const int NTPC=45;
935  double tpcr[NTPC];
936  for(int ipad= 0;ipad< 8; ipad++){tpcr[ipad]=60.0 + 4.8*ipad;}
937  for(int ipad= 8;ipad<13; ipad++){tpcr[ipad]=60.0 + 4.8*7 + 5.2*(ipad-8);}
938  for(int ipad=13;ipad<45; ipad++){tpcr[ipad]=127.950 + 2.0*(ipad-13);}
939  double zppt=209;
940  double zemc=280;
941 
942  for(int itrk=0; itrk<mFakeNtrk; itrk++){
943  double ene=mFakeEmin;
944  double eta=mFakeEtamin;
945  double phi=mFakePhimin;
946  StThreeVectorD v(0,0,0);
947  if(mFakeEmax>mFakeEmin) ene=rand.Uniform(mFakeEmin,mFakeEmax);
948  if(mFakeEtamax>mFakeEtamin) eta=rand.Uniform(mFakeEtamin,mFakeEtamax);
949  if(mFakePhimax>mFakePhimin) phi=rand.Uniform(mFakePhimin,mFakePhimax);
950  if(mFakeVtxSig>0) {
951  double zz=-999;
952  while (zz>60 || zz<-120) {zz=rand.Gaus(0,mFakeVtxSig);}
953  v.setZ(zz);
954  }
955  //Making helix
956  double theta=2*atan(exp(-eta));
957  double pt=ene*sin(theta);
958  double pz=ene*cos(theta);
959  double px=ene*sin(theta)*cos(phi);
960  double py=ene*sin(theta)*sin(phi);
961  StThreeVectorD p(px, py, pz);
962  StPhysicalHelixD pos(p,v,B,1);
963  StPhysicalHelixD neg(p,v,B,-1);
964  if(itrk<3){
965  printf("FAKE pt=%8.3f eta=%8.3f phi=%8.3f ene=%8.3f\n",pt,eta,phi,ene);
966  printf("FAKE pxyz=%8.3f %8.3f %8.3f\n",px,py,pz);
967  printf("FAKE vxyz=%8.3f %8.3f %8.3f\n",v.x(),v.y(),v.z());
968  //print rough trajectory
969  for (double z=0; z<=300; z+=20) {
970  StThreeVectorD zplane(0,0,z);
971  double sp=pos.pathLength(zplane,zvec);
972  double sn=pos.pathLength(zplane,zvec);
973  StThreeVectorD pxyz = pos.at(sp);
974  StThreeVectorD nxyz = neg.at(sn);
975  printf("FAKE z=%8.3f pos=%8.3f %8.3f neg=%8.3f %8.3f\n",z,pxyz.x(),pxyz.y(),nxyz.x(),nxyz.y());
976  }
977  }
978 
979  //make hits
980  StPhysicalHelixD trk;
981  if(itrk%2==0) {trk=pos;}
982  else {trk=neg;}
983  int nhit=0;
984  if(mVtx){
985  StThreeVectorD zplane(0,0,v.z());
986  double s=trk.pathLength(zplane,zvec);
987  StThreeVectorD xyz = trk.at(s);
988  mHit[mQuad][itrk].x[nhit]=xyz.x(); mHit[mQuad][itrk].y[nhit]=xyz.y(); mHit[mQuad][itrk].z[nhit]=xyz.z();
989  mHit[mQuad][itrk].ex[nhit]=mErrVtx; mHit[mQuad][itrk].ey[nhit]=mErrVtx; mHit[mQuad][itrk].ez[nhit]=mErrVtxZ;
990  mHit[mQuad][itrk].det[nhit]=4*6+mQuad;
991  nhit++;
992  if(itrk<5) printf("VTX %8.3f %8.3f %8.3f %8.4f %8.4f\n",xyz.x(),xyz.y(),xyz.z(),mErrVtx,mErrVtxZ);
993  }
994  mHit[mQuad][itrk].dz=0;
995  mHit[mQuad][itrk].trkz=v.z();
996  if(mFgt){
997  int nfgt=0;
998  for(int idisc=0; idisc<6; idisc++){
999  if(idisc==3) continue;
1000  StThreeVectorD zplane(0,0,StFgtGeom::getDiscZ(idisc));
1001  double s=trk.pathLength(zplane,zvec);
1002  StThreeVectorD xyz = trk.at(s);
1003  if(xyz.perp()>10.0 && xyz.perp()<38.3){
1004  mHit[mQuad][itrk].x[nhit]=xyz.x(); mHit[mQuad][itrk].y[nhit]=xyz.y(); mHit[mQuad][itrk].z[nhit]=xyz.z();
1005  mHit[mQuad][itrk].ex[nhit]=mErrFgt; mHit[mQuad][itrk].ey[nhit]=mErrFgt; mHit[mQuad][itrk].ez[nhit]=0.0;
1006  mHit[mQuad][itrk].det[nhit]=4*idisc+mQuad;
1007  nhit++;
1008  nfgt++;
1009  if(itrk<5) printf("FGT%1d %8.3f %8.3f %8.3f %8.4f\n",idisc+1,xyz.x(),xyz.y(),xyz.z(),mErrFgt);
1010  }
1011  }
1012  if(nfgt<3) {continue;}
1013  }
1014  if(mTpc){
1015  for(int ipad=0; ipad<NTPC; ipad++){
1016  double r=tpcr[ipad];
1017  pair <double,double> ss=trk.pathLength(r);
1018  double s=ss.first;
1019  if(s<0) s=ss.second;
1020  StThreeVectorD xyz = trk.at(s);
1021  //double z=r/tan(theta);
1022  if(xyz.z()<198){
1023  mHit[mQuad][itrk].x[nhit]=xyz.x(); mHit[mQuad][itrk].y[nhit]=xyz.y(); mHit[mQuad][itrk].z[nhit]=xyz.z();
1024  if(ipad<13) {mHit[mQuad][itrk].ex[nhit]=mErrTpcI; mHit[mQuad][itrk].ey[nhit]=mErrTpcI;}
1025  else {mHit[mQuad][itrk].ex[nhit]=mErrTpcO; mHit[mQuad][itrk].ey[nhit]=mErrTpcO;}
1026  mHit[mQuad][itrk].ez[nhit]=mErrTpcZ;
1027  mHit[mQuad][itrk].det[nhit]=4*7+mQuad;
1028  nhit++;
1029  if(itrk<5) printf("TPC %8.3f %8.3f %8.3f %8.4f %8.4f\n",xyz.x(),xyz.y(),xyz.z(),mHit[mQuad][itrk].ex[nhit],mErrTpcZ);
1030  }
1031  }
1032  }
1033  if(mPpt){
1034  StThreeVectorD zplane(0,0,zppt);
1035  double s=trk.pathLength(zplane,zvec);
1036  StThreeVectorD xyz = trk.at(s);
1037  if(xyz.perp()>tpcr[0] && xyz.perp()<tpcr[NTPC-1]){
1038  mHit[mQuad][itrk].x[nhit]=xyz.x(); mHit[mQuad][itrk].y[nhit]=xyz.y(); mHit[mQuad][itrk].z[nhit]=xyz.z();
1039  mHit[mQuad][itrk].ex[nhit]=mErrPpt; mHit[mQuad][itrk].ey[nhit]=mErrPpt; mHit[mQuad][itrk].ez[nhit]=0.0;
1040  mHit[mQuad][itrk].det[nhit]=4*8+mQuad;
1041  nhit++;
1042  if(itrk<5) printf("PPT %8.3f %8.3f %8.3f %8.4f\n",xyz.x(),xyz.y(),xyz.z(),mErrPpt);
1043  }
1044  }
1045  if(mEmc){
1046  StThreeVectorD zplane(0,0,zemc);
1047  double s=trk.pathLength(zplane,zvec);
1048  StThreeVectorD xyz = trk.at(s);
1049  if(xyz.perp()>tpcr[0] && xyz.perp()<tpcr[NTPC-1]){
1050  mHit[mQuad][itrk].x[nhit]=xyz.x(); mHit[mQuad][itrk].y[nhit]=xyz.y(); mHit[mQuad][itrk].z[nhit]=xyz.z();
1051  mHit[mQuad][itrk].ex[nhit]=mErrEmc; mHit[mQuad][itrk].ey[nhit]=mErrEmc; mHit[mQuad][itrk].ez[nhit]=0.0;
1052  mHit[mQuad][itrk].det[nhit]=4*11+mQuad;
1053  mHit[mQuad][itrk].tw[nhit]=float(ene);
1054  mHit[mQuad][itrk].nrl[nhit]=0;
1055  mHit[mQuad][itrk].su[nhit]=100.0;
1056  mHit[mQuad][itrk].sv[nhit]=100.0;
1057  nhit++;
1058  if(itrk<5) printf("EMC %8.3f %8.3f %8.3f %8.4f\n",xyz.x(),xyz.y(),xyz.z(),mErrEmc);
1059  }
1060  }
1061  mHit[mQuad][itrk].nhit=nhit;
1062  mHit[mQuad][itrk].dca=0.0;
1063  mHit[mQuad][itrk].chi2=0.0;
1064  mHit[mQuad][itrk].riso=1.0;
1065  //adding gauusian error for each track
1066  for(int ii=0; ii<nhit; ii++){
1067  mHit[mQuad][itrk].x[ii]+= mHit[mQuad][itrk].ex[ii]*rand.Gaus();
1068  mHit[mQuad][itrk].y[ii]+= mHit[mQuad][itrk].ey[ii]*rand.Gaus();
1069  mHit[mQuad][itrk].z[ii]+= mHit[mQuad][itrk].ez[ii]*rand.Gaus();
1070  mHit[mQuad][itrk].r[ii]=sqrt(mHit[mQuad][itrk].x[ii]*mHit[mQuad][itrk].x[ii]+mHit[mQuad][itrk].y[ii]*mHit[mQuad][itrk].y[ii]);
1071  mHit[mQuad][itrk].p[ii]=atan2(mHit[mQuad][itrk].y[ii],mHit[mQuad][itrk].x[ii]);
1072  }
1073  }
1074  mNtrk[mQuad]=mFakeNtrk;
1075  }
1076 }
1077 
1078 void StFgtAlignmentMaker::readFromStraightTrackMaker(){
1079  mFgtInputRPhi=1;
1080  mEventCounter++; // increase counter
1081 
1082  //mudst for vertex
1083  int flagvtx=0;
1084  StThreeVectorF vertex,vtxerr;
1085  const StMuDst* muDst = (const StMuDst*)GetInputDS("MuDst");
1086  if(muDst){
1087  StMuEvent *event = static_cast<StMuEvent*>(muDst->event());
1088  if(event){
1089  int nPrimV=muDst->numberOfPrimaryVertices();
1090  if(nPrimV>0) {
1091  StMuPrimaryVertex* V= muDst->primaryVertex(0); // select highest rank vertex
1092  vertex=V->position();
1093  vtxerr=V->posError();
1094  flagvtx=1;
1095  }
1096  }
1097  }
1098 
1099  //Getting track from Anselm's straight tracker
1100  StFgtStraightTrackMaker *fgtSTracker = static_cast<StFgtStraightTrackMaker * >( GetMaker("fgtStraightTracker"));
1101  if (!fgtSTracker){
1102  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No StFgtStraightTrackMaker" << endm;
1103  return; // if no event, we're done
1104  }
1105  vector<AVTrack>& tracks=fgtSTracker->getTracks();
1106  for(vector<AVTrack>::iterator t=tracks.begin();t!=tracks.end();t++){
1107  if(mDebug>0) cout<<Form("Trk chi2=%8.3f dca=%8.3f",t->chi2,t->dca)<<endl;
1108  vector<AVPoint>* points=t->points;
1109  int quad=-1, nhit=0, ntrk=-1;
1110  for(vector<AVPoint>::iterator p=points->begin(); p!=points->end();p++){
1111  int disc=p->dID;
1112  int quad2=p->quadID;
1113  if(quad==-1) quad=quad2;
1114  ntrk=mNtrk[quad];
1115  mHit[quad][ntrk].det[nhit]=disc*4+quad2;
1116  mHit[quad][ntrk].lr[nhit]=p->r;
1117  mHit[quad][ntrk].lp[nhit]=p->phi;
1118  mHit[quad][ntrk].z[nhit]=StFgtGeom::getDiscZ(disc);
1119  mHit[quad][ntrk].ex[nhit]=mErrFgt;
1120  mHit[quad][ntrk].ey[nhit]=mErrFgt;
1121  mHit[quad][ntrk].ez[nhit]=0;
1122  if(mDebug>0) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1123  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit],
1124  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1125  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1126  <<endl;
1127  nhit++;
1128  }
1129  if(flagvtx==1){
1130  mHit[quad][ntrk].det[nhit]=24+quad;
1131  mHit[quad][ntrk].x[nhit]=vertex.x();
1132  mHit[quad][ntrk].y[nhit]=vertex.y();
1133  mHit[quad][ntrk].z[nhit]=vertex.z();
1134  mHit[quad][ntrk].ex[nhit]=vtxerr.x();
1135  mHit[quad][ntrk].ey[nhit]=vtxerr.y();
1136  mHit[quad][ntrk].ez[nhit]=vtxerr.z();
1137  if(mDebug>0) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1138  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit],
1139  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1140  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1141  <<endl;
1142  nhit++;
1143  }
1144  mHit[quad][ntrk].nhit=nhit;
1145  mHit[quad][ntrk].dca =t->dca;
1146  mHit[quad][ntrk].chi2=t->chi2;
1147  mNtrk[quad]++;
1148  }
1149  if(mEventCounter%100==0) cout << Form("StFgtAlignmentMaker::readFromStraightTrackMaker EVT=%d NTRK=%d %d %d %d",
1150  mEventCounter,mNtrk[0],mNtrk[1],mNtrk[2],mNtrk[3])<<endl;
1151 }
1152 
1153 void StFgtAlignmentMaker::readFromStraightTrackAndStEvent(){
1154  float MaxDrTpc=maxdtpc, MaxDpTpc=maxptpc;
1155  float MaxDrEemc=maxdemc, MaxDpEemc=maxpemc;
1156  float MaxDz=5.0, MaxDca=3.0, MaxChi2=50;
1157 
1158  mFgtInputRPhi=1;
1159  mEventCounter++; // increase counter
1160 
1161  //Get StEvent
1162  StEvent* event;
1163  event = (StEvent *) GetInputDS("StEvent");
1164  if (!event){
1165  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No StEvent" << endm;
1166  return; // if no event, we're done
1167  }
1168  cout << "Event: Run "<< event->runId() << " Event No: " << event->id() << endl;
1169 
1170  //vertex
1171  int flagvtx=0;
1172  StThreeVectorF vertex,vtxerr;
1173  UInt_t NpVX = event->numberOfPrimaryVertices();
1174  cout << "Number of Primary vertex="<<NpVX<<endl;
1175  if(NpVX>0){
1176  const StPrimaryVertex *vx = event->primaryVertex(0); // select highest rank vertex
1177  cout << "Verrex :" << *vx << endl;
1178  vertex=vx->position();
1179  vtxerr=vx->positionError();
1180  flagvtx=1;
1181  }else{
1182  vertex.set(0.0,0.0,-999.0);
1183  vtxerr.set(999.0,999.0,999.0);
1184  //return;
1185  }
1186 
1187  //Getting track from Anselm's straight tracker
1188  StFgtStraightTrackMaker *fgtSTracker = static_cast<StFgtStraightTrackMaker * >( GetMaker("fgtStraightTracker"));
1189  if (!fgtSTracker){
1190  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No SyFgtStraightTrackMaker" << endm;
1191  return; // if no event, we're done
1192  }
1193  vector<AVTrack>& tracks=fgtSTracker->getTracks();
1194  for(vector<AVTrack>::iterator t=tracks.begin();t!=tracks.end();t++){
1195  float dz=-999.0;
1196  if(flagvtx){dz=t->trkZ - vertex.z();}
1197  if(mDebug>0) cout<<Form("Trk chi2=%8.3f dca=%8.3f zvtx=%8.3f dz=%8.3f",t->chi2,t->dca,t->trkZ,dz)<<endl;
1198  if(t->dca>MaxDca) continue;
1199  if(t->chi2>MaxChi2) continue;
1200  vector<AVPoint>* points=t->points;
1201  int quad=-1, nhit=0, ntrk=-1;
1202  //FGT hits on track
1203  for(vector<AVPoint>::iterator p=points->begin(); p!=points->end();p++){
1204  int disc=p->dID;
1205  int quad2=p->quadID;
1206  if(quad==-1) {quad=quad2; mQuad=quad;}
1207  ntrk=mNtrk[quad];
1208  mHit[quad][ntrk].det[nhit]=disc*4+quad2;
1209  mHit[quad][ntrk].lr[nhit]=p->r;
1210  mHit[quad][ntrk].lp[nhit]=p->phi;
1211  mHit[quad][ntrk].z[nhit]=StFgtGeom::getDiscZ(disc);
1212  mHit[quad][ntrk].ex[nhit]=mErrFgt;
1213  mHit[quad][ntrk].ey[nhit]=mErrFgt;
1214  mHit[quad][ntrk].ez[nhit]=0;
1215  mHit[quad][ntrk].tw[nhit]=p->phiCharge;
1216  mHit[quad][ntrk].p1[nhit]=p->rCharge;
1217  mHit[quad][ntrk].p2[nhit]=(p->phiCharge-p->rCharge)/(p->phiCharge+p->rCharge);
1218  mHit[quad][ntrk].po[nhit]=p->fgtHitPhi->getEvenOddChargeAsy();
1219  mHit[quad][ntrk].use[nhit]=1;
1220  if(mDebug>0) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1221  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit],
1222  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1223  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1224  <<endl;
1225  nhit++;
1226  }
1227  //vertex from StEvent
1228  if(flagvtx==1 && fabs(dz)<MaxDz){
1229  mHit[quad][ntrk].det[nhit]=24+quad;
1230  mHit[quad][ntrk].x[nhit]=vertex.x();
1231  mHit[quad][ntrk].y[nhit]=vertex.y();
1232  mHit[quad][ntrk].z[nhit]=vertex.z();
1233  mHit[quad][ntrk].r[nhit]=sqrt(vertex.x()*vertex.x()+vertex.y()*vertex.y());
1234  mHit[quad][ntrk].p[nhit]=atan2(vertex.y(),vertex.x());
1235  mHit[quad][ntrk].ex[nhit]=vtxerr.x();
1236  mHit[quad][ntrk].ey[nhit]=vtxerr.y();
1237  mHit[quad][ntrk].ez[nhit]=vtxerr.z();
1238  mHit[quad][ntrk].use[nhit]=1;
1239  if(mDebug>0) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1240  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit],
1241  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1242  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1243  <<endl;
1244  nhit++;
1245  }
1246  mHit[quad][ntrk].run=event->runId();
1247  mHit[quad][ntrk].evt=event->id();
1248  mHit[quad][ntrk].dca =t->dca;
1249  mHit[quad][ntrk].trkz =t->trkZ;
1250  mHit[quad][ntrk].dz =dz;
1251  mHit[quad][ntrk].chi2=t->chi2;
1252  mHit[quad][ntrk].nhit=nhit;
1253 
1254  //Refit with straight line with FGT and vertex
1255  int refit=1;
1256  double p[4];
1257  if(refit){
1258  fitTrackLine(orig_algpar,ntrk,p);
1259  if(mDebug>0){
1260  printf("FGT ONLY TRA CK : %8.4f %8.4f %8.4f %8.4f\n",t->ax,t->mx,t->ay,t->my);
1261  printf("REFIT with VTX : %8.4f %8.4f %8.4f %8.4f\n",p[0],p[1],p[2],p[3]);
1262  }
1263  }else{ p[0]=t->ax; p[1]=t->mx; p[2]=t->ay; p[3]=t->my;}
1264 
1265  //EEMC near FGT Track
1266  double ddr2=99999999;
1267  mHit[quad][ntrk].nhitEemc=0;
1268  StEEmcIUPointMaker *mEEpoints = (StEEmcIUPointMaker *)GetMaker("mEEpoints");
1269  if(mEEpoints){
1270  StEEmcIUPointVec_t mPoints = mEEpoints->points();
1271  for ( Int_t ipoint=0; ipoint<mEEpoints->numberOfPoints(); ipoint++){
1272  StEEmcIUPoint point=mEEpoints->point(ipoint);
1273  float x=point.position().x();
1274  float y=point.position().y();
1275  float z=point.position().z();
1276  float r=sqrt(x*x + y*y);
1277  float phi=atan2(y,x);
1278  float e[6]; int w[3];
1279  e[0]=point.energy(0);
1280  e[1]=point.energy(1)*1000;
1281  e[2]=point.energy(2)*1000;
1282  e[3]=point.energy(3)*1000;
1283  e[4]=point.cluster(0).energy()*1000;
1284  e[5]=point.cluster(1).energy()*1000;
1285  w[0]=point.numberOfRelatives();
1286  w[1]=point.cluster(0).numberOfStrips();
1287  w[2]=point.cluster(1).numberOfStrips();
1288  //printf("EEMC %d xyz=%6.2f %6.2f %6.2f e=%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n",
1289  // ipoint,x,y,z,e[neemc][0],e[neemc][1],e[neemc][2],e[neemc][3],e[neemc][4],e[neemc][5]);
1290  //double xx = t->ax + t->mx*z;
1291  //double yy = t->ay + t->my*z;
1292  double xx = p[0] + p[1]*z;
1293  double yy = p[2] + p[3]*z;
1294  double rr = sqrt(xx*xx+yy*yy);
1295  double pp = atan2(yy,xx);
1296  double dx = x-xx;
1297  double dy = y-yy;
1298  double dr = r-rr;
1299  double dp = phi-pp;
1300  while(dp>PI) dp-=2*PI;
1301  while(dp<-PI) dp+=2*PI;
1302  if(mDebug>3) printf("EEMC %3d xyrp=%6.1f %6.1f %6.1f %6.2f e=%6.1f %6.1f %6.1f %6.1f %6.1f %6.1f dxyrp=%6.1f %6.1f %6.1f %6.2f\n",
1303  ipoint,x,y,r,phi,e[0],e[1],e[2],e[3],e[4],e[5],dx,dy,dr,dp);
1304  if(dr < 0.5*rr-25.0 && //cut out high eta track
1305  fabs(dr)<MaxDrEemc && fabs(dp)<MaxDpEemc && e[1]>0.05 && e[2]>0.05){
1306  mHit[quad][ntrk].det[nhit]=36+quad;
1307  if(w[0]==0 && e[0]<1.5 && e[1]<4 && e[2]<4 && e[3]>0.05){ //MIP selection
1308  mHit[quad][ntrk].det[nhit]=40+quad;
1309  }
1310  if(e[0]>1 && e[1]>2 && e[2]>3 && e[2]>e[1] && (e[4]+e[5])/e[0]>10){ //Electron selection minus post
1311  mHit[quad][ntrk].det[nhit]=44+quad;
1312  }
1313  if(mHit[quad][ntrk].nhitEemc>0 && (dx*dx+dy*dy)>ddr2) continue;
1314  ddr2=dx*dx+dy*dy;
1315  mHit[quad][ntrk].nhitEemc=1;
1316  mHit[quad][ntrk].x[nhit]=x;
1317  mHit[quad][ntrk].y[nhit]=y;
1318  mHit[quad][ntrk].z[nhit]=z;
1319  mHit[quad][ntrk].r[nhit]=sqrt(x*x+y*y);
1320  mHit[quad][ntrk].p[nhit]=atan2(y,x);
1321  mHit[quad][ntrk].ex[nhit]=0.3;
1322  mHit[quad][ntrk].ey[nhit]=0.3;
1323  mHit[quad][ntrk].ez[nhit]=0.3;
1324  mHit[quad][ntrk].tw[nhit]=e[0];
1325  mHit[quad][ntrk].p1[nhit]=e[1];
1326  mHit[quad][ntrk].p2[nhit]=e[2];
1327  mHit[quad][ntrk].po[nhit]=e[3];
1328  mHit[quad][ntrk].su[nhit]=e[4];
1329  mHit[quad][ntrk].sv[nhit]=e[5];
1330  mHit[quad][ntrk].nrl[nhit]=w[0];
1331  mHit[quad][ntrk].nsu[nhit]=w[1];
1332  mHit[quad][ntrk].nsv[nhit]=w[2];
1333  if(mDebug>0) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1334  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit],
1335  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1336  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1337  <<endl;
1338  }
1339  }
1340  }
1341  if(mHit[quad][ntrk].nhitEemc>0) nhit++;
1342 
1343  //TPC near FGT track
1344  Long64_t ntot=0;
1345  ddr2=99999999;
1346  mHit[quad][ntrk].nhitPrompt=0;
1347  StTpcHitCollection* TpcHitCollection = event->tpcHitCollection();
1348  if(TpcHitCollection){
1349  //printf("got TpcHitCollection\n");
1350  UInt_t numberOfSectors = TpcHitCollection->numberOfSectors();
1351  for (UInt_t i = 0; i< numberOfSectors; i++) {
1352  StTpcSectorHitCollection* sectorCollection = TpcHitCollection->sector(i);
1353  if (!sectorCollection) continue;
1354  //printf("got tpc sector collection\n");
1355  Int_t numberOfPadrows = sectorCollection->numberOfPadrows();
1356  for (int j = 0; j< numberOfPadrows; j++) {
1357  StTpcPadrowHitCollection *rowCollection = sectorCollection->padrow(j);
1358  if (!rowCollection) continue;
1359  //printf("got tpc padrow collection\n");
1360  StSPtrVecTpcHit &hits = rowCollection->hits();
1361  Long64_t NoHits = hits.size();
1362  ntot += NoHits;
1363  //cout << "got nhit = "<<NoHits<<" / "<<ntot<<endl;
1364  for (Long64_t k = 0; k < NoHits; k++) {
1365  const StTpcHit *tpcHit = static_cast<const StTpcHit *> (hits[k]);
1366  const StThreeVectorF& xyz = tpcHit->position();
1367  if(xyz.z()<0) continue;
1368  if(xyz.z()<200.0) continue; //hack only prompt for now
1369  double xx = p[0] + p[1]*xyz.z();
1370  double yy = p[2] + p[3]*xyz.z();
1371  double rr = sqrt(xx*xx+yy*yy);
1372  double pp = atan2(yy,xx);
1373  double dx = xyz.x()-xx;
1374  double dy = xyz.y()-yy;
1375  double dr = xyz.perp()-rr;
1376  double dp = atan2(xyz.y(),xyz.x())-pp;
1377  while(dp> PI) dp-=2*PI;
1378  while(dp<-PI) dp+=2*PI;
1379  if(nhit>=MAXHIT){ printf("NOT ENOUGH SPACE FOR HIT!!!\n"); continue;}
1380  if(fabs(dr)<MaxDrTpc && fabs(dp)<MaxDpTpc && nhit<MAXHIT){
1381  mHit[quad][ntrk].det[nhit]=28+quad;
1382  if(xyz.z()>200) mHit[quad][ntrk].det[nhit]=32+quad; //prompt hit
1383  if(mHit[quad][ntrk].nhitPrompt>0 && (dx*dx+dy*dy)>ddr2) continue;
1384  ddr2=dx*dx+dy*dy;
1385  mHit[quad][ntrk].nhitPrompt=1;
1386  mHit[quad][ntrk].x[nhit]=xyz.x();
1387  mHit[quad][ntrk].y[nhit]=xyz.y();
1388  mHit[quad][ntrk].z[nhit]=xyz.z();
1389  mHit[quad][ntrk].r[nhit]=xyz.perp();
1390  mHit[quad][ntrk].p[nhit]=atan2(xyz.y(),xyz.x());
1391  mHit[quad][ntrk].ex[nhit]=0.2;
1392  mHit[quad][ntrk].ey[nhit]=0.2;
1393  mHit[quad][ntrk].ez[nhit]=0.2;
1394  if(mDebug>0) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1395  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit],
1396  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1397  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1398  <<endl;
1399  }
1400  }
1401  }
1402  }
1403  }
1404  if(mHit[quad][ntrk].nhitPrompt>0) nhit++;
1405  mHit[quad][ntrk].nhit=nhit;
1406  fillETBalance(ntrk);
1407  mNtrk[quad]++;
1408  }
1409  if(mEventCounter%100==0) cout << Form("StFgtAlignmentMaker::readFromStraightTrackMaker EVT=%d NTRK=%d %d %d %d",
1410  mEventCounter,mNtrk[0],mNtrk[1],mNtrk[2],mNtrk[3])<<endl;
1411 }
1412 
1413 void StFgtAlignmentMaker::readFromStraightTrackAndMudst(){
1414  float MaxDrEemc=maxdemc, MaxDpEemc=maxpemc;
1415  float MaxDz=10.0, MaxDca=5.0, MaxChi2=50;
1416 
1417  mFgtInputRPhi=1;
1418  mEventCounter++; // increase counter
1419 
1420  //Get MuDST
1421  const StMuDst* muDst = (const StMuDst*)GetInputDS("MuDst");
1422  if (!muDst){
1423  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No MuDST" << endm;
1424  return; // if no mudst, we're done
1425  }
1426  StMuEvent *event = static_cast<StMuEvent*>(muDst->event());
1427  if (!event){
1428  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No MuDST event" << endm;
1429  return; // if no mudst event, we're done
1430  }
1431  StEventInfo &info=event->eventInfo();
1432 
1433  if(mDebug>0) cout << "Event: Run "<< event->runId() << " Event No: " << event->eventId() << endl;
1434 
1435  //vertex
1436  int flagvtx=0;
1437  StThreeVectorF vertex,vtxerr;
1438  int nPrimV=muDst->numberOfPrimaryVertices();
1439  if(mDebug>2) cout << "Number of Primary vertex="<<nPrimV<<endl;
1440  if(nPrimV>0){
1441  StMuPrimaryVertex* vx= muDst->primaryVertex(0);
1442  vertex=vx->position();
1443  vtxerr=vx->posError();
1444  if(mDebug>2) cout << "Vertex = " << vertex.z() << endl;
1445  flagvtx=1;
1446  }else{
1447  return;
1448  }
1449 
1450  //Getting track from Anselm's straight tracker
1451  StFgtStraightTrackMaker *fgtSTracker = static_cast<StFgtStraightTrackMaker * >( GetMaker("fgtStraightTracker"));
1452  if (!fgtSTracker){
1453  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No SyFgtStraightTrackMaker" << endm;
1454  return; // if no event, we're done
1455  }
1456  vector<AVTrack>& tracks=fgtSTracker->getTracks();
1457  for(vector<AVTrack>::iterator t=tracks.begin();t!=tracks.end();t++){
1458  float dz=t->trkZ - vertex.z();
1459  if(mDebug>0) cout<<Form("Trk chi2=%8.3f dca=%8.3f zvtx=%8.3f dz=%8.3f",t->chi2,t->dca,t->trkZ,dz)<<endl;
1460  if(t->dca>MaxDca) continue;
1461  if(t->chi2>MaxChi2) continue;
1462  histdz->Fill(dz);
1463  if(fabs(dz)>MaxDz) continue;
1464  vector<AVPoint>* points=t->points;
1465  int quad=-1, nhit=0, ntrk=-1;
1466  //FGT hits on track
1467  for(vector<AVPoint>::iterator p=points->begin(); p!=points->end();p++){
1468  int disc=p->dID;
1469  int quad2=p->quadID;
1470  if(quad==-1) {quad=quad2; mQuad=quad;}
1471  ntrk=mNtrk[quad];
1472  mHit[quad][ntrk].run=event->runId();
1473  mHit[quad][ntrk].evt=event->eventId();
1474  mHit[quad][ntrk].det[nhit]=disc*4+quad2;
1475  mHit[quad][ntrk].lr[nhit]=p->r;
1476  mHit[quad][ntrk].lp[nhit]=p->phi;
1477  mHit[quad][ntrk].z[nhit]=StFgtGeom::getDiscZ(disc);
1478  mHit[quad][ntrk].ex[nhit]=mErrFgt;
1479  mHit[quad][ntrk].ey[nhit]=mErrFgt;
1480  mHit[quad][ntrk].ez[nhit]=0;
1481  mHit[quad][ntrk].tw[nhit]=p->phiCharge;
1482  mHit[quad][ntrk].p1[nhit]=p->rCharge;
1483  mHit[quad][ntrk].p2[nhit]=(p->phiCharge-p->rCharge)/(p->phiCharge+p->rCharge);
1484  mHit[quad][ntrk].po[nhit]=p->fgtHitPhi->getEvenOddChargeAsy();
1485  mHit[quad][ntrk].use[nhit]=1;
1486  if(mDebug>0) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Disc=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1487  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit]/4,
1488  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1489  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1490  <<endl;
1491  nhit++;
1492  }
1493  //vertex from mudst
1494  if(flagvtx==1){
1495  mHit[quad][ntrk].det[nhit]=24+quad;
1496  mHit[quad][ntrk].x[nhit]=vertex.x();
1497  mHit[quad][ntrk].y[nhit]=vertex.y();
1498  mHit[quad][ntrk].z[nhit]=vertex.z();
1499  mHit[quad][ntrk].r[nhit]=sqrt(vertex.x()*vertex.x()+vertex.y()*vertex.y());
1500  mHit[quad][ntrk].p[nhit]=atan2(vertex.y(),vertex.x());
1501  mHit[quad][ntrk].ex[nhit]=vtxerr.x();
1502  mHit[quad][ntrk].ey[nhit]=vtxerr.y();
1503  mHit[quad][ntrk].ez[nhit]=vtxerr.z();
1504  mHit[quad][ntrk].use[nhit]=0;
1505  if(mDebug>-1) cout<<Form("Trk=%3d Hit=%3d Quad=%1d Disc=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1506  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit]/4,
1507  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1508  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1509  <<endl;
1510  nhit++;
1511  }
1512  mHit[quad][ntrk].dca =t->dca;
1513  mHit[quad][ntrk].trkz =t->trkZ;
1514  mHit[quad][ntrk].dz =dz;
1515  mHit[quad][ntrk].chi2=t->chi2;
1516  mHit[quad][ntrk].nhit=nhit;
1517 
1518  //mudst Refit with straight line with FGT and vertex
1519  int refit=0;
1520  double p[4];
1521  if(refit){
1522  int flg;
1523  double arg[10],r, e;
1524  TMinuit *m = new TMinuit(4);
1525  m->SetFCN(fitLine);
1526  arg[0] =-1; m->mnexcm("SET PRI", arg, 1,flg);
1527  arg[0] = 1; m->mnexcm("SET ERR", arg, 1,flg);
1528  arg[0] = 500; arg[1] = 1.;
1529  //make "aligned" hits in mHitAlg
1530  getAlign(ntrk,orig_algpar);
1531  //first guess of line fit parameter and set up
1532  double x1 = (mHitAlg.x[mHitAlg.nhit-1]-mHitAlg.x[0])/(mHitAlg.z[mHitAlg.nhit-1]-mHitAlg.z[0]);
1533  double y1 = (mHitAlg.y[mHitAlg.nhit-1]-mHitAlg.y[0])/(mHitAlg.z[mHitAlg.nhit-1]-mHitAlg.z[0]);
1534  double x0 = mHitAlg.x[0] - x1*mHitAlg.z[0];
1535  double y0 = mHitAlg.y[0] - y1*mHitAlg.z[0];
1536  m->mnparm(0,"x0",x0,0.1,0,0,flg);
1537  m->mnparm(1,"x1",x1,0.1,0,0,flg);
1538  m->mnparm(2,"y0",y0,0.1,0,0,flg);
1539  m->mnparm(3,"y1",y1,0.1,0,0,flg);
1540  //do straight line fit
1541  m->mnexcm("MIGRAD", arg ,2, flg);
1542  for(int j=0; j<4; j++){ m->GetParameter(j,r,e); p[j]=r; }
1543  if(mDebug>0){
1544  printf("FGT ONLY TRACK : %8.4f %8.4f %8.4f %8.4f\n",t->ax,t->mx,t->ay,t->my);
1545  printf("REFIT with VTX : %8.4f %8.4f %8.4f %8.4f\n",p[0],p[1],p[2],p[3]);
1546  }
1547  }else{
1548  p[0]=t->ax; p[1]=t->mx; p[2]=t->ay; p[3]=t->mx;
1549  }
1550 
1551  //EEMC near FGT Track
1552  mHit[quad][ntrk].nhitEemc=0;
1553  StEEmcIUPointMaker *mEEpoints = (StEEmcIUPointMaker *)GetMaker("mEEpoints");
1554  if(!mEEpoints) {printf("No EEMC Point\n");}
1555  else{
1556  StEEmcIUPointVec_t mPoints = mEEpoints->points();
1557  for ( Int_t ipoint=0; ipoint<mEEpoints->numberOfPoints(); ipoint++){
1558  StEEmcIUPoint point=mEEpoints->point(ipoint);
1559  float x=point.position().x();
1560  float y=point.position().y();
1561  float z=point.position().z();
1562  float r=sqrt(x*x + y*y);
1563  float phi=atan2(y,x);
1564  float e[6]; int w[3];
1565  e[0]=point.energy(0);
1566  e[1]=point.energy(1)*1000;
1567  e[2]=point.energy(2)*1000;
1568  e[3]=point.energy(3)*1000;
1569  e[4]=point.cluster(0).energy()*1000;
1570  e[5]=point.cluster(1).energy()*1000;
1571  w[0]=point.numberOfRelatives();
1572  w[1]=point.cluster(0).numberOfStrips();
1573  w[2]=point.cluster(1).numberOfStrips();
1574  //double xx = t->ax + t->mx*z;
1575  //double yy = t->ay + t->my*z;
1576  double xx = p[0] + p[1]*z;
1577  double yy = p[2] + p[3]*z;
1578  double rr = sqrt(xx*xx + yy*yy);
1579  double pp = atan2(yy,xx);
1580  double dx = x-xx;
1581  double dy = y-yy;
1582  double dr = r-rr;
1583  double dp = phi-pp;
1584  while(dp>PI) dp-=2*PI;
1585  while(dp<-PI) dp+=2*PI;
1586  if(mDebug>3) printf("EEMC %3d xyrp=%6.1f %6.1f %6.1f %6.2f e=%6.1f %6.1f %6.1f %6.1f %6.1f %6.1f dxyrp=%6.1f %6.1f %6.1f %6.2f\n",
1587  ipoint,x,y,r,phi,e[0],e[1],e[2],e[3],e[4],e[5],dx,dy,dr,dp);
1588  if(dr < 0.5*rr-25.0 && //cut out low _track
1589  fabs(dr)<MaxDrEemc && fabs(dp)<MaxDpEemc && e[1]>0.05 && e[2]>0.05){
1590  mHit[quad][ntrk].det[nhit]=36+quad;
1591  if(w[0]==0 && e[0]<1.5 && e[1]<4 && e[2]<4 && e[3]>0.05){ //MIP selection
1592  mHit[quad][ntrk].det[nhit]=40+quad;
1593  }
1594  if(e[0]>1 && e[1]>2 && e[2]>3 && e[2]>e[1] && (e[4]+e[5])/e[0]>10){ //Electron selection minus post
1595  mHit[quad][ntrk].det[nhit]=44+quad;
1596  }
1597  mHit[quad][ntrk].nhitEemc++;
1598  mHit[quad][ntrk].x[nhit]=x;
1599  mHit[quad][ntrk].y[nhit]=y;
1600  mHit[quad][ntrk].z[nhit]=z;
1601  mHit[quad][ntrk].r[nhit]=sqrt(x*x+y*y);
1602  mHit[quad][ntrk].p[nhit]=atan2(y,x);
1603  mHit[quad][ntrk].ex[nhit]=0.5;
1604  mHit[quad][ntrk].ey[nhit]=0.5;
1605  mHit[quad][ntrk].ez[nhit]=0.5;
1606  mHit[quad][ntrk].tw[nhit]=e[0];
1607  mHit[quad][ntrk].p1[nhit]=e[1];
1608  mHit[quad][ntrk].p2[nhit]=e[2];
1609  mHit[quad][ntrk].po[nhit]=e[3];
1610  mHit[quad][ntrk].su[nhit]=e[4];
1611  mHit[quad][ntrk].sv[nhit]=e[5];
1612  mHit[quad][ntrk].nrl[nhit]=w[0];
1613  mHit[quad][ntrk].nsu[nhit]=w[1];
1614  mHit[quad][ntrk].nsv[nhit]=w[2];
1615  if(mDebug>0) cout<<Form("Trk=%3d EHit=%3d Quad=%1d Disc=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
1616  ntrk,nhit,quad,mHit[quad][ntrk].det[nhit]/4,
1617  mHit[quad][ntrk].x[nhit],mHit[quad][ntrk].y[nhit],mHit[quad][ntrk].z[nhit],
1618  mHit[quad][ntrk].ex[nhit],mHit[quad][ntrk].ey[nhit],mHit[quad][ntrk].ez[nhit])
1619  <<endl;
1620  if(mDebug>0) printf("EEMC %3d q=%1d d=%1d trk=%6.2f %6.2f %6.2f %6.2f fgt=%6.1f %6.1f %6.1f %6.2f xyrp=%6.1f %6.1f %6.1f %6.2f dxyrp=%6.1f %6.1f %6.1f %6.2f\n",
1621  ntrk,quad,mHit[quad][ntrk].det[nhit]/4,p[0],p[1],p[2],p[3],xx,yy,rr,pp,x,y,r,phi,dx,dy,dr,dp);
1622  nhit++;
1623  }
1624  }
1625  }
1626  mHit[quad][ntrk].nhit=nhit;
1627  //save only eemc correlated!!!
1628  //if(mHit[quad][ntrk].nhitEemc==0) continue;
1629  mNtrk[quad]++;
1630  }
1631  if(mEventCounter%100==0) cout << Form("StFgtAlignmentMaker::readFromStraightTrackMaker EVT=%d NTRK=%d %d %d %d",
1632  mEventCounter,mNtrk[0],mNtrk[1],mNtrk[2],mNtrk[3])<<endl;
1633 }
1634 
1635 void StFgtAlignmentMaker::DispFromStraightTrackMaker(){
1636  static int nplot=0;
1637  float zmin=-50, zmin2=100;
1638  float zmax=300;
1639  float maxDist=10.0;
1640  float maxDCA=5.0;
1641  //Getting track from Anselm's straight tracker
1642  StFgtStraightTrackMaker *fgtSTracker = static_cast<StFgtStraightTrackMaker * >( GetMaker("fgtStraightTracker"));
1643  if (!fgtSTracker){
1644  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No SyFgtStraightTrackMaker" << endm;
1645  return; // if no event, we're done
1646  }
1647  vector<AVTrack>& tracks=fgtSTracker->getTracks();
1648  for(vector<AVTrack>::iterator t=tracks.begin();t!=tracks.end();t++){
1649  //fgt track
1650  int quad=-1;
1651  if(t->dca>maxDCA) continue;
1652  if(t->trkZ<-200 || t->trkZ>200) continue;
1653  cout<<Form("Trk chi2=%9.6f dca=%7.2f z=%7.2f mx=%7.2f ax=%7.2f my=%7.2f ay=%7.2f",
1654  t->chi2,t->dca,t->trkZ,t->mx,t->ax,t->my,t->ay)<<endl;
1655  TGraph *hitxz=new TGraph(1); hitxz->SetMarkerStyle(21); hitxz->SetMarkerSize(0.8); hitxz->SetMarkerColor(kRed);
1656  TGraph *hityz=new TGraph(1); hityz->SetMarkerStyle(21); hityz->SetMarkerSize(0.8); hityz->SetMarkerColor(kRed);
1657  TGraph *hitrz=new TGraph(1); hitrz->SetMarkerStyle(21); hitrz->SetMarkerSize(0.8); hitrz->SetMarkerColor(kRed);
1658  TGraph *hitpz=new TGraph(1); hitpz->SetMarkerStyle(21); hitpz->SetMarkerSize(0.8); hitpz->SetMarkerColor(kRed);
1659  TF1 *fxz = new TF1("fxz","[0]+[1]*x",zmin,zmax); fxz->SetParameter(0,t->ax); fxz->SetParameter(1,t->mx);
1660  TF1 *fyz = new TF1("fyz","[0]+[1]*x",zmin,zmax); fyz->SetParameter(0,t->ay); fyz->SetParameter(1,t->my);
1661  TF1 *frz = new TF1("frz","sqrt(([0]+[1]*x)*([0]+[1]*x)+([2]+[3]*x)*([2]+[3]*x))",zmin,zmax);
1662  frz->SetParameter(0,t->ax); frz->SetParameter(1,t->mx); frz->SetParameter(2,t->ay); frz->SetParameter(3,t->my);
1663  TF1 *fpz = new TF1("fpz","atan2([2]+[3]*x,[0]+[1]*x)",zmin,zmax);
1664  fpz->SetParameter(0,t->ax); fpz->SetParameter(1,t->mx); fpz->SetParameter(2,t->ay); fpz->SetParameter(3,t->my);
1665  vector<AVPoint>* points=t->points;
1666  int i=0;
1667  //fgt hits
1668  for(vector<AVPoint>::iterator p=points->begin(); p!=points->end();p++){
1669  quad=p->quadID;
1670  hitxz->SetPoint(i,p->z,p->x);
1671  hityz->SetPoint(i,p->z,p->y);
1672  hitrz->SetPoint(i,p->z,p->r);
1673  hitpz->SetPoint(i,p->z,p->phi);
1674  cout << Form(" %2d quad=%d xyz=%8.3f %8.3f %8.3f",i,quad,p->x,p->y,p->z)<<endl;
1675  i++;
1676  }
1677  //TPC
1678  Long64_t ntot=0;
1679  int ndisp=0, nprompt=0, neemc=0;
1680  TGraph *tpcxz=new TGraph(1); tpcxz->SetMarkerStyle(21); tpcxz->SetMarkerSize(0.8); tpcxz->SetMarkerColor(kBlue);
1681  TGraph *tpcyz=new TGraph(1); tpcyz->SetMarkerStyle(21); tpcyz->SetMarkerSize(0.8); tpcyz->SetMarkerColor(kBlue);
1682  TGraph *tpcrz=new TGraph(1); tpcrz->SetMarkerStyle(21); tpcrz->SetMarkerSize(0.8); tpcrz->SetMarkerColor(kBlue);
1683  TGraph *tpcpz=new TGraph(1); tpcpz->SetMarkerStyle(21); tpcpz->SetMarkerSize(0.8); tpcpz->SetMarkerColor(kBlue);
1684  TGraph *tpcdxz=new TGraph(1); tpcdxz->SetMarkerStyle(21); tpcdxz->SetMarkerSize(0.8); tpcdxz->SetMarkerColor(kBlue);
1685  TGraph *tpcdyz=new TGraph(1); tpcdyz->SetMarkerStyle(21); tpcdyz->SetMarkerSize(0.8); tpcdyz->SetMarkerColor(kBlue);
1686  TGraph *tpcdrz=new TGraph(1); tpcdrz->SetMarkerStyle(21); tpcdrz->SetMarkerSize(0.8); tpcdrz->SetMarkerColor(kBlue);
1687  TGraph *tpcdpz=new TGraph(1); tpcdpz->SetMarkerStyle(21); tpcdpz->SetMarkerSize(0.8); tpcdpz->SetMarkerColor(kBlue);
1688  StEvent* event;
1689  event = (StEvent *) GetInputDS("StEvent");
1690  if(event){
1691  //printf("got StEvent\n");
1692  StTpcHitCollection* TpcHitCollection = event->tpcHitCollection();
1693  if(TpcHitCollection){
1694  //printf("got TpcHitCollection\n");
1695  UInt_t numberOfSectors = TpcHitCollection->numberOfSectors();
1696  for (UInt_t i = 0; i< numberOfSectors; i++) {
1697  StTpcSectorHitCollection* sectorCollection = TpcHitCollection->sector(i);
1698  if (!sectorCollection) continue;
1699  //printf("got tpc sector collection\n");
1700  Int_t numberOfPadrows = sectorCollection->numberOfPadrows();
1701  for (int j = 0; j< numberOfPadrows; j++) {
1702  StTpcPadrowHitCollection *rowCollection = sectorCollection->padrow(j);
1703  if (!rowCollection) continue;
1704  //printf("got tpc padrow collection\n");
1705  StSPtrVecTpcHit &hits = rowCollection->hits();
1706  Long64_t NoHits = hits.size();
1707  ntot += NoHits;
1708  //cout << "got nhit = "<<NoHits<<" / "<<ntot<<endl;
1709  for (Long64_t k = 0; k < NoHits; k++) {
1710  const StTpcHit *tpcHit = static_cast<const StTpcHit *> (hits[k]);
1711  const StThreeVectorF& xyz = tpcHit->position();
1712  if(xyz.z()<0) continue;
1713  double xx = t->ax + t->mx*xyz.z();
1714  double yy = t->ay + t->my*xyz.z();
1715  double rr = sqrt(xx*xx+yy*yy);
1716  double pp = atan2(yy,xx);
1717  double dx = xyz.x()-xx;
1718  double dy = xyz.y()-yy;
1719  double dr = xyz.perp()-rr;
1720  double dp = atan2(xyz.y(),xyz.x())-pp;
1721  static const double PI=3.141592654;
1722  if(dp>PI) dp-=2*PI;
1723  if(dp<-PI) dp+=2*PI;
1724  double dist = sqrt(dx*dx+dy*dy);
1725  if(dist<maxDist){
1726  tpcxz->SetPoint(ndisp,xyz.z(),xyz.x());
1727  tpcyz->SetPoint(ndisp,xyz.z(),xyz.y());
1728  tpcrz->SetPoint(ndisp,xyz.z(),xyz.perp());
1729  tpcpz->SetPoint(ndisp,xyz.z(),atan2(xyz.y(),xyz.x()));
1730  tpcdxz->SetPoint(ndisp,xyz.z(),dx);
1731  tpcdyz->SetPoint(ndisp,xyz.z(),dy);
1732  tpcdrz->SetPoint(ndisp,xyz.z(),dr);
1733  tpcdpz->SetPoint(ndisp,xyz.z(),dp);
1734  ndisp++;
1735  if(xyz.z()>200) nprompt++;
1736  fillHist(6,quad,dx,dy,dr,dp,xx,yy,rr,pp);
1737  }
1738  }
1739  }
1740  }
1741  }
1742  }
1743  cout << "Got TPC hits near fgt track="<<ndisp<<" out of total="<<ntot<<endl;
1744  //EEMC
1745  ndisp=0;
1746  float e[10][6];
1747  TGraph *emcxz=new TGraph(1); emcxz->SetMarkerStyle(21); emcxz->SetMarkerSize(0.8); emcxz->SetMarkerColor(kRed);
1748  TGraph *emcyz=new TGraph(1); emcyz->SetMarkerStyle(21); emcyz->SetMarkerSize(0.8); emcyz->SetMarkerColor(kRed);
1749  TGraph *emcrz=new TGraph(1); emcrz->SetMarkerStyle(21); emcrz->SetMarkerSize(0.8); emcrz->SetMarkerColor(kRed);
1750  TGraph *emcpz=new TGraph(1); emcpz->SetMarkerStyle(21); emcpz->SetMarkerSize(0.8); emcpz->SetMarkerColor(kRed);
1751  TGraph *emcdxz=new TGraph(1); emcdxz->SetMarkerStyle(21); emcdxz->SetMarkerSize(0.8); emcdxz->SetMarkerColor(kRed);
1752  TGraph *emcdyz=new TGraph(1); emcdyz->SetMarkerStyle(21); emcdyz->SetMarkerSize(0.8); emcdyz->SetMarkerColor(kRed);
1753  TGraph *emcdrz=new TGraph(1); emcdrz->SetMarkerStyle(21); emcdrz->SetMarkerSize(0.8); emcdrz->SetMarkerColor(kRed);
1754  TGraph *emcdpz=new TGraph(1); emcdpz->SetMarkerStyle(21); emcdpz->SetMarkerSize(0.8); emcdpz->SetMarkerColor(kRed);
1755  StEEmcIUPointMaker *mEEpoints = (StEEmcIUPointMaker *)GetMaker("mEEpoints");
1756  if(mEEpoints){
1757  StEEmcIUPointVec_t mPoints = mEEpoints->points();
1758  for ( Int_t ipoint=0; ipoint<mEEpoints->numberOfPoints(); ipoint++){
1759  StEEmcIUPoint point=mEEpoints->point(ipoint);
1760  float x=point.position().x();
1761  float y=point.position().y();
1762  float z=point.position().z();
1763  e[neemc][0]=point.energy(0);
1764  e[neemc][1]=point.energy(1)*1000;
1765  e[neemc][2]=point.energy(2)*1000;
1766  e[neemc][3]=point.energy(3)*1000;
1767  e[neemc][4]=point.cluster(0).energy()*1000;
1768  e[neemc][5]=point.cluster(1).energy()*1000;
1769  printf("EEMC %d xyz=%6.2f %6.2f %6.2f e=%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n",
1770  ipoint,x,y,z,e[neemc][0],e[neemc][1],e[neemc][2],e[neemc][3],e[neemc][4],e[neemc][5]);
1771  double xx = t->ax + t->mx*z;
1772  double yy = t->ay + t->my*z;
1773  double rr = sqrt(xx*xx+yy*yy);
1774  double pp = atan2(yy,xx);
1775  double dx = x-xx;
1776  double dy = y-yy;
1777  double dr = sqrt(x*x+y*y)-rr;
1778  double dp = atan2(y,x)-pp;
1779  static const double PI=3.141592654;
1780  if(dp>PI) dp-=2*PI;
1781  if(dp<-PI) dp+=2*PI;
1782  double dist = sqrt(dx*dx+dy*dy);
1783  if(dist<maxDist && e[neemc][1]>0.3 && e[neemc][2]>0.3){
1784  emcxz->SetPoint(ndisp,z,x);
1785  emcyz->SetPoint(ndisp,z,y);
1786  emcrz->SetPoint(ndisp,z,sqrt(x*x+y*y));
1787  emcpz->SetPoint(ndisp,z,atan2(y,x));
1788  emcdxz->SetPoint(ndisp,z,dx);
1789  emcdyz->SetPoint(ndisp,z,dy);
1790  emcdrz->SetPoint(ndisp,z,dr);
1791  emcdpz->SetPoint(ndisp,z,dp);
1792  ndisp++;
1793  neemc++;
1794  fillHist(7,quad,dx,dy,dr,dp,xx,yy,rr,pp);
1795  }
1796  }
1797  }
1798  //Draw
1799  if(ndisp==0) continue;
1800  //if(nprompt==0) continue;
1801  if(neemc==0) continue;
1802  //if(nprompt==0 && neemc==0) continue;
1803  //getting idea where we should plot
1804  double xmin = (t->ax + t->mx*zmin2);
1805  double ymin = (t->ay + t->my*zmin2);
1806  double xmax = (t->ax + t->mx*zmax);
1807  double ymax = (t->ay + t->my*zmax);
1808  double rmin = sqrt(xmin*xmin+ymin*ymin);
1809  double rmax = sqrt(xmax*xmax+ymax*ymax);
1810  double pmax = atan2(ymax,xmax) + 0.3;
1811  double pmin = atan2(ymax,xmax) - 0.3;
1812  if(rmax<70) continue;
1813  //Draw
1814  static TCanvas *c1=0, *c2=0;
1815  if(!c1) c1=new TCanvas("EventDisp","Event Disp",50,50,800,800);
1816  if(!c2) c2=new TCanvas("EventDisp2","Event Disp2",100,100,800,800);
1817  gStyle->SetOptStat(0);
1818  //far view
1819  TH2F *xz = new TH2F("xz","xz",1,zmin,zmax,1,-200,200);
1820  TH2F *yz = new TH2F("yz","yz",1,zmin,zmax,1,-200,200);
1821  TH2F *rz = new TH2F("rz","rz",1,zmin,zmax,1, -10,200);
1822  TH2F *pz = new TH2F("pz","pz",1,zmin,zmax,1,-3.2,3.2);
1823  c1->Clear();
1824  c1->Divide(2,2);
1825  c1->cd(1); xz->Draw(); hitxz->Draw("P"); fxz->Draw("LSAME"); tpcxz->Draw("P"); emcxz->Draw("P");
1826  c1->cd(2); yz->Draw(); hityz->Draw("P"); fyz->Draw("LSAME"); tpcyz->Draw("P"); emcyz->Draw("P");
1827  c1->cd(3); rz->Draw(); hitrz->Draw("P"); frz->Draw("LSAME"); tpcrz->Draw("P"); emcrz->Draw("P");
1828  c1->cd(4); pz->Draw(); hitpz->Draw("P"); fpz->Draw("LSAME"); tpcpz->Draw("P"); emcpz->Draw("P");
1829  c1->Update();
1830  //close view for TPC
1831  //if(xmax<xmin) {double t=xmin; xmin=xmax; xmax=t;}
1832  //if(ymax<ymin) {double t=ymin; ymin=ymax; ymax=t;}
1833  //float ddx=fabs(xmax-xmin), avx=(xmax+xmin)/2.0;
1834  //float ddy=fabs(ymax-ymin), avy=(ymax+ymin)/2.0;
1835  //if(ddx<20) {xmin=avx-10.0; xmax=avx+10.0;}
1836  //if(ddy<20) {ymin=avy-10.0; ymax=avy+10.0;}
1837  //TH2F *xz2 = new TH2F("xz","xz",1,zmin2,zmax,1,xmin,xmax);
1838  //TH2F *yz2 = new TH2F("yz","yz",1,zmin2,zmax,1,ymin,ymax);
1839  //TH2F *rz2 = new TH2F("rz","rz",1,zmin2,zmax,1,rmin,rmax);
1840  //TH2F *pz2 = new TH2F("pz","pz",1,zmin2,zmax,1,pmin,pmax);
1841  TH2F *xz2 = new TH2F("dxz","dxz",1,zmin2,zmax,1,-12,12);
1842  TH2F *yz2 = new TH2F("dyz","dyz",1,zmin2,zmax,1,-12,12);
1843  TH2F *rz2 = new TH2F("drz","drz",1,zmin2,zmax,1,-12,12);
1844  TH2F *pz2 = new TH2F("dpz","dpz",1,zmin2,zmax,1,-0.3,0.3);
1845  c2->Clear();
1846  c2->Divide(2,2);
1847  c2->cd(1); xz2->Draw(); /* hitxz->Draw("P"); fxz->Draw("LSAME"); */ tpcdxz->Draw("P"); emcdxz->Draw("P");
1848  c2->cd(2); yz2->Draw(); /* hityz->Draw("P"); fyz->Draw("LSAME"); */ tpcdyz->Draw("P"); emcdyz->Draw("P");
1849  c2->cd(3); rz2->Draw(); /* hitrz->Draw("P"); frz->Draw("LSAME"); */ tpcdrz->Draw("P"); emcdrz->Draw("P");
1850  c2->cd(4); pz2->Draw(); /* hitpz->Draw("P"); fpz->Draw("LSAME"); */ tpcdpz->Draw("P"); emcdpz->Draw("P");
1851  for(int i=0; i<neemc; i++){
1852  char t[100];
1853  sprintf(t,"E=%4.1f P1=%4.2f P2=%4.2f T=%4.2f U=%4.2f V=%4.2f",e[i][0],e[i][1],e[i][2],e[i][3],e[i][4],e[i][5]);
1854  printf("%s\n",t);
1855  TText *tt = new TText(0.1,0.9-0.07*i,t); tt->SetNDC(); tt->SetTextSize(0.04); tt->Draw();
1856  }
1857  c2->Update();
1858  if(nplot<50){
1859  char tmp[100];
1860  sprintf(tmp,"disp/evdisp1_%d.png",nplot); c1->SaveAs(tmp);
1861  sprintf(tmp,"disp/evdisp2_%d.png",nplot); c2->SaveAs(tmp);
1862  nplot++;
1863  }
1864  //printf("Hit something -->");
1865  //cin.get();
1866  //printf(" Contineing...\n");
1867  }
1868 }
1869 
1870 void StFgtAlignmentMaker::calcETBalanceFromStEvent() {
1871  StEvent* event;
1872  event = (StEvent *) GetInputDS("StEvent");
1873  if (!event){
1874  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No StEvent" << endm;
1875  return; // if no event, we're done
1876  }
1877  ele.set(0,0,0); tracks.set(0,0,0); bemcs.set(0,0,0); eemcs.set(0,0,0); tot.set(0,0,0); rec.set(0,0,0); iso.set(0,0,0);
1878  int nele=0,ntrack=0, nbemc=0, neemc=0, ntot=0, nrec=0, niso=0;
1879  UInt_t NpVX = event->numberOfPrimaryVertices();
1880  if (!NpVX) return;
1881  NpVX=1; //only 1st vertex for now
1882  for (UInt_t i = 0; i < NpVX; i++) {
1883  const StPrimaryVertex *vtx = event->primaryVertex(i);
1884  if(mDebug>2) cout << Form("Vertex: %3i ",i) << *vtx << endl;
1885  //EEMC find W candidate
1886  StEEmcIUPointMaker *mEEpoints = (StEEmcIUPointMaker *)GetMaker("mEEpoints");
1887  if(mEEpoints){
1888  StEEmcIUPointVec_t mPoints = mEEpoints->points();
1889  for ( Int_t ipoint=0; ipoint<mEEpoints->numberOfPoints(); ipoint++){
1890  StEEmcIUPoint point=mEEpoints->point(ipoint);
1891  float e=point.energy(0);
1892  StThreeVectorF pp(point.position().x(),point.position().y(),point.position().z());
1893  StThreeVectorF d=pp - vtx->position();
1894  d=d/d.mag();
1895  StThreeVectorF ee=e*d;
1896  eemcs+=ee;
1897  if(e>10.0 &&
1898  // point.numberOfRelatives()==0 &&
1899  point.cluster(0).energy()>0.02 &&
1900  point.cluster(1).energy()>0.02){
1901  if(ee.perp() > ele.perp()) ele=ee;
1902  nele++;
1903  }
1904  }
1905  }
1906  //EEMC isolation cone around W candidate
1907  if(mEEpoints && nele>0){
1908  StEEmcIUPointVec_t mPoints = mEEpoints->points();
1909  for ( Int_t ipoint=0; ipoint<mEEpoints->numberOfPoints(); ipoint++){
1910  StEEmcIUPoint point=mEEpoints->point(ipoint);
1911  float e=point.energy(0);
1912  StThreeVectorF pp(point.position().x(),point.position().y(),point.position().z());
1913  StThreeVectorF d=pp - vtx->position();
1914  d=d/d.mag();
1915  neemc++;
1916  double deta = ele.pseudoRapidity() - d.pseudoRapidity();
1917  double dphi = ele.phi() - d.phi(); if(dphi>PI) {dphi-=2*PI;} if(dphi<-PI) {dphi+=2*PI;}
1918  if(fabs(deta)<ISOCUT && fabs(dphi)<ISOCUT) { StThreeVectorF ee=e*d; iso+=ee; niso++;}
1919  }
1920  }
1921  //Tracks
1922  UInt_t nDaughters = vtx->numberOfDaughters();
1923  for (UInt_t j = 0; j < nDaughters; j++) {
1924  StPrimaryTrack* pTrack = (StPrimaryTrack*) vtx->daughter(j);
1925  if (! pTrack) continue;
1926  tracks+=pTrack->geometry()->momentum();
1927  ntrack++;
1928  double deta = ele.pseudoRapidity() - pTrack->geometry()->momentum().pseudoRapidity();
1929  double dphi = ele.phi() - pTrack->geometry()->momentum().phi(); if(dphi>PI) {dphi-=2*PI;} if(dphi<-PI) {dphi+=2*PI;}
1930  if(fabs(deta)<ISOCUT && fabs(dphi)<ISOCUT) { iso+=pTrack->geometry()->momentum(); }
1931  }
1932  //BEMC
1933  StEmcDetector *btow=event->emcCollection()->detector(kBarrelEmcTowerId);
1934  if(btow){
1935  for(int i=0; i<btow->numberOfModules(); i++){
1936  StEmcModule* mod=btow->module(i+1);
1937  StSPtrVecEmcRawHit& hits = mod->hits();
1938  for(UInt_t k=0;k<hits.size();k++) {
1939  float e=hits[k]->energy();
1940  if(e<0.2) continue;
1941  int mod=hits[k]->module();
1942  int eta=hits[k]->eta();
1943  int sub=hits[k]->sub();
1944  float x,y,z;
1945  StEmcGeom::instance("bemc")->getXYZ(mod,eta,sub,x,y,z);
1946  //printf("mod=%3d eta=%3d sub=%3d ene=%8.2f xyz=%8.2f %8.2f %8.2f\n",mod,eta,sub,e,x,y,z);
1947  StThreeVectorF d = StThreeVectorF(x,y,z) - vtx->position();
1948  d=d/d.mag();
1949  bemcs+=e*d;
1950  nbemc++;
1951  double deta = ele.pseudoRapidity() - d.pseudoRapidity();
1952  double dphi = ele.phi() - d.phi(); if(dphi>PI) {dphi-=2*PI;} if(dphi<-PI) {dphi+=2*PI;}
1953  if(fabs(deta)<ISOCUT && fabs(dphi)<ISOCUT) { StThreeVectorF ee=e*d; iso+=ee; }
1954  }
1955  }
1956  }
1957  }
1958  tot=tracks+bemcs+eemcs; ntot=ntrack+nbemc+neemc;
1959  rec=tot-ele; if(nele>0) {nrec=ntot-1;} else {nrec=ntot;}
1960  float dphi=ele.phi() - tot.phi();
1961  ptbal = tot.perp()*cos(dphi);
1962  riso = ele.mag()/iso.mag();
1963  printf("NELE=%4d NTRK=%4d NBEMC=%4d NEEMC=%4d NTOT=%4d NREC=%4d NISO=%4d\n",nele,ntrack,nbemc,neemc,ntot,nrec,niso);
1964  printf("EELE=%6.1f ETRK=%6.1f EBEMC=%6.1f EEEMC=%6.1f ETOT=%6.1f EREC=%6.1f EISO=%6.1f\n",
1965  ele.mag(),tracks.mag(),bemcs.mag(),eemcs.mag(),tot.mag(),rec.mag(),iso.mag());
1966  printf("PELE=%6.1f PTRK=%6.1f PBEMC=%6.1f PEEMC=%6.1f PTOT=%6.1f PREC=%6.1f PISO=%6.1f\n",
1967  ele.perp(),tracks.perp(),bemcs.perp(),eemcs.perp(),tot.perp(),rec.perp(),iso.perp());
1968  printf("Dphi=%6.1f ETBal=%6.1f R(ele/iso)=%6.3f\n",dphi,ptbal,riso);
1969 }
1970 
1971 void StFgtAlignmentMaker::fillETBalance(int itrk){
1972  mHit[mQuad][itrk].ele[0]=ele.x();
1973  mHit[mQuad][itrk].ele[1]=ele.y();
1974  mHit[mQuad][itrk].ele[2]=ele.z();
1975  mHit[mQuad][itrk].ele[3]=ele.mag();
1976  mHit[mQuad][itrk].trk[0]=tracks.x();
1977  mHit[mQuad][itrk].trk[1]=tracks.y();
1978  mHit[mQuad][itrk].trk[2]=tracks.z();
1979  mHit[mQuad][itrk].trk[3]=tracks.mag();
1980  mHit[mQuad][itrk].bemc[0]=bemcs.x();
1981  mHit[mQuad][itrk].bemc[1]=bemcs.y();
1982  mHit[mQuad][itrk].bemc[2]=bemcs.z();
1983  mHit[mQuad][itrk].bemc[3]=bemcs.mag();
1984  mHit[mQuad][itrk].eemc[0]=eemcs.x();
1985  mHit[mQuad][itrk].eemc[1]=eemcs.y();
1986  mHit[mQuad][itrk].eemc[2]=eemcs.z();
1987  mHit[mQuad][itrk].eemc[3]=eemcs.mag();
1988  mHit[mQuad][itrk].tot[0]=tot.x();
1989  mHit[mQuad][itrk].tot[1]=tot.y();
1990  mHit[mQuad][itrk].tot[2]=tot.z();
1991  mHit[mQuad][itrk].tot[3]=tot.mag();
1992  mHit[mQuad][itrk].rec[0]=rec.x();
1993  mHit[mQuad][itrk].rec[1]=rec.y();
1994  mHit[mQuad][itrk].rec[2]=rec.z();
1995  mHit[mQuad][itrk].rec[3]=rec.mag();
1996  mHit[mQuad][itrk].iso[0]=iso.x();
1997  mHit[mQuad][itrk].iso[1]=iso.y();
1998  mHit[mQuad][itrk].iso[2]=iso.z();
1999  mHit[mQuad][itrk].iso[3]=iso.mag();
2000  mHit[mQuad][itrk].ptbal=ptbal;
2001  mHit[mQuad][itrk].riso=riso;
2002 }
2003 
2004 void StFgtAlignmentMaker::readFromStEvent() {
2005  mFgtInputRPhi=1;
2006  mEventCounter++; // increase counter
2007  StEvent* event;
2008  event = (StEvent *) GetInputDS("StEvent");
2009  if (!event){
2010  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No StEvent" << endm;
2011  return; // if no event, we're done
2012  }
2013  cout << "Event: Run "<< event->runId() << " Event No: " << event->id() << endl;
2014 
2015  if (event->fgtCollection()) {
2016  cout << "# of FGT hits: " << event->fgtCollection()->getNumHits() << endl;
2017  }else{
2018  cout << "No FGT collection" << endl;
2019  }
2020 
2021  UInt_t NpVX = event->numberOfPrimaryVertices();
2022  if(mDebug>2) cout << "Number of Primary vertex="<<NpVX<<endl;
2023  if (!NpVX) return;
2024 
2025  for (UInt_t i = 0; i < NpVX; i++) {
2026  const StPrimaryVertex *vx = event->primaryVertex(i);
2027  if(mDebug>2) cout << Form("Vertex: %3i ",i) << *vx << endl;
2028  UInt_t nDaughters = vx->numberOfDaughters();
2029  for (UInt_t j = 0; j < nDaughters; j++) {
2030  StPrimaryTrack* pTrack = (StPrimaryTrack*) vx->daughter(j);
2031  if (! pTrack) continue;
2032  //if (pTrack->geometry()->momentum().pseudoRapidity()<1.0) continue;
2033  int nfgt=0, ntpc=0, nprompt=0, quad=-1;
2034  StPtrVecHit &hits=pTrack->detectorInfo()->hits();
2035  for (StPtrVecHitConstIterator iter=hits.begin(); iter != hits.end(); iter++){
2036  StDetectorId det=(*iter)->detector();
2037  if (det == kTpcId) {
2038  ntpc++;
2039  StThreeVectorF xyz=(*iter)->position();
2040  if(xyz.z()>200.0) {nprompt++;}
2041  }else if (det == kFgtId) {
2042  nfgt++;
2043  StFgtPoint* point=(StFgtPoint*)(*iter);
2044  quad=point->getQuad();
2045  }
2046  }
2047  cout << Form("Track: Tpc=%2d Pmp=%1d Fgt=%1d Quad=%1d ",ntpc,nprompt,nfgt,quad) << *pTrack << endl;
2048  if(nfgt==0) continue;
2049  int itrk=mNtrk[quad];
2050  int ihit=0;
2051  //get vertex
2052  mHit[quad][itrk].det[ihit]=24+quad;
2053  StThreeVectorF vxyz=vx->position();
2054  mHit[quad][itrk].x[ihit]=vxyz.x();
2055  mHit[quad][itrk].y[ihit]=vxyz.y();
2056  mHit[quad][itrk].z[ihit]=vxyz.z();
2057  StThreeVectorF evxyz=vx->positionError();
2058  mHit[quad][itrk].ex[ihit]=evxyz.x();
2059  mHit[quad][itrk].ey[ihit]=evxyz.y();
2060  mHit[quad][itrk].ez[ihit]=evxyz.z();
2061  ihit++;
2062  for (StPtrVecHitConstIterator iter=hits.begin(); iter != hits.end(); iter++){
2063  StDetectorId det=(*iter)->detector();
2064  if (det == kTpcId){
2065  mHit[quad][itrk].det[ihit]=28+quad;
2066  StThreeVectorF xyz=(*iter)->position();
2067  mHit[quad][itrk].x[ihit]=xyz.x();
2068  mHit[quad][itrk].y[ihit]=xyz.y();
2069  mHit[quad][itrk].z[ihit]=xyz.z();
2070  if(mHit[quad][itrk].z[ihit]>200.0){mHit[quad][itrk].det[ihit]=32+quad;}
2071  StThreeVectorF exyz=(*iter)->positionError();
2072  mHit[quad][itrk].ex[ihit]=exyz.x();
2073  mHit[quad][itrk].ey[ihit]=exyz.y();
2074  mHit[quad][itrk].ez[ihit]=exyz.z();
2075  }else if (det == kFgtId){
2076  StFgtPoint* point=(StFgtPoint*)(*iter);
2077  int disc=point->getDisc();
2078  int quad2=point->getQuad();
2079  mHit[quad][itrk].det[ihit]=disc*4+quad2;
2080  StThreeVectorF xyz=(*iter)->position();
2081  mHit[quad][itrk].lr[ihit]=point->getPositionR();
2082  mHit[quad][itrk].lp[ihit]=point->getPositionPhi();
2083  mHit[quad][itrk].z[ihit]=StFgtGeom::getDiscZ(disc);
2084  StThreeVectorF exyz=(*iter)->positionError();
2085  if(mErrFgt>=0.0) {mHit[quad][itrk].ex[ihit]=mErrFgt;} else {mHit[quad][itrk].ex[ihit]=exyz.x();}
2086  if(mErrFgt>=0.0) {mHit[quad][itrk].ey[ihit]=mErrFgt;} else {mHit[quad][itrk].ey[ihit]=exyz.y();}
2087  if(mErrFgt>=0.0) {mHit[quad][itrk].ez[ihit]=0.0;} else {mHit[quad][itrk].ez[ihit]=exyz.z();}
2088  }
2089  ihit++;
2090  }
2091  mHit[quad][itrk].nhit=ihit;
2092  for(int i=0; i<mHit[quad][itrk].nhit; i++){
2093  cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
2094  itrk,i,quad,mHit[quad][itrk].det[i],
2095  mHit[quad][itrk].x[i],mHit[quad][itrk].y[i],mHit[quad][itrk].z[i],
2096  mHit[quad][itrk].ex[i],mHit[quad][itrk].ey[i],mHit[quad][itrk].ez[i])
2097  <<endl;
2098  }
2099  mNtrk[quad]++;
2100  }//loop over tracks
2101  }//loop over vertex
2102 }
2103 
2104 void StFgtAlignmentMaker::readFromStEventGlobal() {
2105  mFgtInputRPhi=1;
2106  mEventCounter++; // increase counter
2107  StEvent* event;
2108  event = (StEvent *) GetInputDS("StEvent");
2109  if (!event){
2110  gMessMgr->Warning() << "StFgtAlignmentMaker::Make : No StEvent" << endm;
2111  return; // if no event, we're done
2112  }
2113  cout << "Event: Run "<< event->runId() << " Event No: " << event->id() << endl;
2114 
2115  if (event->fgtCollection()) {
2116  cout << "# of FGT point: " << event->fgtCollection()->getNumPoints() << endl;
2117  }else{
2118  cout << "No FGT collection" << endl;
2119  }
2120 
2121  StSPtrVecTrackNode& trackNode = event->trackNodes();
2122  UInt_t nTracks = trackNode.size();
2123  StTrackNode *node = 0;
2124  cout << "# of tracks "<<nTracks<<endl;
2125  for (UInt_t i=0; i < nTracks; i++) {
2126  node = trackNode[i]; if (!node) continue;
2127  StGlobalTrack* gTrack = static_cast<StGlobalTrack*>(node->track(global));
2128  int nfgt=0, ntpc=0, nprompt=0, quad=-1;
2129  StPtrVecHit &hits=gTrack->detectorInfo()->hits();
2130  for (StPtrVecHitConstIterator iter=hits.begin(); iter != hits.end(); iter++){
2131  StDetectorId det=(*iter)->detector();
2132  if (det == kTpcId) {
2133  ntpc++;
2134  StThreeVectorF xyz=(*iter)->position();
2135  if(xyz.z()>200.0) {nprompt++;}
2136  }else if (det == kFgtId) {
2137  nfgt++;
2138  StFgtPoint* point=(StFgtPoint*)(*iter);
2139  quad=point->getQuad();
2140  }
2141  }
2142  cout << Form("Track: Tpc=%2d Pmp=%1d Fgt=%1d Quad=%1d ",ntpc,nprompt,nfgt,quad) << *gTrack << endl;
2143  if(nfgt==0) continue;
2144  int itrk=mNtrk[quad];
2145  int ihit=0;
2146  for (StPtrVecHitConstIterator iter=hits.begin(); iter != hits.end(); iter++){
2147  StDetectorId det=(*iter)->detector();
2148  if (det == kTpcId){
2149  mHit[quad][itrk].det[ihit]=28+quad;
2150  StThreeVectorF xyz=(*iter)->position();
2151  mHit[quad][itrk].x[ihit]=xyz.x();
2152  mHit[quad][itrk].y[ihit]=xyz.y();
2153  mHit[quad][itrk].z[ihit]=xyz.z();
2154  if(mHit[quad][itrk].z[ihit]>200.0){mHit[quad][itrk].det[ihit]=32+quad;}
2155  StThreeVectorF exyz=(*iter)->positionError();
2156  mHit[quad][itrk].ex[ihit]=exyz.x();
2157  mHit[quad][itrk].ey[ihit]=exyz.y();
2158  mHit[quad][itrk].ez[ihit]=exyz.z();
2159  }else if (det == kFgtId){
2160  StFgtPoint* point=(StFgtPoint*)(*iter);
2161  int disc=point->getDisc();
2162  int quad2=point->getQuad();
2163  mHit[quad][itrk].det[ihit]=disc*4+quad2;
2164  StThreeVectorF xyz=(*iter)->position();
2165  mHit[quad][itrk].lr[ihit]=point->getPositionR();
2166  mHit[quad][itrk].lp[ihit]=point->getPositionPhi();
2167  mHit[quad][itrk].z[ihit]=StFgtGeom::getDiscZ(disc);
2168  StThreeVectorF exyz=(*iter)->positionError();
2169  if(mErrFgt>0.0) {mHit[quad][itrk].ex[ihit]=mErrFgt;} else {mHit[quad][itrk].ex[ihit]=exyz.x();}
2170  if(mErrFgt>0.0) {mHit[quad][itrk].ey[ihit]=mErrFgt;} else {mHit[quad][itrk].ey[ihit]=exyz.y();}
2171  if(mErrFgt>0.0) {mHit[quad][itrk].ez[ihit]=0;} else {mHit[quad][itrk].ez[ihit]=exyz.z();}
2172  }
2173  ihit++;
2174  }
2175  mHit[quad][itrk].nhit=ihit;
2176  for(int i=0; i<mHit[quad][itrk].nhit; i++){
2177  cout<<Form("Trk=%3d Hit=%3d Quad=%1d Det=%2d XYZ=%8.4f %8.4f %8.4f err=%8.4f %8.4f %8.4f",
2178  itrk,i,quad,mHit[quad][itrk].det[i],
2179  mHit[quad][itrk].x[i],mHit[quad][itrk].y[i],mHit[quad][itrk].z[i],
2180  mHit[quad][itrk].ex[i],mHit[quad][itrk].ey[i],mHit[quad][itrk].ez[i])
2181  <<endl;
2182  }
2183  mNtrk[quad]++;
2184  }//loop over tracks
2185 }
2186 
2187 void StFgtAlignmentMaker::writeTree(){
2188  char fname[200];
2189  sprintf(fname,"alignment_trkout_step%d.root",mStep+1);
2190  if(mDay>0){
2191  sprintf(fname,"%d/alignment_trkout_step%d_day%d.root",mDay,mStep+1,mDay);
2192  } else if(mRunNumber>0) {
2193  int yearday=mRunNumber/1000;
2194  sprintf(fname,"%d/alignment_trkout_step%d.%d.root",yearday,mStep+1,mRunNumber);
2195  if(mSeqNumber>0){
2196  sprintf(fname,"%d/alignment_trkout_step%d.%d.%d.root",yearday,mStep+1,mRunNumber,mSeqNumber);
2197  }
2198  }
2199  cout << "Creating "<<fname<<endl;
2200  TFile *f=new TFile(fname,"RECREATE","tracks");
2201  for(int quad=0; quad<kFgtNumQuads; quad++){
2202  char c[10];
2203  sprintf(c,"q%1d",quad);
2204  TTree* t = new TTree(c,c);
2205  t->Branch("run", &(mHitAlg.run), "run/I");
2206  t->Branch("evt", &(mHitAlg.evt), "evt/I");
2207  t->Branch("nhit",&(mHitAlg.nhit), "nhit/I");
2208  t->Branch("nuse",&(mHitAlg.nhitUse), "nuse/I");
2209  t->Branch("nfgt",&(mHitAlg.nhitFgt), "nfgt/I");
2210  t->Branch("nvtx",&(mHitAlg.nhitVtx), "nvtx/I");
2211  t->Branch("ntpc",&(mHitAlg.nhitTpc), "ntpc/I");
2212  t->Branch("nppt",&(mHitAlg.nhitPrompt),"nppt/I");
2213  t->Branch("nemc",&(mHitAlg.nhitEemc), "nemc/I");
2214  t->Branch("used",&(mHitAlg.used), "used/I");
2215  t->Branch("dca" ,&(mHitAlg.dca), "dca/F");
2216  t->Branch("chi2",&(mHitAlg.chi2), "chi2/F");
2217  t->Branch("trkz",&(mHitAlg.trkz), "trkz/F");
2218  t->Branch("dz" , &(mHitAlg.dz), "dz/F");
2219  t->Branch("eta", &(mHitAlg.eta), "eta/F");
2220  t->Branch("phi", &(mHitAlg.phi), "phi/F");
2221  t->Branch("opt", &(mHitAlg.opt), "opt/F");
2222  t->Branch("det", mHitAlg.det, "det[nhit]/I");
2223  t->Branch("x", mHitAlg.x, "x[nhit]/F");
2224  t->Branch("y", mHitAlg.y, "y[nhit]/F");
2225  t->Branch("z", mHitAlg.z, "z[nhit]/F");
2226  t->Branch("ex", mHitAlg.ex, "ex[nhit]/F");
2227  t->Branch("ey", mHitAlg.ey, "ey[nhit]/F");
2228  t->Branch("ez", mHitAlg.ez, "ez[nhit]/F");
2229  t->Branch("lr", mHitAlg.lr, "lr[nhit]/F");
2230  t->Branch("lp", mHitAlg.lp, "lp[nhit]/F");
2231  t->Branch("r", mHitAlg.r, "r[nhit]/F");
2232  t->Branch("p", mHitAlg.p, "p[nhit]/F");
2233  t->Branch("dx", mHitAlg.dx, "dx[nhit]/F");
2234  t->Branch("dy", mHitAlg.dy, "dy[nhit]/F");
2235  t->Branch("dr", mHitAlg.dr, "dr[nhit]/F");
2236  t->Branch("dp", mHitAlg.dp, "dp[nhit]/F");
2237  t->Branch("tw", mHitAlg.tw, "tw[nhit]/F");
2238  t->Branch("p1", mHitAlg.p1, "p1[nhit]/F");
2239  t->Branch("p2", mHitAlg.p2, "p2[nhit]/F");
2240  t->Branch("po", mHitAlg.po, "po[nhit]/F");
2241  t->Branch("su", mHitAlg.su, "su[nhit]/F");
2242  t->Branch("sv", mHitAlg.sv, "sv[nhit]/F");
2243  t->Branch("nrl", mHitAlg.nrl, "nrl[nhit]/I");
2244  t->Branch("nsu", mHitAlg.nsu, "nsu[nhit]/I");
2245  t->Branch("nsv", mHitAlg.nsv, "nsv[nhit]/I");
2246  t->Branch("ele", mHitAlg.ele, "ele[4]/F");
2247  t->Branch("trk", mHitAlg.trk, "trk[4]/F");
2248  t->Branch("bemc",mHitAlg.bemc, "bemc[4]/F");
2249  t->Branch("eemc",mHitAlg.eemc, "eemc[4]/F");
2250  t->Branch("tot", mHitAlg.tot, "tot[4]/F");
2251  t->Branch("rec", mHitAlg.rec, "rec[4]/F");
2252  t->Branch("iso", mHitAlg.iso, "iso[4]/F");
2253  t->Branch("ptbal",&(mHitAlg.ptbal), "ptbal/F");
2254  t->Branch("riso", &(mHitAlg.riso), "riso/F");
2255  t->Branch("use", mHitAlg.use, "use[nhit]/O");
2256 
2257  int itrk=0,ntrk=0;
2258  for(itrk=0; itrk<mNtrk[quad]; itrk++){
2259  if(mHit[quad][itrk].used==0) continue;
2260  memcpy(&mHitAlg,&mHit[quad][itrk],sizeof(mHitAlg));
2261  //printf("WTREE %1d %8d %2d %6.2f %6.2f %6.2f\n",quad,itrk,mHitAlg.nhit,mHitAlg.z[0],mHitAlg.z[1],mHitAlg.z[2]);
2262  t->Fill();
2263  ntrk++;
2264  }
2265  cout<<Form(" Wrote %6d/%6d tracks for quad=%1d",ntrk,itrk,quad)<<endl;
2266  }
2267  f->Write();
2268  f->Close();
2269 }
2270 
2271 void StFgtAlignmentMaker::readFromTree(){
2272  cout <<"ReadFromTree mDb="<<mDb<<endl;
2273  mFgtInputRPhi=1;
2274  if(!mInTreeFile) {
2275  cout<<"StFgtAlignment: Please set input TTree file name using setReadTree()" << endl;
2276  return;
2277  }
2278  cout << "Reading "<<mInTreeFile<<endl;
2279  TFile *f=new TFile(mInTreeFile);
2280  for(int quad=0; quad<kFgtNumQuads; quad++){
2281  char c[10];
2282  sprintf(c,"q%1d",quad);
2283  TTree *t = (TTree*) f->Get(c);
2284  t->SetBranchAddress("run", &(mHitAlg.run));
2285  t->SetBranchAddress("evt", &(mHitAlg.evt));
2286  t->SetBranchAddress("nhit",&(mHitAlg.nhit));
2287  t->SetBranchAddress("nuse",&(mHitAlg.nhitUse));
2288  t->SetBranchAddress("nfgt",&(mHitAlg.nhitFgt));
2289  t->SetBranchAddress("ntpc",&(mHitAlg.nhitTpc));
2290  t->SetBranchAddress("nppt",&(mHitAlg.nhitPrompt));
2291  t->SetBranchAddress("nemc",&(mHitAlg.nhitEemc));
2292  t->SetBranchAddress("used",&(mHitAlg.used));
2293  t->SetBranchAddress("dca" ,&(mHitAlg.dca));
2294  t->SetBranchAddress("chi2",&(mHitAlg.chi2));
2295  t->SetBranchAddress("trkz",&(mHitAlg.trkz));
2296  t->SetBranchAddress("dz" , &(mHitAlg.dz));
2297  t->SetBranchAddress("eta", &(mHitAlg.eta));
2298  t->SetBranchAddress("phi", &(mHitAlg.phi));
2299  t->SetBranchAddress("opt", &(mHitAlg.opt));
2300  t->SetBranchAddress("det", mHitAlg.det);
2301  t->SetBranchAddress("x", mHitAlg.x);
2302  t->SetBranchAddress("y", mHitAlg.y);
2303  t->SetBranchAddress("z", mHitAlg.z);
2304  t->SetBranchAddress("ex", mHitAlg.ex);
2305  t->SetBranchAddress("ey", mHitAlg.ey);
2306  t->SetBranchAddress("ez", mHitAlg.ez);
2307  t->SetBranchAddress("lr", mHitAlg.lr);
2308  t->SetBranchAddress("lp", mHitAlg.lp);
2309  t->SetBranchAddress("r", mHitAlg.r);
2310  t->SetBranchAddress("p", mHitAlg.p);
2311  t->SetBranchAddress("dx", mHitAlg.dx);
2312  t->SetBranchAddress("dy", mHitAlg.dy);
2313  t->SetBranchAddress("dr", mHitAlg.dr);
2314  t->SetBranchAddress("dp", mHitAlg.dp);
2315  t->SetBranchAddress("tw", mHitAlg.tw);
2316  t->SetBranchAddress("p1", mHitAlg.p1);
2317  t->SetBranchAddress("p2", mHitAlg.p2);
2318  t->SetBranchAddress("po", mHitAlg.po);
2319  t->SetBranchAddress("su", mHitAlg.su);
2320  t->SetBranchAddress("sv", mHitAlg.sv);
2321  t->SetBranchAddress("nrl", mHitAlg.nrl);
2322  t->SetBranchAddress("nsu", mHitAlg.nsu);
2323  t->SetBranchAddress("nsv", mHitAlg.nsv);
2324  t->SetBranchAddress("ele", mHitAlg.ele);
2325  t->SetBranchAddress("trk", mHitAlg.trk);
2326  t->SetBranchAddress("bemc",mHitAlg.bemc);
2327  t->SetBranchAddress("eemc",mHitAlg.eemc);
2328  t->SetBranchAddress("tot", mHitAlg.tot);
2329  t->SetBranchAddress("rec", mHitAlg.rec);
2330  t->SetBranchAddress("iso", mHitAlg.iso);
2331  t->SetBranchAddress("ptbal",&(mHitAlg.ptbal));
2332  t->SetBranchAddress("riso", &(mHitAlg.riso));
2333  t->SetBranchAddress("use", mHitAlg.use);
2334  mNtrk[quad]=t->GetEntries();
2335  int itrk;
2336  for(itrk=0; itrk<mNtrk[quad]; itrk++){
2337  if(itrk>=MAXTRK){
2338  cout<<Form(" MAXTRK=%d reached\n",MAXTRK)<<endl;
2339  mNtrk[quad]=MAXTRK;
2340  break;
2341  }
2342  t->GetEntry(itrk);
2343  memcpy(&mHit[quad][itrk],&mHitAlg,sizeof(mHitAlg));
2344  //cout << Form("q=%1d itrk=%6d ",quad,itrk)<<" mDb="<<mDb<<endl;
2345  }
2346  cout<<Form(" Read %6d tracks for quad=%1d",itrk,quad)<<endl;
2347  }
2348  f->Close();
2349 }
2350 
2351 void StFgtAlignmentMaker::setPar(TMinuit* m, fgtAlignment_st* algpar,int discmask,int quadmask, int parmask){
2352  char nam[20];
2353  int iflag;
2354  double *st,*mn,*mx;
2355  double stp[2]={ 0.01, 0.1};
2356  double min[2]={-0.10,-5.0};
2357  double max[2]={ 0.10, 5.0};
2358  double zero[2]={0.00, 0.0};
2359  const char* name[6]={"phi ","theta","psi ","xoff ","yoff ","zoff "};
2360  double *par = &(algpar->phi[0]);
2361  for(int disc=0; disc<6; disc++){
2362  for(int quad=0; quad<4; quad++){
2363  int i=disc*4+quad;
2364  for(int p=0; p<6; p++){
2365  int j= i + p*24;
2366  if( (discmask & 1<<disc) && (quadmask & 1<<quad) && (parmask & 1<<p) ){st=stp; mn=min; mx=max;}
2367  else {st=zero; mn=zero; mx=zero;}
2368  sprintf(nam,"D%1dQ%1d-%5s" ,disc,quad,name[p]);
2369  int k=p/3;
2370  m->mnparm(j, nam, par[j], st[k], mn[k], mx[k], iflag);
2371  }
2372  }
2373  }
2374 }
2375 
2376 void StFgtAlignmentMaker::getPar(TMinuit* m, fgtAlignment_st* algpar){
2377  double r,e;
2378  double *par= &(algpar->phi[0]);
2379  for(int disc=0; disc<6; disc++){
2380  for(int quad=0; quad<4; quad++){
2381  int i=disc*4+quad;
2382  for(int p=0; p<6; p++){
2383  int j= i + p*24;
2384  m->GetParameter(j,r,e);
2385  par[j]=r;
2386  }
2387  }
2388  }
2389 }
2390 
2391 void StFgtAlignmentMaker::writePar(fgtAlignment_st* algpar){
2392  char fname[100]="fgt_alignment.dat";
2393  int yearday=mRunNumber/1000;
2394  if(mRunNumber>0) {
2395  sprintf(fname,"%d/fgt_alignment_%d.dat",yearday,mRunNumber);
2396  if(mSeqNumber>0) sprintf(fname,"%d/fgt_alignment_%d.%07d.dat",yearday,mRunNumber,mSeqNumber);
2397  }
2398  printf("Writing %s\n",fname);
2399  FILE *f=fopen(fname,"w");
2400  for(int disc=0; disc<6; disc++){
2401  for(int quad=0; quad<4; quad++){
2402  int i=disc*4+quad;
2403  printf("%1d %1d %3d %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2404  disc,quad,i,
2405  algpar->phi[i], algpar->theta[i],algpar->psi[i],
2406  algpar->xoff[i],algpar->yoff[i],algpar->zoff[i]);
2407  fprintf(f,"%1d %1d %3d %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2408  disc,quad,i,
2409  algpar->phi[i], algpar->theta[i],algpar->psi[i],
2410  algpar->xoff[i],algpar->yoff[i],algpar->zoff[i]);
2411  }
2412  }
2413  fclose(f);
2414 }
2415 
2416 void StFgtAlignmentMaker::readPar(fgtAlignment_st* algpar){
2417  algpar=orig_algpar=new fgtAlignment_st;
2418  printf("Reading %s\n",mReadParFile);
2419  FILE *f=fopen(mReadParFile,"r");
2420  if(f==0) {printf("Could not open %s\n",mReadParFile); return;}
2421  for(int disc=0; disc<6; disc++){
2422  for(int quad=0; quad<4; quad++){
2423  int ii, i=disc*4+quad;
2424  fscanf(f,"%d %d %d %lf %lf %lf %lf %lf %lf\n",
2425  &disc,&quad,&ii,
2426  &(algpar->phi[i]), &(algpar->theta[i]),&(algpar->psi[i]),
2427  &(algpar->xoff[i]),&(algpar->yoff[i]),&(algpar->zoff[i]));
2428  printf("%1d %1d %3d %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2429  disc,quad,i,
2430  algpar->phi[i], algpar->theta[i],algpar->psi[i],
2431  algpar->xoff[i],algpar->yoff[i],algpar->zoff[i]);
2432  }
2433  }
2434  fclose(f);
2435 }
2436 
2437 void StFgtAlignmentMaker::bookHist(){
2438  const float maxd[NAXIS][kFgtNumDiscs+6]={ {maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdvtx,maxdtpc,maxdtpc,maxdemc,maxdemc,maxdemc},
2439  {maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdvtx,maxdtpc,maxdtpc,maxdemc,maxdemc,maxdemc},
2440  {maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdfgt,maxdvtx,maxdtpc,maxdtpc,maxdemc,maxdemc,maxdemc},
2441  {maxpfgt,maxpfgt,maxpfgt,maxpfgt,maxpfgt,maxpfgt,maxpvtx,maxptpc,maxptpc,maxpemc,maxpemc,maxpemc}};
2442  float xmin[4]={ 0,-10,-40,-40};
2443  float xmax[4]={ 40, 40, 0, 10};
2444  float ymin[4]={-10,-40,-40, 0};
2445  float ymax[4]={ 40, 0, 10, 40};
2446  float rmin=10.0, rmax=40.0;
2447  float phimin[4],phimax[4];
2448  phimin[0]=StFgtGeom::phiQuadXaxis(0); phimax[0]=StFgtGeom::phiQuadXaxis(3);
2449  phimin[1]=StFgtGeom::phiQuadXaxis(1); phimax[1]=StFgtGeom::phiQuadXaxis(0);
2450  phimin[2]=StFgtGeom::phiQuadXaxis(2); phimax[2]=StFgtGeom::phiQuadXaxis(1)+2*PI;
2451  phimin[3]=StFgtGeom::phiQuadXaxis(3); phimax[3]=StFgtGeom::phiQuadXaxis(2);
2452  const int nbin=50;
2453  //1d residuals
2454  const char* caxis1[NAXIS]={"dx","dy","dr","dphi"};
2455  //2d residuals
2456  const char* caxis2[NAXIS*2]={"dx.vs.x","dy.vs.x", "dr.vs.r", "dphi.vs.r",
2457  "dx.vs.y","dy.vs.y", "dr.vs.phi","dphi.vs.phi"};
2458  const char* cquad[kFgtNumQuads]={"A","B","C","D"};
2459  const char* cdisc[kFgtNumDiscs+6]={"1","2","3","4","5","6","V","T","P","E","m","e"};
2460 
2461  for(int disc=0; disc<kFgtNumDiscs+6; disc++){
2462  for(int quad=0; quad<kFgtNumQuads; quad++){
2463  for(int axis=0; axis<NAXIS; axis++){
2464  char c[50];
2465  //1d hist
2466  sprintf(c,"%1s%1s-%s",cdisc[disc],cquad[quad],caxis1[axis]);
2467  hist1[disc][quad][axis]=new TH1F(c,c,nbin,-maxd[axis][disc],maxd[axis][disc]);
2468 
2469  //2d hist
2470  float xmin1,xmax1,xmin2,xmax2;
2471  if(axis<2){
2472  xmin1=xmin[quad]; xmax1=xmax[quad]; xmin2=ymin[quad]; xmax2=ymax[quad];
2473  if(disc==kFgtNumDiscs){ xmin1=-10;xmax1=10; xmin2=-10; xmax2=10;}
2474  if(disc>kFgtNumDiscs) { xmin1=-160; xmax1=160; xmin2=-160; xmax2=160; }
2475  }else{
2476  xmin1=rmin; xmax1=rmax; xmin2=phimin[quad]; xmax2=phimax[quad];
2477  if(disc==kFgtNumDiscs) { xmin1=0; xmax1=10; xmin2=-PI; xmax2=PI;}
2478  if(disc>kFgtNumDiscs) { xmin1*=4; xmax1*=4; xmin2=-PI; xmax2=PI;}
2479  }
2480  sprintf(c,"%1s%1s-%s",cdisc[disc],cquad[quad],caxis2[axis]);
2481  hist2[disc][quad][axis]=new TH2F(c,c,nbin, xmin1,xmax1,nbin,-maxd[axis][disc],maxd[axis][disc]);
2482  sprintf(c,"%1s%1s-%s",cdisc[disc],cquad[quad],caxis2[axis+NAXIS]);
2483  hist2[disc][quad][axis+NAXIS]=new TH2F(c,c,nbin, xmin2,xmax2,nbin,-maxd[axis][disc],maxd[axis][disc]);
2484  }
2485  }
2486  }
2487  histdz=new TH1F("Zfgt-Ztpc","Zfgt-Ztpc",60,-30,30);
2488 }
2489 
2490 void StFgtAlignmentMaker::resetHist(){
2491  for(int disc=0; disc<kFgtNumDiscs+6; disc++){
2492  for(int quad=0; quad<kFgtNumQuads; quad++){
2493  for(int axis=0; axis<NAXIS; axis++){
2494  hist1[disc][quad][axis]->Reset("ICES");
2495  hist2[disc][quad][axis]->Reset("ICES");
2496  hist2[disc][quad][axis+NAXIS]->Reset("ICES");
2497  }
2498  }
2499  }
2500 }
2501 
2502 void StFgtAlignmentMaker::saveHist(){
2503  char fname[100];
2504  sprintf(fname,"alignment_step%d.root",mStep+1);
2505  if(mDay>0){
2506  sprintf(fname,"%d/alignment_step%d_day%d.root",mDay,mStep+1,mDay);
2507  } else if(mRunNumber>0) {
2508  int yearday=mRunNumber/1000;
2509  sprintf(fname,"%d/alignment_step%d.%d.root",yearday,mStep+1,mRunNumber);
2510  if(mSeqNumber>0){
2511  sprintf(fname,"%d/alignment_step%d.%d.%d.root",yearday,mStep+1,mRunNumber,mSeqNumber);
2512  }
2513  }
2514  cout << "Writing " << fname << endl;
2515  TFile *hfile = new TFile(fname,"RECREATE");
2516  for(int disc=0; disc<kFgtNumDiscs+6; disc++){
2517  for(int quad=0; quad<kFgtNumQuads; quad++){
2518  for(int axis=0; axis<NAXIS; axis++){
2519  hist1[disc][quad][axis]->Write();
2520  hist2[disc][quad][axis]->Write();
2521  hist2[disc][quad][axis+NAXIS]->Write();
2522  }
2523  }
2524  }
2525  histdz->Write();
2526  hfile->Close();
2527 }
2528 
2529 void StFgtAlignmentMaker::setHitMask(int hitmask_disc){
2530  cout << Form("Alignment Step %d for Quad=%1d using Hits from=",mStep,mQuad);
2531  for(int i=0; i<6; i++) {if(hitmask_disc & (1<<i) ) {cout << Form("FgtD%1d ",i+1);}}
2532  if(hitmask_disc & 0x40) cout << "Vertex ";
2533  if(hitmask_disc & 0x80) cout << "TPC ";
2534  if(hitmask_disc & 0x100) cout << "Prompt ";
2535  if(hitmask_disc & 0x200) cout << "Eemc ";
2536  if(hitmask_disc & 0x400) cout << "Mip ";
2537  if(hitmask_disc & 0x800) cout << "Ele ";
2538  cout << Form("(hitmask_disc = %03x)",hitmask_disc);
2539  mHitMaskDisc=hitmask_disc;
2540  mNtrkUse[mQuad]=0;
2541  for(int itrk=0; itrk<mNtrk[mQuad]; itrk++){
2542  mHit[mQuad][itrk].nhitUse=0;
2543  mHit[mQuad][itrk].nhitFgt=0;
2544  mHit[mQuad][itrk].nhitVtx=0;
2545  mHit[mQuad][itrk].nhitTpc=0;
2546  mHit[mQuad][itrk].nhitPrompt=0;
2547  mHit[mQuad][itrk].nhitEemc=0;
2548  mHit[mQuad][itrk].used=0;
2549  for(int ihit=0; ihit<mHit[mQuad][itrk].nhit; ihit++){
2550  mHit[mQuad][itrk].use[ihit]=false;
2551  int det=mHit[mQuad][itrk].det[ihit];
2552  int disc=det/4; //fgt disc0-5
2553  if((hitmask_disc & (1<<disc)) ==0) continue;
2554  if (disc<6) {
2555  if(mFgtRCut[mStep]>0 && fabs(mHit[mQuad][itrk].dr[ihit])>mFgtRCut[mStep]) continue;
2556  if(mFgtPCut[mStep]>0 && fabs(mHit[mQuad][itrk].dp[ihit])>mFgtPCut[mStep]) continue;
2557  mHit[mQuad][itrk].nhitFgt++;
2558  }else if(disc==6) { //vertex
2559  if(mDzCut[mStep]>0 && fabs(mHit[mQuad][itrk].dz )>mDzCut[mStep] ) continue;
2560  if(mDcaCut[mStep]>0 && fabs(mHit[mQuad][itrk].dca)>mDcaCut[mStep]) continue;
2561  mHit[mQuad][itrk].nhitVtx++;
2562  }else if(disc==7) { //tpc
2563  if(mTpcRCut[mStep]>0 && fabs(mHit[mQuad][itrk].dr[ihit])>mTpcRCut[mStep]) continue;
2564  if(mTpcPCut[mStep]>0 && fabs(mHit[mQuad][itrk].dp[ihit])>mTpcPCut[mStep]) continue;
2565  mHit[mQuad][itrk].nhitTpc++;
2566  }else if(disc==8) { //tpc prompt
2567  if(mTpcRCut[mStep]>0 && fabs(mHit[mQuad][itrk].dr[ihit])>mTpcRCut[mStep]) continue;
2568  if(mTpcPCut[mStep]>0 && fabs(mHit[mQuad][itrk].dp[ihit])>mTpcPCut[mStep]) continue;
2569  mHit[mQuad][itrk].nhitTpc++;
2570  mHit[mQuad][itrk].nhitPrompt++;
2571  }else if(disc==9 || disc==10 || disc==11) { //eemc, mip, electron
2572  if(mEmcRCut[mStep]>0 && fabs(mHit[mQuad][itrk].dr[ihit])>mEmcRCut[mStep]) continue;
2573  if(mEmcPCut[mStep]>0 && fabs(mHit[mQuad][itrk].dp[ihit])>mEmcPCut[mStep]) continue;
2574  mHit[mQuad][itrk].nhitEemc++;
2575  //if(mHit[mQuad][itrk].dr[itrk]!=0)
2576  // printf("EEMCCUT %10.6f %10.6f %d\n",mHit[mQuad][itrk].dr[ihit],mEmcRCut[mStep],mHit[mQuad][itrk].nhitEemc);
2577  }
2578  mHit[mQuad][itrk].use[ihit]=true;
2579  mHit[mQuad][itrk].nhitUse++;
2580  }
2581  //printf("SetHitMask %d n=%d %d %d %d %d %f %f\n",itrk,mHit[mQuad][itrk].nhitUse,mHit[mQuad][itrk].nhitFgt,mHit[mQuad][itrk].nhitTpc,mHit[mQuad][itrk].nhitPrompt,mHit[mQuad][itrk].nhitEemc,mHit[mQuad][itrk].dca,mHit[mQuad][itrk].chi2);
2582  if(mHit[mQuad][itrk].nhitUse>=mReqHit &&
2583  mHit[mQuad][itrk].nhitFgt>=mReqFgtHit &&
2584  mHit[mQuad][itrk].nhitVtx>=mReqVtx &&
2585  mHit[mQuad][itrk].nhitTpc>=mReqTpcHit &&
2586  mHit[mQuad][itrk].nhitPrompt>=mReqPromptHit &&
2587  mHit[mQuad][itrk].nhitEemc>=mReqEemcHit){
2588  mHit[mQuad][itrk].used=1;
2589  mNtrkUse[mQuad]++;
2590  }
2591  }
2592  cout << Form(" usable ntrack=%d with nhit>=%d nFgt>=%d vVtx>=%d nTpc>=%d nPrompt>=%d nEemc>=%d",
2593  mNtrkUse[mQuad],mReqHit,mReqFgtHit,mReqVtx,mReqTpcHit,mReqPromptHit,mReqEemcHit)<<endl;
2594 }
2595 
2596 void StFgtAlignmentMaker::overWriteError(){
2597  printf("Over Writing Hit Errors\n");
2598  for(int iquad=0; iquad<kFgtNumQuads; iquad++){
2599  printf(" quad=%d ntrk=%d\n",iquad,mNtrk[iquad]);
2600  for(int itrk=0; itrk<mNtrk[iquad]; itrk++){
2601  for(int ihit=0; ihit<mHit[iquad][itrk].nhit; ihit++){
2602  // printf("BBB quad=%d itrk=%d nhit=%d ihit=%d 1928=%d\n",iquad,itrk,mHit[iquad][itrk].nhit,ihit,mHit[iquad][1928].nhit);
2603  if(ihit>100) return; //hack
2604  int det=mHit[iquad][itrk].det[ihit];
2605  int disc=det/4;
2606  if (disc<6) { //fgt disc0-5
2607  mHit[iquad][itrk].ex[ihit]=mErrFgt;
2608  mHit[iquad][itrk].ey[ihit]=mErrFgt;
2609  mHit[iquad][itrk].ez[ihit]=0.0;
2610  }else if(disc==6) { //vertex
2611  mHit[iquad][itrk].ex[ihit]=mErrVtx;
2612  mHit[iquad][itrk].ey[ihit]=mErrVtx;
2613  mHit[iquad][itrk].ez[ihit]=mErrVtxZ;
2614  }else if(disc==7) { //tpc
2615  float x= mHit[iquad][itrk].x[ihit];
2616  float y= mHit[iquad][itrk].y[ihit];
2617  float r=sqrt(x*x+y*y);
2618  if(r<127.950){ //inner
2619  mHit[iquad][itrk].ex[ihit]=mErrTpcI;
2620  mHit[iquad][itrk].ey[ihit]=mErrTpcI;
2621  mHit[iquad][itrk].ez[ihit]=mErrTpcZ;
2622  }else{ //outer
2623  mHit[iquad][itrk].ex[ihit]=mErrTpcO;
2624  mHit[iquad][itrk].ey[ihit]=mErrTpcO;
2625  mHit[iquad][itrk].ez[ihit]=mErrTpcZ;
2626  }
2627  }else if(disc==8) { //tpc prompt
2628  mHit[iquad][itrk].ex[ihit]=mErrPpt;
2629  mHit[iquad][itrk].ey[ihit]=mErrPpt;
2630  mHit[iquad][itrk].ez[ihit]=0;
2631  }else if(disc==9 || disc==10 || disc==11) { //eemc, mip, electron
2632  mHit[iquad][itrk].ex[ihit]=mErrEmc;
2633  mHit[iquad][itrk].ey[ihit]=mErrEmc;
2634  mHit[iquad][itrk].ez[ihit]=0;
2635  }
2636  }
2637  }
2638  }
2639 }
2640 
2641 void StFgtAlignmentMaker::doAlignment(fgtAlignment_st* input,
2642  int discmask,int quadmask, int parmask, int hitmask_disc, int residmask,
2643  int trackType, int minHit, int minFgtHit, int minVtx, int minTpcHit, int minPromptHit, int minEemcHit,
2644  fgtAlignment_st* result){
2645  double arg[10];
2646  int iflag, min=0;
2647 
2648  TMinuit *m = new TMinuit(NPAR);
2649  arg[0] =-1; m->mnexcm("SET PRI", arg ,1,iflag);
2650  arg[0] =-1; m->mnexcm("SET NOWarning", arg, 0,iflag);
2651  arg[0] = 1; m->mnexcm("SET ERR", arg ,1,iflag);
2652  if(trackType==0) { m->SetFCN(funcLine); min=2;}
2653  else if(trackType==1) { m->SetFCN(funcHelix); min=3;}
2654  mReqHit=minHit;
2655  if(mReqHit<min) mReqHit=min;
2656  mReqFgtHit=minFgtHit;
2657  mReqVtx=minVtx;
2658  mReqTpcHit=minTpcHit;
2659  mReqPromptHit=minPromptHit;
2660  mReqEemcHit=minEemcHit;
2661  mResidMaskDisc=residmask;
2662 
2663  setHitMask(hitmask_disc);
2664  if(mNtrkUse[mQuad]==0) return;
2665  setPar(m,input,discmask,quadmask,parmask);
2666 
2667  if(parmask>0) {
2668  arg[0] = 0; m->mnexcm("SET PRI", arg ,1,iflag);
2669  //arg[0] =-1; m->mnexcm("SET NOWarning", arg, 0,iflag);
2670  }
2671 
2672  if(mFillHist>0){
2673  int np=NPAR;
2674  double *gin=0, f=0;
2675  double *par = (double *) input;
2676  if (trackType==0) { funcLine (np,gin,f,par,iflag);}
2677  else if(trackType==1) { funcHelix(np,gin,f,par,iflag);}
2678  }else{
2679  arg[0] = 20000; arg[1] = 1.; m->mnexcm("MIGRAD", arg ,2, iflag);
2680  if(mFillHist==0) getPar(m,result);
2681  }
2682 }
2683 
StMuDst::primaryVertex
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
Line
STAR includes.
Definition: Line.hh:23
StHelix::x
double x(double s) const
coordinates of helix at point s
Definition: StHelix.hh:182
StEEmcIUPoint::energy
void energy(Float_t e, Int_t layer=0)
Set the energy of this point.
Definition: StEEmcIUPoint.h:25
StEEmcIUPoint::position
void position(TVector3 p)
Set the position of this point at the SMD plane.
Definition: StEEmcIUPoint.h:23
StEEmcIUPointMaker::points
StEEmcIUPointVec_t points()
Return a unique key assigned by the cluster maker.
Definition: StEEmcIUPointMaker.h:182
StEEmcIUPointMaker
Class for building points from smd clusters.
Definition: StEEmcIUPointMaker.h:21
StEEmcIUPointMaker::numberOfPoints
Int_t numberOfPoints()
Return number of points.
Definition: StEEmcIUPointMaker.h:183
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StEEmcIUPoint::cluster
void cluster(StEEmcIUSmdCluster c, Int_t plane)
Add an smd cluster to this point.
Definition: StEEmcIUPoint.h:31
StFgtAlignmentMaker
A typical Analysis Class.
Definition: StFgtAlignmentMaker.h:68
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StEEmcIUPointMaker::point
StEEmcIUPoint point(Int_t ipoint)
Return a specified point.
Definition: StEEmcIUPointMaker.h:184
StMaker::GetName
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
kStOK
Definition: Stypes.h:40
StEEmcIUPoint
Base class for representing EEMC points.
Definition: StEEmcIUPoint.h:12
St_vertexSeedC
Definition: St_vertexSeedC.h:7
StEEmcIUPoint::numberOfRelatives
void numberOfRelatives(Int_t r)
Set the number of other points which share tower energy.
Definition: StEEmcIUPoint.h:34
StFgtPoint
Represents a point in the FGT.
Definition: StFgtPoint.h:49
Generated by   doxygen 1.8.5