StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StiCATpcTrackerInterface.cxx
1 #include "StiCATpcTrackerInterface.h"
2 #include "TPCCATracker/AliHLTTPCCAGBHit.h"
3 #include "TPCCATracker/AliHLTTPCCAGBTrack.h"
4 #include "TPCCATracker/AliHLTTPCCAParam.h"
5  // need for hits data
6 #include "StTpcHit.h"
7 #include "StTpcDb/StTpcDb.h"
8 #include "StDbUtilities/StTpcCoordinateTransform.hh"
9 #include "StDbUtilities/StTpcLocalSectorCoordinate.hh"
10  // for MCdata
11 #include "tables/St_g2t_track_Table.h"
12 #include "tables/St_g2t_tpc_hit_Table.h"
13 #include "TDatabasePDG.h"
14 #include "StBFChain.h"
15 #include "Sti/StiDetectorBuilder.h"
16 #include "Sti/StiDetectorGroups.h"
17 #include "Sti/StiGenericDetectorGroup.h"
18 #include "Sti/StiToolkit.h"
19  //to obtain error coefficients
20 #include "StDetectorDbMaker/StiTpcInnerHitErrorCalculator.h"
21 #include "StDetectorDbMaker/StiTpcOuterHitErrorCalculator.h"
22  //to get Magnetic Field
23 #include "StarMagField/StarMagField.h"
24 
25  // for sti perfo
26 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
27 #include "TPCCATrackerPerformance/AliHLTTPCCAStiPerformance.h"
28 #include "TPCCATrackerPerformance/AliHLTTPCCAMergerPerformance.h"
29 #include "StDetectorDbMaker/St_tpcPadConfigC.h"
30 #include "Sti/StiKalmanTrack.h"
31 #include "Sti/StiKalmanTrackNode.h"
32 #endif /* DO_TPCCATRACKER_EFF_PERFORMANCE */
33 #include <vector>
34 #include <algorithm>
35 using std::vector;
36 
37 #include <string>
38 using std::string;
39 
40 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
41 //#define STORE_STANDALONE_DATA // write data in files for Standalone
42 #endif
43 
44 StiCATpcTrackerInterface &StiCATpcTrackerInterface::Instance()
45 {
46  // reference to static object
47  static StiCATpcTrackerInterface g;
48  return g;
49 }
50 
51 StiCATpcTrackerInterface::StiCATpcTrackerInterface()
52 {
53 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
54  fOutFile = StMaker::GetChain()->GetTFile();
55  if(!fOutFile) cout << "W StiCATpcTrackerInterface: Warning - There isn't any tag file, so histograms won't be saved!" << endl;
56 
57  fPerformance = &(AliHLTTPCCAPerformance::Instance());
58  fPerformance->SetOutputFile(fOutFile);
59 #endif
60  //yf SetNewEvent();
61 } // StiCATpcTrackerInterface::StiCATpcTrackerInterface()
62 
63 StiCATpcTrackerInterface::~StiCATpcTrackerInterface( )
64 { // never called for static object
65 } // StiCATpcTrackerInterface::StiCATpcTrackerInterface()
66 
67 void StiCATpcTrackerInterface::SetNewEvent()
68 {
69  fHitsMap = 0;
70  fSeeds.clear();
71 
72  fSeedFinder = 0;
73  fStiTracks = 0;
74 
75  fIdTruth.clear(); // id of the Track, which has created CaHit
76 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
77  assert(fPerformance != 0);
78  fMCTracks.clear();
79  fMCPoints.clear();
80  fHitLabels.clear();
81 #endif
82  fCaParam.clear();// settings for all sectors to give CATracker
83  fCaHits.clear(); // hits to give CATracker
84  fSeedHits.clear(); // hits to make seeds
85 
86 //VP if (!fSeedFinder) fSeedFinder = new StiCATpcSeedFinder;
87 
88  if (fTracker) delete fTracker;
89  fTracker = new AliHLTTPCCAGBTracker;
90  if (fStiTracker) delete fStiTracker;
91  fStiTracker = new AliHLTTPCCAGBTracker;
92 
93 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
94  fPerformance->SetTracker(fTracker);
95 #endif
96 }
97 
98 
100 void StiCATpcTrackerInterface::Run()
101 {
102  assert(fHitsMap != 0);
103 
104  TStopwatch timer;
105  timer.Start();
106 
107  MakeSettings();
108  MakeHits();
109 
110 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
111  // put data in performance
112  FillPerformance(fCaHits,fIdTruth, fMCTracks,fMCPoints,fHitLabels);
113  fPerformance->SetMCTracks(fMCTracks);
114  fPerformance->SetMCPoints(fMCPoints);
115  fPerformance->SetHitLabels(fHitLabels);
116 #endif // DO_TPCCATRACKER_EFF_PERFORMANCE
117 
118 
119  // run tracker
120  fTracker->SetSettings(fCaParam);
121  fTracker->SetHits(fCaHits);
122 
123 #ifdef STORE_STANDALONE_DATA // write data in files for Standalone
124  static int iEvent = -1;
125  iEvent++;
126  TString name = "./data/";
127  if (iEvent == 0) fTracker->SaveSettingsInFile(string(name));
128  name += "event";
129  name += iEvent;
130  name += "_";
131  fTracker->SaveHitsInFile(string(name));
132 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
133  fPerformance->SaveDataInFiles(string(name));
134 #endif
135 // check
136  if(1){
137  if (fTracker) delete fTracker;
138  fTracker = new AliHLTTPCCAGBTracker;
139 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
140  fPerformance->SetTracker(fTracker);
141 #endif
142  TString name = "./data/";
143  fTracker->ReadSettingsFromFile(string(name));
144  name += "event";
145  name += iEvent;
146  name += "_";
147  fTracker->ReadHitsFromFile(string(name));
148 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
149  fPerformance->ReadDataFromFiles(string(name));
150 #endif
151  fTracker->SetSettings(fCaParam);
152  fTracker->SetHits(fCaHits);
153 
154  }
155 #endif // STORE_STANDALONE_DATA
156 
157 
158 
159  timer.Stop();
160  fPreparationTime_real = timer.RealTime();
161  fPreparationTime_cpu = timer.CpuTime();
162 
163  fTracker->FindTracks();
164 
165  // copy hits
166  timer.Start();
167 
168  MakeSeeds();
169 
170  timer.Stop();
171  fPreparationTime_real += timer.RealTime();
172  fPreparationTime_cpu += timer.CpuTime();
173 } // void StiCATpcTrackerInterface::Run()
174 
175 void StiCATpcTrackerInterface::RunPerformance()
176 {
177  cout << " ---- CA TPC Tracker ---- " << endl;
178 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
179  TStopwatch timer;
180  timer.Start();
181 
182  assert(fStiTracks != 0);
183  assert(fStiTracker != 0);
184 
185  // ----- Efficiency ---------
186 
187  // init all performances
188  fPerformance->InitSubPerformances();
189 
190  // -- correct init of stiPerformance --
191 
192  FillStiPerformance(); // get fStiCaHits fStiIdTruth fStiCaTracks
193 
194  // prepare tracker
195  fStiTracker->SetHits(fStiCaHits);
196  fStiTracker->SetSettings(fCaParam);
197 
198  AliHLTTPCCAStiPerformance* stiPerf = (AliHLTTPCCAStiPerformance*)fPerformance->GetSubPerformance("Sti Performance");
199  stiPerf->SetNewEvent(fStiTracker,&fStiHitLabels,&fStiMCTracks,&fStiMCPoints);
200  stiPerf->SetRecoTracks(fStiCaTracks);
201  fPerformance->ExecPerformance();
202 
203  timer.Stop();
204  fPreparationTime_real += timer.RealTime();
205  fPreparationTime_cpu += timer.CpuTime();
206 #endif // DO_TPCCATRACKER_EFF_PERFORMANCE
207 
208  // ----- Timing ---------
209  // TODO saparete in procedure in TPCCATracker/ ...
210 #if 1
211  const int printTime = 3; // 0 - show nothing; 1 - show avarege; 2 - show cur event; 3 - show both
212  const int fullTiming = 1; // show info about each stage of tracking; same as printTime
213  if ((printTime == 2) || (printTime == 3)){
214 
215  std::cout << "Sector reconstruction Time"
216  << " Real = " << std::setw( 10 ) << fTracker->SliceTrackerTime() * 1.e3 << " ms,"
217  << " CPU = " << std::setw( 10 ) << fTracker->SliceTrackerCpuTime() * 1.e3 << " ms,";
218  if (fTracker->SliceTrackerTime() > 0)
219  std::cout << " parallelization speedup: " << fTracker->SliceTrackerCpuTime() / fTracker->SliceTrackerTime();
220  std::cout << std::endl;
221  if ((fullTiming == 2) || (fullTiming == 3)) {
222  std::cout
223  << " | sum slice trackers: " << std::setw( 10 ) << fTracker->StatTime( 0 ) * 1000. << " ms\n"
224  << " | NeighboursFinder: " << std::setw( 10 ) << fTracker->StatTime( 1 ) * 1000. << " ms, " << std::setw( 12 ) << fTracker->StatTime( 5 ) << " cycles\n"
225  << " | StartHitsFinder: " << std::setw( 10 ) << fTracker->StatTime( 4 ) * 1000. << " ms\n"
226  << " | TrackletConstructor: " << std::setw( 10 ) << fTracker->StatTime( 2 ) * 1000. << " ms, " << std::setw( 12 ) << fTracker->StatTime( 7 ) << " cycles\n"
227  << " | TrackletSelector: " << std::setw( 10 ) << fTracker->StatTime( 3 ) * 1000. << " ms, " << std::setw( 12 ) << fTracker->StatTime( 8 ) << " cycles\n"
228  << " | WriteOutput: " << std::setw( 10 ) << fTracker->StatTime( 6 ) * 1000. << " ms\n"
229  << " | merge: " << std::setw( 10 ) << fTracker->StatTime( 9 ) * 1000. << " ms" << std::endl;
230  }
231  std::cout << "Total (sector+merge) reconstuction time"
232  << " Real = " << std::setw( 10 ) << (fTracker->SliceTrackerTime()+fTracker->StatTime( 9 )) * 1.e3 << " ms,"
233  << " CPU = " << std::setw( 10 ) << (fTracker->SliceTrackerCpuTime()+fTracker->StatTime( 10 )) * 1.e3 << " ms"
234  << std::endl;
235  std::cout << "Preparation time"
236  << " Real = " << std::setw( 10 ) << fPreparationTime_real * 1.e3 << " ms,"
237  << " CPU = " << std::setw( 10 ) << fPreparationTime_cpu * 1.e3 << " ms"
238  << std::endl;
239  };
240  if ((printTime == 1) || (printTime == 3)){
241  const int NTimers = 20;
242  static int statIEvent = 0;
243  static double statTime[NTimers];
244  static double statTime_SliceTrackerTime = 0;
245  static double statTime_SliceTrackerCpuTime = 0;
246  static double statPreparationTime_real = 0;
247  static double statPreparationTime_cpu = 0;
248 
249  if (!statIEvent){
250  for (int i = 0; i < NTimers; i++){
251  statTime[i] = 0;
252  }
253  }
254 
255  statIEvent++;
256  for (int i = 0; i < NTimers; i++){
257  statTime[i] += fTracker->StatTime( i );
258  }
259  statTime_SliceTrackerTime += fTracker->SliceTrackerTime();
260  statTime_SliceTrackerCpuTime += fTracker->SliceTrackerCpuTime();
261  statPreparationTime_real += fPreparationTime_real;
262  statPreparationTime_cpu += fPreparationTime_cpu;
263 
264  std::cout << "Average sector reconstruction Time"
265  << " Real = " << std::setw( 10 ) << 1./statIEvent*statTime_SliceTrackerTime * 1.e3 << " ms,"
266  << " CPU = " << std::setw( 10 ) << 1./statIEvent*statTime_SliceTrackerCpuTime * 1.e3 << " ms,";
267  if (statTime_SliceTrackerTime > 0)
268  std::cout << " parallelization speedup: " << statTime_SliceTrackerCpuTime / statTime_SliceTrackerTime;
269  std::cout << std::endl;
270  if ((fullTiming == 1) || (fullTiming == 3)) {
271  std::cout
272  << " | sum slice trackers: " << std::setw( 10 ) << 1./statIEvent*statTime[ 0 ] * 1000. << " ms\n"
273  << " | NeighboursFinder: " << std::setw( 10 ) << 1./statIEvent*statTime[ 1 ] * 1000. << " ms, " << std::setw( 12 ) << 1./statIEvent*statTime[ 5 ] << " cycles\n"
274  << " | StartHitsFinder: " << std::setw( 10 ) << 1./statIEvent*statTime[ 4 ] * 1000. << " ms\n"
275  << " | TrackletConstructor: " << std::setw( 10 ) << 1./statIEvent*statTime[ 2 ] * 1000. << " ms, " << std::setw( 12 ) << 1./statIEvent*statTime[ 7 ] << " cycles\n"
276  << " | TrackletSelector: " << std::setw( 10 ) << 1./statIEvent*statTime[ 3 ] * 1000. << " ms, " << std::setw( 12 ) << 1./statIEvent*statTime[ 8 ] << " cycles\n"
277  << " | WriteOutput: " << std::setw( 10 ) << 1./statIEvent*statTime[ 6 ] * 1000. << " ms\n"
278  << " | merge: " << std::setw( 10 ) << 1./statIEvent*statTime[ 9 ] * 1000. << " ms" << std::endl;
279  }
280  std::cout << "Total (sector+merge) avarage reconstuction time"
281  << " Real = " << std::setw( 10 ) << (statTime_SliceTrackerTime+statTime[ 9 ])/statIEvent * 1.e3 << " ms,"
282  << " CPU = " << std::setw( 10 ) << (statTime_SliceTrackerCpuTime+statTime[ 10 ])/statIEvent * 1.e3 << " ms"
283  << std::endl;
284  std::cout << "Avarage preparation time"
285  << " Real = " << std::setw( 10 ) << statPreparationTime_real/statIEvent * 1.e3 << " ms,"
286  << " CPU = " << std::setw( 10 ) << statPreparationTime_cpu/statIEvent * 1.e3 << " ms."
287  << std::endl;
288  }
289 #endif // 0 timing
290 
291 #if 0//def DO_TPCCATRACKER_EFF_PERFORMANCE outdated
292  ((AliHLTTPCCAMergerPerformance*)(AliHLTTPCCAPerformance::Instance().GetSubPerformance("Merger")))->FillTree();
293 #endif //DO_TPCCATRACKER_EFF_PERFORMANCE
294 } // void StiCATpcTrackerInterface::Run()
295 
296 
297 void StiCATpcTrackerInterface::MakeSettings()
298 {
299 
300  const int NSlices = 24; //TODO initialize from StRoot
301  for ( int iSlice = 0; iSlice < NSlices; iSlice++ ) {
302  AliHLTTPCCAParam SlicePar;
303  // memset(&SlicePar, 0, sizeof(AliHLTTPCCAParam));
304 
305  Int_t sector = iSlice+1;
306  // Int_t sector = iSlice;
307  const int NoOfInnerRows = St_tpcPadConfigC::instance()->innerPadRows(sector);
308  const int NRows = St_tpcPadConfigC::instance()->padRows(sector);
309  SlicePar.SetISlice( iSlice );
310  SlicePar.SetNRows ( NRows );
311  SlicePar.SetNInnerRows ( NoOfInnerRows );
312  Double_t beta = 0;
313  if (sector > 12) beta = (24-sector)*2.*TMath::Pi()/12.;
314  else beta = sector *2.*TMath::Pi()/12.;
315  SlicePar.SetAlpha ( beta );
316  SlicePar.SetDAlpha ( 30*TMath::DegToRad() ); //TODO initialize from StRoot
317  SlicePar.SetCosAlpha ( TMath::Cos(SlicePar.Alpha()) );
318  SlicePar.SetSinAlpha ( TMath::Sin(SlicePar.Alpha()) );
319  SlicePar.SetAngleMin ( SlicePar.Alpha() - 0.5*SlicePar.DAlpha() );
320  SlicePar.SetAngleMax ( SlicePar.Alpha() + 0.5*SlicePar.DAlpha() );
321  SlicePar.SetRMin ( 51. ); //TODO initialize from StRoot
322  SlicePar.SetRMax ( 194. ); //TODO initialize from StRoot
323  SlicePar.SetErrX ( 0. ); //TODO initialize from StRoot
324  SlicePar.SetErrY ( 0.12 ); // 0.06 for Inner //TODO initialize from StRoot
325  SlicePar.SetErrZ ( 0.16 ); // 0.12 for Inner NodePar->fitPars() //TODO initialize from StRoot
326  // SlicePar.SetPadPitch ( 0.675 );// 0.335 -"-
327  if (! StiKalmanTrackNode::IsLaser()) {
328  float x[3]={0,0,0},b[3];
329  StarMagField::Instance()->BField(x,b);
330  SlicePar.SetBz ( - b[2] ); // change sign because change z
331  } else {
332  static const double HZERO=2e-6;
333  SlicePar.SetBz (HZERO);
334  }
335  if (sector <= 12) {
336  SlicePar.SetZMin ( 0. ); //TODO initialize from StRoot
337  SlicePar.SetZMax ( 210. ); //TODO initialize from StRoot
338  } else {
339  SlicePar.SetZMin (-210. ); //TODO initialize from StRoot
340  SlicePar.SetZMax ( 0. ); //TODO initialize from StRoot
341  }
342  for( int iR = 0; iR < NRows; iR++){
343  SlicePar.SetRowX(iR, St_tpcPadConfigC::instance()->radialDistanceAtRow(sector,iR+1));
344  }
345 
346  Double_t *coeffInner = StiTpcInnerHitErrorCalculator::instance()->coeff();
347  for(int iCoef=0; iCoef<6; iCoef++)
348  {
349  SlicePar.SetParamS0Par(0, 0, iCoef, (float)coeffInner[iCoef] );
350  }
351 
352  SlicePar.SetParamS0Par(0, 0, 6, 0.0f );
353 
354  Double_t *coeffOuter =StiTpcOuterHitErrorCalculator::instance()->coeff();
355  for(int iCoef=0; iCoef<6; iCoef++)
356  {
357  SlicePar.SetParamS0Par(0, 1, iCoef, (float)coeffOuter[iCoef] );
358  }
359 
360  SlicePar.SetParamS0Par(0, 1, 6, 0.0f );
361  SlicePar.SetParamS0Par(0, 2, 0, 0.0f );
362  SlicePar.SetParamS0Par(0, 2, 1, 0.0f );
363  SlicePar.SetParamS0Par(0, 2, 2, 0.0f );
364  SlicePar.SetParamS0Par(0, 2, 3, 0.0f );
365  SlicePar.SetParamS0Par(0, 2, 4, 0.0f );
366  SlicePar.SetParamS0Par(0, 2, 5, 0.0f );
367  SlicePar.SetParamS0Par(0, 2, 6, 0.0f );
368  SlicePar.SetParamS0Par(1, 0, 0, 0.0f );
369  SlicePar.SetParamS0Par(1, 0, 1, 0.0f );
370  SlicePar.SetParamS0Par(1, 0, 2, 0.0f );
371  SlicePar.SetParamS0Par(1, 0, 3, 0.0f );
372  SlicePar.SetParamS0Par(1, 0, 4, 0.0f );
373  SlicePar.SetParamS0Par(1, 0, 5, 0.0f );
374  SlicePar.SetParamS0Par(1, 0, 6, 0.0f );
375  SlicePar.SetParamS0Par(1, 1, 0, 0.0f );
376  SlicePar.SetParamS0Par(1, 1, 1, 0.0f );
377  SlicePar.SetParamS0Par(1, 1, 2, 0.0f );
378  SlicePar.SetParamS0Par(1, 1, 3, 0.0f );
379  SlicePar.SetParamS0Par(1, 1, 4, 0.0f );
380  SlicePar.SetParamS0Par(1, 1, 5, 0.0f );
381  SlicePar.SetParamS0Par(1, 1, 6, 0.0f );
382  SlicePar.SetParamS0Par(1, 2, 0, 0.0f );
383  SlicePar.SetParamS0Par(1, 2, 1, 0.0f );
384  SlicePar.SetParamS0Par(1, 2, 2, 0.0f );
385  SlicePar.SetParamS0Par(1, 2, 3, 0.0f );
386  SlicePar.SetParamS0Par(1, 2, 4, 0.0f );
387  SlicePar.SetParamS0Par(1, 2, 5, 0.0f );
388  SlicePar.SetParamS0Par(1, 2, 6, 0.0f );
389 
390  fCaParam.push_back(SlicePar);
391  } // for iSlice
392 } // void StiCATpcTrackerInterface::MakeSettings()
393 
394 
395 void StiCATpcTrackerInterface::MakeHits()
396 {
397  StTpcCoordinateTransform tran(gStTpcDb);
398  StTpcLocalSectorCoordinate loc;
399  for (HitMapToVectorAndEndType::iterator it= fHitsMap->begin(); it!= fHitsMap->end(); ++it){
400  vector<StiHit*>& tempvec = (*it).second.hits();
401  vector<StiHit*>::iterator start = tempvec.begin();
402  vector<StiHit*>::iterator stop = tempvec.end();
403  for (vector<StiHit*>::iterator cit = start; cit != stop; cit++) {
404 
405  // get local coordinates. take into account distortion
406  StiHit *hit = *cit;
407  if ( !(hit->detector()->isActive()) ) continue;
408  if (! hit->stHit()) continue;
409  if ( hit->timesUsed()) continue;//VP
410 
411  const StTpcHit *tpcHit = dynamic_cast<const StTpcHit*>(hit->stHit());
412  if ( ! tpcHit) continue;
413  StGlobalCoordinate glob(tpcHit->position());
414  tran(glob,loc,tpcHit->sector(),tpcHit->padrow());
415 
416  // obtain seed Hit
417  SeedHit_t hitc;
418 #if 0
419  hitc.mMinPad = tpcHit->minPad();
420  hitc.mMaxPad = tpcHit->maxPad();
421  hitc.mMinTmbk = tpcHit->minTmbk();
422  hitc.mMaxTmbk = tpcHit->maxTmbk();
423 #endif
424  hitc.padrow = tpcHit->padrow()-1;
425  hitc.x = loc.position().x();
426  hitc.y = loc.position().y();
427  hitc.z = loc.position().z();
428  hitc.status=0;
429  hitc.taken=0;
430  hitc.track_key=tpcHit->idTruth();
431  hitc.hit = hit;
432  fSeedHits.push_back(hitc);
433 
434  // convert to CA Hit
435  AliHLTTPCCAGBHit caHit;
436  caHit.SetIRow( hitc.padrow );
437 // caHit.SetX( hit->x() );
438  caHit.SetX( hit->position() ); // take position of the row
439  caHit.SetY( - hit->y() );
440  caHit.SetZ( - hit->z() );
441  // caHit.SetErrX( );
442  caHit.SetErrY( 0.12 );// TODO: read parameters from somewhere
443  caHit.SetErrZ( 0.16 );
444  caHit.SetISlice( tpcHit->sector() - 1 );
445  caHit.SetIRow( hitc.padrow );
446  caHit.SetID( fCaHits.size() );
447  fIdTruth.push_back( hitc.track_key );
448 
449  fCaHits.push_back(caHit);
450  }
451  }
452 
453 } // void StiCATpcTrackerInterface::MakeHits()
454 
455 void StiCATpcTrackerInterface::ConvertPars(const AliHLTTPCCATrackParam& caPar, double _alpha, StiNodePars& nodePars, StiNodeErrs& nodeErrs)
456 {
457  // set jacobian integral coef
458  double JI[5];
459  JI[0] = -1.; // y
460  JI[1] = -1.; // z
461  JI[2] = -1.; // eta
462  JI[3] = 1.; // ptin
463  JI[4] = -1.; // tanl
464  // get parameters
465  nodePars.x() = caPar.GetX();
466  nodePars.y() = JI[0] * caPar.GetY();
467  nodePars.z() = JI[1] * caPar.GetZ();
468  nodePars.eta() = JI[2] * asin(caPar.GetSinPhi()); // (signed curvature)*(local Xc of helix axis - X current point on track)
469  nodePars.ptin() = JI[3] * caPar.GetQPt(); // signed invert pt [sign = sign(-qB)]
470  nodePars.tanl() = JI[4] * caPar.GetDzDs(); // tangent of the track momentum dip angle
471 
472  // get h & signB
473  static const double EC = 2.99792458e-4, ZEROHZ = 2e-6;
474 // #define USE_CA_FIELD
475 #ifndef USE_CA_FIELD // use field in the point. Need this because of check in StiTrackNodeHelper::set()
476  const double ca = cos(_alpha), sa = sin(_alpha); // code has been taken from StiKalmanTrackNode::getHz()
477  double globalXYZ[3] = {
478  ca * nodePars.x() - sa * nodePars.y(),
479  sa * nodePars.x() + ca * nodePars.y(),
480  nodePars.z()
481  };
482 
483  double h2=ZEROHZ;
484  if (! StiKalmanTrackNode::IsLaser()) {
485  double b[3];
486  StarMagField::Instance()->BField(globalXYZ,b);
487  h2 = b[2];
488  }
489 
490 #else // these parameters have been obtained with that MF, so let's use it.
491  double h2 = - fTracker->Slice(0).Param().Bz(); // change sign because change z
492 #endif // 1
493  h2 *= EC;
494  // get parameters. continue
495  nodePars.hz() = h2; // Z component magnetic field in units Pt(Gev) = Hz * RCurv(cm)
496  nodePars.ready(); // set cosCA, sinCA & curv
497 //std::cout << nodePars._ptin << std::endl;
498  // set jacobian integral coef
499  double J[5];
500  J[0] = JI[0]; // y
501  J[1] = JI[1]; // z
502  J[2] = JI[2]/cos(nodePars.eta()); // eta
503  J[3] = JI[3]; // ptin
504  J[4] = JI[4]; // tanl
505 
506  // get cov matrises
507  const float *caCov = caPar.GetCov();
508 // double nodeCov[15];
509 // for (int i1 = 0, i = 0; i1 < 5; i1++){
510 // for (int i2 = 0; i2 <= i1; i2++, i++){
511 // nodeCov[i] = J[i1]*J[i2]*caCov[i];
512 // }
513 // }
514  // if ( (caCov[0] <= 0) || (caCov[2] <= 0) || (caCov[5] <= 0) || (caCov[9] <= 0) || (caCov[14] <= 0))
515  // cout << "Warrning: Bad CA Cov Matrix." << endl;
516  // if ( (nodeCov[0] <= 0) || (nodeCov[2] <= 0) || (nodeCov[5] <= 0) || (nodeCov[9] <= 0) || (nodeCov[14] <= 0))
517  // cout << "Warrning: Bad Node Cov Matrix." << endl;
518 
519  double *A = nodeErrs.G();
520 /* for (int i1 = 0, i = 0; i1 < 5; i1++){
521  for (int i2 = 0; i2 <= i1; i2++, i++){
522  A[i+i1+2] = caCov[i];
523  }
524  }*/
525 
526  nodeErrs._cYY = caCov[ 0];
527  nodeErrs._cZY = caCov[ 1]*J[0]*J[1];
528  nodeErrs._cZZ = caCov[ 2]*J[0]*J[1];
529  nodeErrs._cEY = caCov[ 3]*J[0]*J[2];
530  nodeErrs._cEZ = caCov[ 4]*J[1]*J[2];
531  nodeErrs._cEE = caCov[ 5]*J[2]*J[2];
532  nodeErrs._cTY = caCov[ 6]*J[0]*J[4];
533  nodeErrs._cTZ = caCov[ 7]*J[1]*J[4];
534  nodeErrs._cTE = caCov[ 8]*J[2]*J[4];
535  nodeErrs._cTT = caCov[ 9]*J[4]*J[4];
536  nodeErrs._cPY = caCov[10]*J[0]*J[3];
537  nodeErrs._cPZ = caCov[11]*J[1]*J[3];
538  nodeErrs._cPE = caCov[12]*J[2]*J[3];
539  nodeErrs._cTP = caCov[13]*J[4]*J[3];
540  nodeErrs._cPP = caCov[14]*J[3]*J[3];
541 #if 1
542  A[0] = 1; // don't use parameter X
543  A[1] = 0;
544  A[3] = 0;
545  A[6] = 0;
546  A[10] = 0;
547  A[15] = 0;
548 #endif
549 }
550 
551 void StiCATpcTrackerInterface::MakeSeeds()
552 {
553  const int NRecoTracks = fTracker->NTracks();
554  for ( int iTr = 0; iTr < NRecoTracks; iTr++ ) {
555  // get seed
556  const AliHLTTPCCAGBTrack tr = fTracker->Track( iTr );
557  Seed_t seed;
558 
559  const int NHits = tr.NHits();
560 //VP float last_x = 1e10; // for check
561  for ( int iHit = NHits-1; iHit >= 0; iHit-- ){
562  const int index = fTracker->TrackHit( tr.FirstHitRef() + iHit );
563  const int hId = fTracker->Hit( index ).ID();
564 // if ( last_x == fSeedHits[hId].hit->position() ) continue; // track can have 2 hits on 1 row because of track segments merger.
565  seed.vhit.push_back(&(fSeedHits[hId]));
566 // assert( last_x >= fSeedHits[hId].hit->position() ); // should be back order - from outer to inner.
567 //VP last_x = fSeedHits[hId].hit->position();
568  }
569 
570  seed.total_hits = seed.vhit.size();
571 
572  ConvertPars( tr.OuterParam(), tr.Alpha(), seed.firstNodePars, seed.firstNodeErrs );
573  ConvertPars( tr.InnerParam(), tr.Alpha(), seed.lastNodePars, seed.lastNodeErrs );
574 
575  fSeeds.push_back(seed);
576  }
577 } // void StiCATpcTrackerInterface::MakeSeeds()
578 
579 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
580 bool myfunction (AliHLTTPCCALocalMCPoint i,AliHLTTPCCALocalMCPoint j) {
581  return (i.TrackI() < j.TrackI()) || ( i.TrackI()==j.TrackI() && i.IRow() < j.IRow() );
582 }
584 void StiCATpcTrackerInterface::FillPerformance(const vector<AliHLTTPCCAGBHit>& caHits, const vector<int>& idsTruth, vector<AliHLTTPCCAMCTrack>& mcTracks, vector<AliHLTTPCCALocalMCPoint>& mcPoints, vector<AliHLTTPCCAHitLabel>& hitLabels)
585 {
586  mcTracks.clear();
587  mcPoints.clear();
588  hitLabels.clear();
589  hitLabels.resize(caHits.size());
590 
591  StBFChain *chain = (StBFChain *) StMaker::GetChain();
592  St_g2t_track *Track = (St_g2t_track *) chain->FindObject("g2t_track");
593  St_g2t_tpc_hit *mcTpcHits = (St_g2t_tpc_hit *) chain->FindObject("g2t_tpc_hit");
594  if (Track && mcTpcHits) {
595  static StTpcLocalSectorCoordinate localSect;
596  static StTpcLocalSectorDirection dirSect;
597  StTpcCoordinateTransform transform(gStTpcDb);
598  int nMCTracks = Track->GetNRows();
599  const g2t_track_st *track = Track->GetTable();
600 
601  vector<int> truthToIndex; // need for finding MCPoints to MCTracks correspondence
602  truthToIndex.resize(nMCTracks+1,-1);
603 
604  const g2t_tpc_hit_st *mcTpcHit = mcTpcHits->GetTable();
605  const Int_t nMcTpcHits = mcTpcHits->GetNRows();
606 
607  AliHLTTPCCAMCTrack mcTrack;
608  for (Int_t l = 0; l < nMCTracks; l++, track++) {
609  if (! track->n_tpc_hit) continue;
610  if (track->ptot < 1e-3) continue;
611  memset(&mcTrack, 0, sizeof(AliHLTTPCCAMCTrack));
612  mcTrack.SetPDG( TDatabasePDG::Instance()->ConvertGeant3ToPdg(track->ge_pid) );
613  TParticlePDG *part = TDatabasePDG::Instance()->GetParticle( mcTrack.PDG() );
614  Float_t q = part->Charge();
615  mcTrack.SetP ( track->ptot );
616  mcTrack.SetPt( track->pt );
617  // mcTrack.SetNHits( track->n_tpc_hit ); // it is not really what we need. See below
618  for (Int_t i = 0; i < 3; i++) { // TODO
619  mcTrack.SetPar( i, ((i>=1)?-1:1) * 0 );
620  }
621  for (Int_t i = 0; i < 3; i++) {
622  mcTrack.SetPar( 3+i, ((i>=1)?-1:1) * track->p[i]/mcTrack.P() );
623  }
624  mcTrack.SetPar(6, q/mcTrack.P()/3 ); // q=3q
625  mcTrack.SetNTurns( 1 ); // TODO: read from somewhere
626  const int idTruth = l+1;
627 
628  truthToIndex[idTruth] = mcTracks.size();
629  mcTracks.push_back(mcTrack);
630  } // for iMCTrack
631 
632  // get hits
633  const int nHits = caHits.size();
634  for (int iHit = 0; iHit < nHits; iHit++){
635  const int iTrack = truthToIndex[idsTruth[iHit]];
636  if (iTrack == -1) continue;
637  AliHLTTPCCAHitLabel &l = hitLabels[iHit];
638  l.fLab[0] = iTrack;
639  l.fLab[1] = -1;
640  l.fLab[2] = -1;
641  mcTracks[iTrack].SetNHits( 1 + mcTracks[iTrack].NHits() );
642  } // j
643 
644  // get MCPoints
645  for (Int_t iP = 0; iP < nMcTpcHits; iP++, mcTpcHit++) {
646  AliHLTTPCCALocalMCPoint mcPoint;
647  Int_t is = mcTpcHit->volume_id/100000;
648  if (is) continue;
649  const int iTrack = truthToIndex[mcTpcHit->track_p];
650  if (iTrack == -1) continue;
651  AliHLTTPCCAMCTrack& mcTrack = mcTracks[iTrack];
652 
653  mcPoint.SetIRow( mcTpcHit->volume_id%100 - 1 );
654  mcPoint.SetISlice( (mcTpcHit->volume_id/100)%100 - 1 );
655 
656  Int_t sector = (mcTpcHit->volume_id%100000)/100;
657  Int_t row = mcTpcHit->volume_id%100;
658 
659  // get detector
660  StiDetectorBuilder* detectorBuilder = 0;
661  StiDetectorGroups *groups=StiToolkit::instance()->getDetectorGroups();
662  vector<StiGenericDetectorGroup *>::iterator it = groups->begin();
663  for (; it != groups->end(); ++it) {
664  StiGenericDetectorGroup *group = *it;
665  detectorBuilder = group->getDetectorBuilder();
666  TString builderName = (const char*)(detectorBuilder->getName().c_str());
667  if (builderName == "TpcBuilder")
668  break;
669  }
670  int sectorTpm = sector - 1;
671  if (sectorTpm >= 12) sectorTpm = 11 - (sectorTpm - 11)%12;
672  StiDetector* detector = detectorBuilder->getDetector(row-1,sectorTpm);
673 
674  static StiHit tmpStiHit; // used for convert only
675 
676  // Get XYZ
677  StTpcHit* mpTmp = new StTpcHit;
678  double energyTmp = 0.;
679  tmpStiHit.setGlobal(detector, mpTmp, mcTpcHit->x[0], mcTpcHit->x[1], mcTpcHit->x[2], energyTmp);
680  mcPoint.SetX( tmpStiHit.x() );
681  mcPoint.SetY( - tmpStiHit.y() );
682  mcPoint.SetZ( - tmpStiHit.z() );
683 
684  // Get P
685  StGlobalDirection gDir(mcTpcHit->p[0],mcTpcHit->p[1],mcTpcHit->p[2]);
686  transform(gDir,dirSect,sector,row);
687  double px = mcTpcHit->p[0], py = mcTpcHit->p[1], pz_new = mcTpcHit->p[2], px_new, py_new;
688  double angle = detector->getPlacement()->getNormalRefAngle();
689  px_new = cos(angle)*px + sin(angle)*py;
690  py_new = -sin(angle)*px + cos(angle)*py;
691 
692  mcPoint.SetPx( px_new );
693  mcPoint.SetPy( - py_new );
694  mcPoint.SetPz( - pz_new );
695  const int qpSign = fabs(mcTrack.Par(6))/mcTrack.Par(6);
696  mcPoint.SetQP( double( qpSign ) / sqrt(px_new*px_new + py_new*py_new + pz_new*pz_new) );
697  mcPoint.SetTrackI( iTrack );
698  mcPoint.SetTrackID( mcTpcHit->track_p );
699 
700  mcPoints.push_back(mcPoint);
701  } // for mcPoints
702 
703  std::sort( mcPoints.begin(), mcPoints.end(), myfunction );
704  mcTracks[0].SetFirstMCPointID(0);
705  int iPLast = 0;
706  int iTr = 0;
707  for( UInt_t iP = 0; iP < mcPoints.size(); ++iP ) {
708  if ( mcPoints[iP].TrackI() != iTr ) {
709  mcTracks[iTr].SetNMCPoints(iP - iPLast);
710  iPLast = iP;
711  iTr = mcPoints[iP].TrackI();
712  mcTracks[iTr].SetFirstMCPointID(iP);
713  }
714  }
715  mcTracks[iTr].SetNMCPoints(mcPoints.size() - iPLast);
716  }
717 } // void StiCATpcTrackerInterface::FillPerformance()
718 
719 void StiCATpcTrackerInterface::FillStiPerformance()
720 {
721  fStiCaHits.clear();
722  fStiIdTruth.clear();
723  fStiCaTracks.clear();
724 
725  AliHLTTPCCAStiPerformance* stiPerf = (AliHLTTPCCAStiPerformance*)fPerformance->GetSubPerformance("Sti Performance");
726  int NGBHits = 0;
727 
728 
729  StTpcCoordinateTransform tran(gStTpcDb);
730  StTpcLocalSectorCoordinate loc;
731 
732  for(int iTr=0; iTr<fStiTracks->getTrackCount(0); iTr++)
733  {
734  auto * track = (StiKalmanTrack*) fStiTracks->at(iTr);
735  vector<StiHit*> hits_v = track-> getHits();
736 
737  AliHLTTPCCAGBTrack GBTrack;
738  int nTrackHits = 0;
739 
740  // get local coordinates. take into account distortion
741 
742  for(unsigned iH=0; iH<hits_v.size(); iH++)
743  {
744  StiHit *hit = hits_v[iH];
745  if (! hit->stHit()) continue;
746  if ( hit->timesUsed()) continue;
747  const StTpcHit *tpcHit = dynamic_cast<const StTpcHit*>(hit->stHit());
748  if ( ! tpcHit) continue;
749  StGlobalCoordinate glob(tpcHit->position());
750  tran(glob,loc,tpcHit->sector(),tpcHit->padrow());
751 
752  // convert to CA Hit
753  AliHLTTPCCAGBHit caHit;
754  caHit.SetIRow( tpcHit->padrow()-1 );
755  caHit.SetX( hit->position() ); // take position of the row
756  caHit.SetY( - hit->y() );
757  caHit.SetZ( - hit->z() );
758  caHit.SetErrY( 0.12 );// TODO: read parameters from somewhere
759  caHit.SetErrZ( 0.16 );
760  caHit.SetISlice( tpcHit->sector() - 1 );
761  caHit.SetID( fStiCaHits.size() );
762  fStiIdTruth.push_back( tpcHit->idTruth() );
763 
764  int iHit = NGBHits+iH;
765  stiPerf->AddHit(iHit);
766  nTrackHits++;
767 
768  fStiCaHits.push_back(caHit);
769  }
770 
771  StiKalmanTrackNode *NodePar = track->getInnerMostDetHitNode(kTpcId);
772  if (!NodePar) continue;
773 
774  double JI[5];
775  JI[0] = -1.; // y
776  JI[1] = -1.; // z
777  JI[2] = -1.; // eta
778  JI[3] = 1.; // ptin
779  JI[4] = -1.; // tanl
780 
781  double J[5];
782  J[0] = JI[0]; // y
783  J[1] = JI[1]; // z
784  J[2] = JI[2]*cos(NodePar->fitPars().eta()); // eta
785  J[3] = JI[3]; // ptin
786  J[4] = JI[4]; // tanl
787 
788  double Cos = cos(NodePar->fitPars().eta());
789  float sgnCos = fabs(Cos)/Cos;
790 
791  AliHLTTPCCATrackParam GBTrackParam;
792  GBTrackParam.SetX ( NodePar->fitPars().x() );
793  GBTrackParam.SetY ( JI[0] * NodePar->fitPars().y() );
794  GBTrackParam.SetZ ( JI[1] * NodePar->fitPars().z() );
795  GBTrackParam.SetSinPhi( JI[2] * sin(NodePar->fitPars().eta()));
796  GBTrackParam.SetQPt ( JI[3] * NodePar->fitPars().ptin() );
797  GBTrackParam.SetDzDs ( JI[4] * NodePar->fitPars().tanl() );
798 
799  GBTrackParam.SetCov( 0, NodePar->fitErrs()._cYY);
800  GBTrackParam.SetCov( 1, NodePar->fitErrs()._cZY*J[0]*J[1]);
801  GBTrackParam.SetCov( 2, NodePar->fitErrs()._cZZ*J[0]*J[1]);
802  GBTrackParam.SetCov( 3, NodePar->fitErrs()._cEY*J[0]*J[2]);
803  GBTrackParam.SetCov( 4, NodePar->fitErrs()._cEZ*J[1]*J[2]);
804  GBTrackParam.SetCov( 5, NodePar->fitErrs()._cEE*J[2]*J[2]);
805  GBTrackParam.SetCov( 6, NodePar->fitErrs()._cTY*J[0]*J[4]);
806  GBTrackParam.SetCov( 7, NodePar->fitErrs()._cTZ*J[1]*J[4]);
807  GBTrackParam.SetCov( 8, NodePar->fitErrs()._cTE*J[2]*J[4]);
808  GBTrackParam.SetCov( 9, NodePar->fitErrs()._cTT*J[4]*J[4]);
809  GBTrackParam.SetCov(10, NodePar->fitErrs()._cPY*J[0]*J[3]);
810  GBTrackParam.SetCov(11, NodePar->fitErrs()._cPZ*J[1]*J[3]);
811  GBTrackParam.SetCov(12, NodePar->fitErrs()._cPE*J[2]*J[3]);
812  GBTrackParam.SetCov(13, NodePar->fitErrs()._cTP*J[4]*J[3]);
813  GBTrackParam.SetCov(14, NodePar->fitErrs()._cPP*J[3]*J[3]);
814 
815  GBTrackParam.SetSignCosPhi(sgnCos);
816 
817  GBTrack.SetInnerParam(GBTrackParam);
818 
819  GBTrack.SetNHits( nTrackHits );
820  GBTrack.SetFirstHitRef( NGBHits );
821 
822  NGBHits += nTrackHits;
823 
824  fStiCaTracks.push_back(GBTrack);
825  }
826 
827  // prepare mc
828  vector<AliHLTTPCCAMCTrack> stiMCTracks;
829  vector<AliHLTTPCCALocalMCPoint> stiMCPoints;
830  vector<AliHLTTPCCAHitLabel> stiHitLabels;
831  FillPerformance(fStiCaHits,fStiIdTruth, stiMCTracks,stiMCPoints,stiHitLabels);
832 
833  fStiMCTracks.Resize(stiMCTracks.size());
834  for (unsigned i = 0; i < stiMCTracks.size(); i++) {
835  fStiMCTracks[i] = stiMCTracks[i];
836  }
837  fStiMCPoints.Resize(stiMCPoints.size());
838  for (unsigned i = 0; i < stiMCPoints.size(); i++) {
839  fStiMCPoints[i] = stiMCPoints[i];
840  }
841  fStiHitLabels.Resize(stiHitLabels.size());
842  for (unsigned i = 0; i < stiHitLabels.size(); i++) {
843  fStiHitLabels[i] = stiHitLabels[i];
844  }
845 
846 
847 } // void StiCATpcTrackerInterface::FillStiPerformance()
848 #endif // DO_TPCCATRACKER_EFF_PERFORMANCE
Seed_t
Definition: StiCATpcSeedFinder.h:17
StiKalmanTrack
Definition of Kalman Track.
Definition: StiKalmanTrack.h:82
StiHit::setGlobal
void setGlobal(const StiDetector *detector, const StMeasuredPoint *stHit, Float_t x, Float_t y, Float_t z, Float_t energy)
Definition: StiHit.cxx:133
SeedHit_t
Definition: StiCATpcSeedFinder.h:9
StiHit
Definition: StiHit.h:51
StiHit::x
const Float_t & x() const
Return the local x, y, z values.
Definition: StiHit.h:67
StiHit::detector
const StiDetector * detector() const
Definition: StiHit.h:96
StiCATpcTrackerInterface::StiCATpcTrackerInterface
StiCATpcTrackerInterface()
constructor - destructor
Definition: StiCATpcTrackerInterface.cxx:51
AliHLTTPCCAGBTracker::SetHits
void SetHits(const std::vector< AliHLTTPCCAGBHit > &hits)
Try to group close hits in row formed by one track. After sort hits.
Definition: AliHLTTPCCAGBTracker.cxx:957
StiKalmanTrack.h
Definition of Kalman Track.
StiCATpcTrackerInterface::Run
void Run()
Copy data to CATracker. Run CATracker. Copy tracks in fSeeds.
Definition: StiCATpcTrackerInterface.cxx:100
StiHit::timesUsed
UInt_t timesUsed() const
Return the number of times this hit was assigned to a track.
Definition: StiHit.h:108
StiToolkit.h
Abstract interface for a STI toolkit.
StiTrackContainer::getTrackCount
int getTrackCount(Filter< StiTrack > *filter) const
Definition: StiTrackContainer.cxx:48
StiGenericDetectorGroup::getDetectorBuilder
virtual StiDetectorBuilder * getDetectorBuilder()
Get a detector builder appropriate for this detector group.
Definition: StiGenericDetectorGroup.cxx:33
StiHit::position
Float_t position() const
Return the position of the detector plane from which the hit arose.
Definition: StiHit.h:93
StiHit::stHit
const StMeasuredPoint * stHit() const
Definition: StiHit.h:104
Track
C++ STL includes.
Definition: AgUStep.h:47
StiCATpcTrackerInterface::Instance
static StiCATpcTrackerInterface & Instance()
Instance.
Definition: StiCATpcTrackerInterface.cxx:44
Named::getName
const string & getName() const
Get the name of the object.
Definition: Named.cxx:22
Generated by   doxygen 1.8.5