StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StLaserAnalysisMaker.cxx
1 //
2 // $Id: StLaserAnalysisMaker.cxx,v 1.23 2015/02/10 20:27:16 fisyak Exp $
3 // $Log: StLaserAnalysisMaker.cxx,v $
4 // Revision 1.23 2015/02/10 20:27:16 fisyak
5 // Adjust split style for ROOT_VERSION_CODE
6 //
7 // Revision 1.22 2015/02/06 00:31:46 perev
8 // defence agains zero dcaGeometry pointer added
9 //
10 // Revision 1.21 2014/09/17 19:11:49 fisyak
11 // Fix bug in Fit logic, bug #2901
12 //
13 // Revision 1.20 2014/08/19 17:23:31 fisyak
14 // Activate CORRECT_RAFT_DIRECTION
15 //
16 // Revision 1.19 2014/03/14 13:11:20 fisyak
17 // comment out gStTpcDb->ScaleY()
18 //
19 // Revision 1.18 2014/03/13 21:59:45 fisyak
20 // add cluster position in Local Sector Coordinate System
21 //
22 // Revision 1.16 2009/03/11 15:49:40 fisyak
23 // Remove hits not beloging to primary tracks
24 //
25 // Revision 1.15 2008/06/02 13:48:03 fisyak
26 // Add t0 handlers for Tpx/Tpc time offsets
27 //
28 // Revision 1.14 2007/10/16 15:26:03 fisyak
29 // Freeze version for Run VII
30 //
31 // Revision 1.13 2007/07/05 14:37:04 fisyak
32 // Freeze version for run VII
33 //
34 // Revision 1.12 2007/05/09 13:36:44 fisyak
35 // Freeze before intruducing Bundle coordinate system
36 //
37 // Revision 1.11 2007/04/17 05:12:06 perev
38 // GetTFile()==>StMaker. Jerome request
39 //
40 // Revision 1.10 2007/03/06 16:31:55 fisyak
41 // Before Selection of good runs
42 //
43 //#define CORRECT_LASER_POSITIONS
44 #define CORRECT_RAFT_DIRECTION
45 #define __TRACKHITS__
46 #ifndef __TRACKHITS__
47 #define ADDPRIMTRACKHITS
48 #endif
49 #include <assert.h>
50 #include "TFile.h"
51 #include "StEventTypes.h"
52 #include "StLaserAnalysisMaker.h"
53 #include "laserino.h"
54 #include "LaserEvent.h"
55 #include "StTpcDb/StTpcDb.h"
56 #include "TGeoMatrix.h"
57 #include "TRVector.h"
58 #include "TRMatrix.h"
59 #include "TRSymMatrix.h"
60 #include "StDetectorDbMaker/StDetectorDbClock.h"
61 #define PrPP(A,B) if (Debug()) {LOG_INFO << "::StLaserAnalysisMaker" << (#A) << "\t" << (#B) << " = \t" << (B) << endm;}
62 ClassImp(StLaserAnalysisMaker);
63 static LaserEvent *event = 0;
64 static const Int_t NS = 12;
65 static const Int_t NB = 6;
66 static const Int_t NM = 7;
67 static LaserB *Lasers[NS][NB][NM];
68 #if 1
69 static Int_t NoBeams = 0;
70 static LaserRaft *LaserBeams[NS*NB*NM];
71 #endif
72 static TGeoHMatrix *Traft[14];
73 static TGeoHMatrix *Raft2Tpc[14];
74 static TGeoHMatrix *Bundles2Tpc[14][6];
75 static TGeoHMatrix *Mirrors2Tpc[14][6][7];
76 //________________________________________________________________________________
77 Int_t StLaserAnalysisMaker::Init(){
78  if (! IsActive()) return kStOk;
79  event = new LaserEvent();
80  TFile *f = GetTFile();
81 #if 0
82  if (! f) {
83  f = new TFile("StLaserAnalysisMaker.root","recreate");
84  }
85 #endif
86  assert(f);
87  if (f) {
88  f->cd();
89  m_laser = new TTree("laser","Tpc laser track tree");
90  m_laser->SetAutoSave(100000000); //Save every 100 MB
91  Int_t bufsize= 64000;
92 #if ROOT_VERSION_CODE <= ROOT_VERSION(5,34,10)
93  Int_t split = 99;
94 #else
95  // Int_t split = -2; // by default, split Event in sub branches << old style
96  Int_t split = 99;
97 #endif
98  if (split) bufsize /= 4;
99  Int_t branchStyle = 1; //new style by default
100  if (split < 0) {branchStyle = 0; split = -1-split;}
101  TTree::SetBranchStyle(branchStyle);
102  m_laser->Branch("event", "LaserEvent",&event, bufsize, split);
103  }
104  return StMaker::Init();
105 }
106 //_____________________________________________________________________________
107 Int_t StLaserAnalysisMaker::InitRun(Int_t run){
108  const Int_t numSectors = 24;
109  Double_t gamma = 0;
110  Double_t dGamma = 720./numSectors;
111  Int_t sector = 0;
112  TGeoHMatrix TpcHalf[2];
113  Double_t rotHalfs[18] = {
114  0, 1, 0, -1, 0, 0, 0, 0, 1, // sector 13-24
115  0, -1, 0, -1, 0, 0, 0, 0, -1 // sector 1-12
116  };
117  for (Int_t half = east; half <= west; half++) TpcHalf[half].SetRotation(&rotHalfs[9*half]);
118  TGeoHMatrix RotSec[24];
119  for (sector = 1; sector <= numSectors; sector++) {
120  if (sector > 12) gamma = (numSectors-sector)*dGamma;
121  else gamma = sector*dGamma;
122  RotSec[sector-1].RotateZ(-gamma);
123  }
124 #ifdef CORRECT_RAFT_DIRECTION
125  Double_t zWheel = (229.71+1.7780);
126  // East West
127  Double_t ZWheel[2] = {-zWheel, zWheel};
128  Double_t RDWheel[2] = { 190.5802, 190.5705};
129  StThreeVectorD xyzDSE(RDWheel[0], 0, ZWheel[0]);
130  StThreeVectorD xyzDSW(RDWheel[1], 0, ZWheel[1]);
131  TGeoHMatrix R;
132  StThreeVectorD xyzTE, xyzTW;
133  if (Debug()) {
134  cout << " TpcHalf(east) "; StTpcDb::instance()->TpcHalf(east).Print();
135  cout << " TpcHalf(west) "; StTpcDb::instance()->TpcHalf(west).Print();
136  }
137  for (sector = 1; sector <= 12; sector++) {
138  if (Debug()) {
139  cout << "Sector " << sector << " ===========" << endl;
140  }
141  // S2R = Rot^-1 * Half * Rot
142  TGeoHMatrix ET = RotSec[sector+12-1].Inverse()*StTpcDb::instance()->TpcHalf(east)*RotSec[sector+12-1]; if (Debug()) ET.Print();
143  ET.LocalToMaster(xyzDSE.xyz(), xyzTE.xyz()); PrPP(InitRun,xyzTE);
144  TGeoHMatrix EW = RotSec[sector -1].Inverse()*StTpcDb::instance()->TpcHalf(west)*RotSec[sector -1]; if (Debug()) EW.Print();
145  EW.LocalToMaster(xyzDSW.xyz(), xyzTW.xyz()); PrPP(InitRun,xyzTW);
146  // Survey line direction
147  StThreeVectorD dif = (xyzTW - xyzTE)/2.; PrPP(InitRun,dif);
148  // Survey line center
149  StThreeVectorD sum = (xyzTW + xyzTE)/2.; PrPP(InitRun,sum);
150  StThreeVectorD unit = dif.unit(); PrPP(InitRun,unit);
151  Double_t alpha = - unit.y();
152  Double_t beta = unit.x();
153  StThreeVectorD tra = sum; tra.xyz()[0] -= 0.5*(RDWheel[0]+RDWheel[1]); PrPP(InitRun,tra);
154  // Transformation from survey line to raft
155  TGeoHMatrix S2R(Form("S2R_%i",sector));
156  S2R.RotateX(alpha*180/TMath::Pi());
157  S2R.RotateY(beta*180/TMath::Pi());
158  S2R.SetTranslation(tra.xyz()); if (Debug()) {cout << "S2R "; S2R.Print();}
159  StThreeVectorD unitR;
160  S2R.MasterToLocalVect(unit.xyz(),unitR.xyz()); PrPP(InitRun,unitR);
161  if (Debug()) {cout << "RotSecO[" << sector-1 << "]:"; RotSec[sector-1].Print();}
162  R = RotSec[sector-1]*S2R; if (Debug()) {cout << "R:"; R.Print();}
163  RotSec[sector-1] = R; if (Debug()) {cout << "RotSecM[" << sector-1 << "]:"; RotSec[sector-1].Print();}
164  if (Debug()) {cout << "RotSecO[" << sector+12-1 << "]:"; RotSec[sector+12-1].Print();}
165  // TGeoHMatrix R = RotSec[sector+12-1]*S2R.Inverse(); if (Debug()) {cout << "R:"; R.Print();}
166  TGeoHMatrix R = RotSec[sector+12-1]*S2R; if (Debug()) {cout << "R:"; R.Print();}
167  RotSec[sector+12-1] = R; if (Debug()) {cout << "RotSecM[" << sector+12-1 << "]:"; RotSec[sector+12-1].Print();}
168  }
169 #endif /* CORRECT_RAFT_DIRECTION */
170  memset (LaserBeams, 0, NS*NB*NM*sizeof(LaserRaft*));
171  NoBeams = 0;
172  memset(Traft, 0, 14*sizeof(TGeoHMatrix *));
173  memset(Raft2Tpc, 0, 14*sizeof(TGeoHMatrix *));
174  memset(Bundles2Tpc, 0, 14*6*sizeof(TGeoHMatrix *));
175  memset(Mirrors2Tpc, 0, 14*6*7*sizeof(TGeoHMatrix *));
176  Double_t y0[12] = { 0.0373, -0.0104, -0.0081, -0.0092, 0.0000, 0.0492, 0.0008, -0.0123, 0.0281, 0.0210, -0.0102, -0.0627};
177  Double_t y1[12] = { 0.0088, -0.0033, 0.0000, -0.0045, 0.0000, 0.0079, 0.0006, -0.0013, 0.0068, 0.0052, -0.0033, -0.0168};
178  Double_t z0[12] = {-0.0414, 0.0363, -0.1394, 0.0508, 0.0000, 0.0241, -0.0331, 0.0689, -0.1474, -0.0469, 0.1104, 0.0203};
179  Double_t z1[12] = {-0.0002, -0.0001, -0.0089, -0.0002, 0.0000, -0.0000, -0.0002, -0.0000, 0.0003, 0.0031, 0.0007, 0.0002};
180  for (Int_t i = 0; i < NoRaftPositions; i++) {
181  if (! RaftPositions[i].Sector) continue;
182  Int_t raft = RaftPositions[i].Raft;
183  Traft[raft-1] = new TGeoHMatrix(Form("Raft%0i",raft));
184  Traft[raft-1]->RotateX(RaftPositions[i].alpha*180/TMath::Pi());
185  Traft[raft-1]->RotateY(RaftPositions[i].beta*180/TMath::Pi());
186  Traft[raft-1]->RotateZ(RaftPositions[i].gamma*180/TMath::Pi());
187  StThreeVectorD xyzRaft(RaftPositions[i].X,RaftPositions[i].Y,RaftPositions[i].Z - 0.05465 + 0.1022-0.0530);
188  if (RaftPositions[i].Sector <= 12) xyzRaft.setY(xyzRaft.y() - 0.0480 - 0.0095 - 0.0028);
189  else xyzRaft.setY(xyzRaft.y() + 0.0328 + 0.0118 + 0.0032);
190  Int_t s = (RaftPositions[i].Sector-1)/2;
191  xyzRaft.setY(xyzRaft.y() + y0[s] + y1[s]);
192  xyzRaft.setZ(xyzRaft.z() + z0[s] + z1[s]);
193  Traft[raft-1]->SetTranslation(xyzRaft.xyz());
194  if (Debug()) {
195  RaftPositions[i].Print();
196  Traft[raft-1]->Print();
197  }
198  }
199  for (Int_t r = 0; r < 14; r++) {
200  Int_t sector = Bundles[r][0].Sector;
201  if (! sector) continue;
202  Int_t raft = Bundles[r][0].Raft;
203  Int_t half = west;
204  if (sector > 12) half = east;
205  Raft2Tpc[raft-1] = new TGeoHMatrix(Form("Raft%iToTpc",raft));
206  *Raft2Tpc[raft-1] = RotSec[sector-1] * TpcHalf[half] * (*Traft[raft-1]);
207  for (Int_t b = 0; b < 6; b++) {
208  Int_t bundle = Bundles[r][b].Bundle;
209  TGeoHMatrix rotB;
210  rotB.SetTranslation(&Bundles[r][b].X);
211  Bundles2Tpc[raft-1][bundle-1] = new TGeoHMatrix(Form("BundleR%iB%i",raft,bundle));
212  *Bundles2Tpc[raft-1][bundle-1] = *Raft2Tpc[raft-1] * rotB;
213  for (Int_t m = 0; m < 7; m++) {
214  if (Mirrors[r][b][m].Sector < 2) continue;
215  LASERINO_t *local = &Mirrors[r][b][m];
216  if (! local) continue;
217  Int_t mirror = Mirrors[r][b][m].Mirror;
218  LaserBeams[NoBeams] = new LaserRaft();
219  LaserBeams[NoBeams]->Sector = local->Sector;
220  LaserBeams[NoBeams]->Raft = local->Raft;
221  LaserBeams[NoBeams]->Bundle = local->Bundle;
222  LaserBeams[NoBeams]->Mirror = local->Mirror;
223 #ifdef CORRECT_LASER_POSITIONS
224  LaserBeams[NoBeams]->XyzB = StThreeVectorD(Mirrors[r][b][m].X+Mirrors[r][b][m].dX,
225  Mirrors[r][b][m].Y+Mirrors[r][b][m].dY,
226  Mirrors[r][b][m].Z+Mirrors[r][b][m].dZ);
227 #else
228  LaserBeams[NoBeams]->XyzB = StThreeVectorD(Mirrors[r][b][m].X,
229  Mirrors[r][b][m].Y,
230  Mirrors[r][b][m].Z);
231 #endif
232  Double_t theta = Bundles[r][b].ThetaZ + Mirrors[r][b][m].ThetaZ;
233  Double_t phi = Bundles[r][b].Phi + Mirrors[r][b][m].Phi;
234  TGeoHMatrix rotM;
235  rotM.SetTranslation(LaserBeams[NoBeams]->XyzB.xyz());
236  Mirrors2Tpc[raft-1][bundle-1][mirror-1] = new TGeoHMatrix(Form("MirrorR%iB%iM%i",raft,bundle,mirror));
237  *Mirrors2Tpc[raft-1][bundle-1][mirror-1] = *Bundles2Tpc[raft-1][bundle-1] * rotM;
238 
239  LaserBeams[NoBeams]->dirB = StThreeVectorD(-TMath::Cos(phi)*TMath::Cos(theta),
240  -TMath::Sin(phi)*TMath::Cos(theta),
241  -TMath::Sin(theta));
242  rotB.LocalToMaster(LaserBeams[NoBeams]->XyzB.xyz(), LaserBeams[NoBeams]->XyzU.xyz());
243  rotB.LocalToMasterVect(LaserBeams[NoBeams]->dirB.xyz(), LaserBeams[NoBeams]->dirU.xyz());
244  Raft2Tpc[raft-1]->LocalToMaster(LaserBeams[NoBeams]->XyzU.xyz(),LaserBeams[NoBeams]->XyzL.xyz());
245  Raft2Tpc[raft-1]->LocalToMasterVect(LaserBeams[NoBeams]->dirU.xyz(),LaserBeams[NoBeams]->dirL.xyz());
246  LaserBeams[NoBeams]->Theta = LaserBeams[NoBeams]->dirL.theta();
247  LaserBeams[NoBeams]->Phi = LaserBeams[NoBeams]->dirL.phi();
248  if (Debug()) {
249  local->Print();
250  Raft2Tpc[raft-1]->Print();
251  LaserBeams[NoBeams]->Print();
252  }
253  NoBeams++;
254  }
255  }
256  }
257  // average Z for membrane = -3.6 cm
258  memset(Lasers, 0, NS*NB*NM*sizeof(LaserB *));
259  const TGeoHMatrix &Tpc2Global = gStTpcDb->Tpc2GlobalMatrix();
260  LaserB *theLaser = 0;
261  for (Int_t i = 0; i < NoBeams; i++) {
262  if (! LaserBeams[i]) continue;
263  Int_t s = LaserBeams[i]->Sector/2 - 1;
264  if (s < 0 || s >= NS) continue;
265  Int_t b = LaserBeams[i]->Bundle - 1;
266  if (b < 0 || b >= NB) continue;
267  Int_t m = LaserBeams[i]->Mirror - 1;
268  if (m < 0 || m >= NM) continue;
269  theLaser = new LaserB(*LaserBeams[i]);
270  Lasers[s][b][m] = theLaser;
271  Tpc2Global.LocalToMaster(theLaser->XyzL.xyz(),theLaser->XyzG.xyz());
272  Tpc2Global.LocalToMasterVect(theLaser->dirL.xyz(),theLaser->dirG.xyz());
273  gStTpcDb->SupS2Tpc(theLaser->Sector).MasterToLocal(theLaser->XyzL.xyz(),theLaser->XyzS.xyz());
274  gStTpcDb->SupS2Tpc(theLaser->Sector).MasterToLocalVect(theLaser->dirL.xyz(),theLaser->dirS.xyz());
275  theLaser->PhiG = theLaser->dirG.phi();
276  theLaser->ThetaG = theLaser->dirG.theta();
277  theLaser->IsValid = 0;
278 #if 0
279  // From Laser Z
280  if (theLaser->Sector == 2 && theLaser->Bundle == 3 && theLaser->Mirror == 5) continue; // misplaced
281  if (theLaser->Sector == 4 && theLaser->Bundle == 4) continue; //dead
282  if (theLaser->Sector == 8) continue; // dead
283  if (theLaser->Sector == 10) continue; // dead
284  if (theLaser->Sector == 12 && theLaser->Bundle == 4 && theLaser->Mirror == 4) continue; // dead
285  if (theLaser->Sector == 12 && theLaser->Bundle == 4 && theLaser->Mirror == 5) continue; // dead
286  if (theLaser->Sector == 14 && theLaser->Bundle == 4 && theLaser->Mirror == 4) continue; // misplaced swap ?
287  if (theLaser->Sector == 16 && theLaser->Bundle == 4 && theLaser->Mirror == 5) continue; // misplaced
288  if (theLaser->Sector == 16 && theLaser->Bundle == 4 && theLaser->Mirror == 5) continue; // misplaced
289  if (theLaser->Sector == 18 && theLaser->Bundle == 2) continue; // very strange pattern
290  if (theLaser->Sector == 22 && theLaser->Bundle == 3) continue; // strange pattern
291  if (theLaser->Sector == 22 && theLaser->Bundle == 4) continue; // missing
292  if (theLaser->Sector == 24 && theLaser->Bundle == 2 && theLaser->Mirror == 6) continue; // missing
293 #endif
294  theLaser->IsValid = 1;
295  if (Debug()) {
296  LaserBeams[i]->Print();
297  theLaser->Print();
298  }
299  }
300  return kStOk;
301 }
302 //_____________________________________________________________________________
303 void StLaserAnalysisMaker::Clear(const Option_t *option){
304  if (event) event->Clear("C");
305  StMaker::Clear(option);
306 }
307 //_____________________________________________________________________________
308 Int_t StLaserAnalysisMaker::Make(){
309  StEvent* pEvent = dynamic_cast<StEvent*> (GetInputDS("StEvent"));
310  // Fill the event header.
311 
312  StEvtHddr *EvtHddr = (StEvtHddr*)GetDataSet("EvtHddr");
313  if (! EvtHddr) return kStWarn;
314  event->SetHeader(EvtHddr->GetEventNumber(), EvtHddr->GetRunNumber(), EvtHddr->GetDate(), EvtHddr->GetTime(),
315  gStTpcDb->Electronics()->tZero(), gStTpcDb->DriftVelocity(24,0), gStTpcDb->Electronics()->samplingFrequency(),
316  EvtHddr->GetInputTriggerMask());
317  event->SetDVWest(gStTpcDb->DriftVelocity(1,0));
318  event->SetDVEast(gStTpcDb->DriftVelocity(13,0));
319 #if 0
320  event->SetScaleY(gStTpcDb->ScaleY());
321 #endif
322  event->GetHeader()->SetDriftDistance(gStTpcDb->Dimensions()->gatingGridZ());
323  event->GetHeader()->SetInnerSectorzOffset(gStTpcDb->Dimensions()->zInnerOffset());
324  event->GetHeader()->SetOuterSectorzOffset(gStTpcDb->Dimensions()->zOuterOffset());
325  event->GetHeader()->SettriggerTimeOffset(gStTpcDb->triggerTimeOffset());
326  event->GetHeader()->SetBField(pEvent->runInfo()->magneticField());
327  StDetectorDbClock* dbclock = StDetectorDbClock::instance();
328  double freq = dbclock->getCurrentFrequency()/1000000.0;
329  event->GetHeader()->SetOnlClock(freq);
330  if (! pEvent) return kStWarn;
331  UInt_t nvtx = pEvent->numberOfPrimaryVertices();
332  for (UInt_t i = 0; i < nvtx; i++) {
333  event->AddVertex(pEvent->primaryVertex(i));
334  }
335  StSPtrVecTrackNode& trackNode = pEvent->trackNodes();
336  UInt_t nTracks = trackNode.size();
337  StTrackNode *node=0;
338  for (unsigned int i=0; i < nTracks; i++) {
339  node = trackNode[i];
340  if (!node) continue;
341  StGlobalTrack *gTrack = static_cast<StGlobalTrack *>(node->track(global));
342  if (! gTrack) continue;
343  if (!gTrack->dcaGeometry()) continue;
344  Int_t key = gTrack->key();
345  // if (gTrack->numberOfPossiblePoints(kTpcId) < 25) continue;
346  StPrimaryTrack *pTrack = static_cast<StPrimaryTrack*>(node->track(primary));
347  StThreeVectorD g3 = gTrack->outerGeometry()->momentum();
348  if (g3.mag() < 10) continue;
349  StThreeVectorD unit = g3.unit();
350  StPhysicalHelixD helixO = gTrack->outerGeometry()->helix();
351  // find a match using last hit
352  StPtrVecHit hvec = gTrack->detectorInfo()->hits();
353  Double_t rMax = 0;
354  Int_t jmax = -1;
355  // StThreeVectorD *pred = 0;
356  StTpcHit *tpcHit = 0;
357  Int_t sector = -1, s = -1;
358  Int_t bundle = -1;
359  Double_t dZ, dZmin = 9999;
360  Int_t b, m;
361  Double_t dPhi, dPhimin = 999;
362  Int_t mmax = -1;
363  StThreeVectorD Pred, Diff;
364  Track *t = 0;
365  for (unsigned int j=0; j<hvec.size(); j++) {// hit loop
366  if (hvec[j]->detector() != kTpcId) continue;
367  tpcHit = static_cast<StTpcHit *> (hvec[j]);
368 #ifndef __TRACKHITS__
369 #ifdef ADDPRIMTRACKHITS
370  if (pTrack)
371  event->AddHit(tpcHit, key);
372 #endif
373 #else
374  event->AddHit(tpcHit,key);
375 #endif
376  if (tpcHit->position().perp() > rMax) {
377  rMax = tpcHit->position().perp();
378  jmax = j;
379  }
380  }
381  LaserB *theLaser = 0;
382  if (jmax < 0 || rMax < 120) goto ADDTRACK;
383  tpcHit = static_cast<StTpcHit *> (hvec[jmax]);
384  // sector
385  sector = tpcHit->sector();
386  if (pTrack) goto ADDTRACK;
387  if (2*(sector/2) != sector) goto ADDTRACK;
388  s = sector/2 - 1;
389  if (s < 0 || s >= NS) goto ADDTRACK;
390  // bundle
391  dZmin = 9999;
392  for (b = 0; b < NB; b++) {
393  if (! Lasers[s][b][0]) continue;
394  dZ = TMath::Abs(tpcHit->position().z() - Lasers[s][b][0]->XyzG.z());
395  if (dZ < dZmin) {bundle = b; dZmin = dZ;}
396  }
397  if (bundle < 0 || dZmin > 15) goto ADDTRACK;
398  // mirror
399  // minimum distance in XY plane from laser spot
400  dPhimin = 999;
401  for (m = 0; m < NM; m++) {
402  if (! Lasers[s][bundle][m]) continue;
403  dPhi = TMath::Abs(Lasers[s][bundle][m]->PhiG - g3.phi());
404  if (dPhi < dPhimin) {
405  dPhimin = dPhi;
406  mmax = m;
407  }
408  }
409  if (mmax < 0 || dPhimin > 0.1) goto ADDTRACK;
410  theLaser = Lasers[s][bundle][mmax];
411  ADDTRACK:
412  t = event->AddTrack(sector,gTrack,theLaser,tpcHit->position().z());
413  if (theLaser) {
414  Int_t raft = theLaser->Raft;
415  Int_t bundle = theLaser->Bundle;
416  Int_t mirror = theLaser->Mirror;
417  if (raft > 0 && raft <= 14 &&
418  bundle > 0 && bundle <= 6 &&
419  mirror > 0 && mirror <= 7) {
420  t->SetPredictions(Raft2Tpc[raft-1], Bundles2Tpc[raft-1][bundle-1], Mirrors2Tpc[raft-1][bundle-1][mirror-1]);
421  if (theLaser->IsValid) event->AddTrackFit(t);
422  }
423  }
424  if (Debug()) {
425  t->Print();
426  }
427  }
428  static const Double_t sigma = 0.0385; // run 8090018 Full Field.
429  for (Int_t k = 0; k < 11; k++) {
430  FitDV *fit = (FitDV *) (*event->Fit())[k];
431  Bool_t ok = kTRUE;
432  Int_t N = fit->N;
433  if (N > 3) {
434  for (;;) {
435  TRVector Y;
436  TRMatrix A(0,2);
437  for (Int_t i = 0; i < N; i++) {
438  if (! fit->Flag[i]) {
439  Double_t dev = fit->Y[i]/sigma;
440  Y.AddRow(&dev);
441  Double_t a[2] = {1./sigma, fit->X[i]/sigma};
442  A.AddRow(a);
443  }
444  }
445  Int_t ndf = A.GetNrows() - 2;
446  if (ndf < 2) break;
447  TRSymMatrix S(A,TRArray::kATxA); if (Debug()) cout << "S: " << S << endl;
448  TRVector B(A,TRArray::kATxB,Y); if (Debug()) cout << "B: " << B << endl;
449  TRSymMatrix SInv(S,TRArray::kInvertedA);if (Debug()) cout << "SInv: " << SInv << endl;
450  if (! SInv.IsValid()) {break;}
451  TRVector X(SInv,TRArray::kSxA,B); if (Debug()) cout << "X: " << X << endl;
452  TRVector R(Y);
453  R -= TRVector(A,TRArray::kAxB,X); if (Debug()) cout << "R: " << R << endl;
454  Double_t chisq = R*R;
455  Double_t prob = TMath::Prob(chisq,ndf);
456  if (prob > 0.01) {
457  fit->offset = X[0];
458  fit->slope = X[1];
459  fit->chisq = chisq;
460  fit->dslope = SInv[2];
461  fit->doffset = SInv[0];
462  fit->Prob = prob;
463  fit->ndf = ndf;
464  if (Debug()) fit->Print();
465 #if 0
466  TClonesArray &tracks = *event->Tracks();
467  Int_t Ntrack = event->GetNtrack();
468  for (Int_t i = 0; i < Ntrack; i++) {
469  Track *t = (Track *) tracks[i];
470  if (t->Flag == 2 && t->Laser.Sector == 2*(k+1)) {
471  t->fit = *fit;
472  }
473  }
474 #endif
475  break;
476  } else {
477  Int_t j = -1;
478  Int_t imax = -1;
479  Double_t Rmax = 0;
480  for (Int_t i = 0; i < N; i++) {
481  if (! fit->Flag[i]) {
482  j++;
483  Double_t RR = TMath::Abs(R[j]);
484  if (RR > Rmax) {
485  imax = i;
486  Rmax = RR;
487  }
488  }
489  }
490  if (imax < 0) break;
491  fit->Flag[imax] = 1;
492  }
493  }
494  }
495  if (! ok) break;
496  }
497  m_laser->Fill(); //Fill the Tree
498  return kStOK;
499 }
500 //________________________________________________________________________________
501 Int_t StLaserAnalysisMaker::Finish() {
502  for (Int_t s = 0; s < NS; s++)
503  for (Int_t b = 0; b < NB; b++)
504  for (Int_t m = 0; m < NM; m++) SafeDelete(Lasers[s][b][m]);
505 #if 0
506  if (! GetTFile()) {
507  TSeqCollection *files = gROOT->GetListOfFiles();
508  if (files) {
509  TFile *f = 0;
510  TIter next(files);
511  while ((f = (TFile *) next())) {
512  TString name(gSystem->BaseName(f->GetName()));
513  if (name == "StLaserAnalysis.root") {
514  delete f;
515  break;
516  }
517  }
518  }
519  }
520 #endif
521  return StMaker::Finish();
522 }
St_starClockOnlC
Definition: St_starClockOnlC.h:7
Track::fit
bool fit()
Perform linear fits in zx- and zy-plane.
Definition: Track.cc:26
StMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
LASERINO_t
Definition: laserino.h:39
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
Track
C++ STL includes.
Definition: AgUStep.h:47
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5