Update of the HLT CA tracker
[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliTPCtrackerCA.cxx
1 // $Id$
2 //***************************************************************************
3 // This file is property of and copyright by the ALICE HLT Project          * 
4 // ALICE Experiment at CERN, All rights reserved.                           *
5 //                                                                          *
6 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
7 //                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
8 //                  for The ALICE HLT Project.                              *
9 //                                                                          *
10 // Permission to use, copy, modify and distribute this software and its     *
11 // documentation strictly for non-commercial purposes is hereby granted     *
12 // without fee, provided that the above copyright notice appears in all     *
13 // copies and that both the copyright notice and this permission notice     *
14 // appear in the supporting documentation. The authors make no claims       *
15 // about the suitability of this software for any purpose. It is            *
16 // provided "as is" without express or implied warranty.                    *
17 //***************************************************************************
18
19 #include "AliTPCtrackerCA.h"
20
21 #include "TTree.h"
22 #include "Riostream.h"
23 //#include "AliCluster.h"
24 #include "AliTPCClustersRow.h"
25 #include "AliTPCParam.h"
26 #include "AliTPCClusterParam.h"
27
28 #include "AliRun.h"
29 #include "AliRunLoader.h"
30 #include "AliStack.h"
31
32 #include "AliHLTTPCCAGBTracker.h"
33 #include "AliHLTTPCCAGBHit.h"
34 #include "AliHLTTPCCAGBTrack.h"
35 #include "AliHLTTPCCAPerformance.h"
36 #include "AliHLTTPCCAParam.h"
37 #include "AliHLTTPCCATrackConvertor.h"
38 #include "AliHLTTPCCATracker.h"
39
40 #include "TMath.h"
41 #include "AliTPCLoader.h"
42 #include "AliTPC.h"
43 #include "AliTPCclusterMI.h"
44 #include "AliTPCTransform.h"
45 #include "AliTPCcalibDB.h"
46 #include "AliTPCtrack.h"
47 #include "AliTPCseed.h"
48 #include "AliESDtrack.h"
49 #include "AliESDEvent.h"
50 #include "AliTrackReference.h"
51 #include "TStopwatch.h"
52
53 //#include <fstream.h>
54
55 ClassImp(AliTPCtrackerCA)
56
57 AliTPCtrackerCA::AliTPCtrackerCA()
58   :AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
59 {
60   //* default constructor
61 }
62
63 AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCtrackerCA &):
64   AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
65 {
66   //* dummy
67 }
68
69 AliTPCtrackerCA & AliTPCtrackerCA::operator=(const AliTPCtrackerCA& )
70 {
71   //* dummy 
72   return *this;
73 }
74
75
76 AliTPCtrackerCA::~AliTPCtrackerCA() 
77 {
78   //* destructor
79   if( fClusters ) delete[] fClusters;
80   if( fHLTTracker ) delete fHLTTracker;
81 }
82
83 //#include "AliHLTTPCCADisplay.h"
84
85 AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCParam *par): 
86   AliTracker(),fParam(par), fClusters(0), fNClusters(0), fHLTTracker(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
87 {
88   //* constructor
89  
90   fDoHLTPerformance = 0;
91   fDoHLTPerformanceClusters = 0;
92
93   fHLTTracker = new AliHLTTPCCAGBTracker;
94   fHLTTracker->SetNSlices( fParam->GetNSector()/2 );
95
96   if( fDoHLTPerformance ){
97     AliHLTTPCCAPerformance::Instance().SetTracker( fHLTTracker );
98   }
99
100   for( Int_t iSlice=0; iSlice<fHLTTracker->NSlices(); iSlice++ ){
101   
102     Float_t bz = AliTracker::GetBz();
103
104     Float_t inRmin = fParam->GetInnerRadiusLow();
105     //Float_t inRmax = fParam->GetInnerRadiusUp();
106     //Float_t outRmin = fParam->GetOuterRadiusLow(); 
107     Float_t outRmax = fParam->GetOuterRadiusUp();
108     Float_t plusZmin = 0.0529937; 
109     Float_t plusZmax = 249.778; 
110     Float_t minusZmin = -249.645; 
111     Float_t minusZmax = -0.0799937; 
112     Float_t dalpha = 0.349066;
113     Float_t alpha = 0.174533 + dalpha*iSlice;
114     
115     Bool_t zPlus = (iSlice<18 );
116     Float_t zMin =  zPlus ?plusZmin :minusZmin;
117     Float_t zMax =  zPlus ?plusZmax :minusZmax;
118     //TPCZmin = -249.645, ZMax = 249.778    
119     //Float_t rMin =  inRmin;
120     //Float_t rMax =  outRmax;
121         
122     Float_t padPitch = 0.4;
123     Float_t sigmaZ = 0.228808;
124
125     Int_t nRows = fParam->GetNRowLow()+fParam->GetNRowUp();
126     Float_t rowX[200];
127     for( Int_t irow=0; irow<fParam->GetNRowLow(); irow++){
128       rowX[irow] = fParam->GetPadRowRadiiLow(irow);
129     }     
130     for( Int_t irow=0; irow<fParam->GetNRowUp(); irow++){
131       rowX[fParam->GetNRowLow()+irow] = fParam->GetPadRowRadiiUp(irow);
132     }
133     AliHLTTPCCAParam param;
134     param.Initialize( iSlice, nRows, rowX, alpha, dalpha,
135                       inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, bz );
136     param.SetHitPickUpFactor( 3. );
137     param.SetMaxTrackMatchDRow( 5 );
138     param.SetTrackConnectionFactor( 3.5 );
139
140     AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
141     for( Int_t iRow=0; iRow<nRows; iRow++ ){
142       Int_t    type = (iRow<63) ? 0: ( (iRow>126) ? 1:2 );
143       for( int iyz=0; iyz<2; iyz++ ){
144         for( int k=0; k<7; k++ ){
145           //std::cout<<param.fParamS0Par[iyz][type][k]<<" "<<clparam->fParamS0Par[iyz][type][k] - param.fParamS0Par[iyz][type][k]<<std::endl;
146           param.SetParamS0Par(iyz,type,k,clparam->fParamS0Par[iyz][type][k]);
147         }
148       }
149     }
150     fHLTTracker->Slices()[iSlice].Initialize( param ); 
151   }
152 }
153
154
155
156 Int_t AliTPCtrackerCA::LoadClusters (TTree * fromTree)
157
158   // load clusters to the local arrays
159   fNClusters = 0;
160   if( fClusters ) delete[] fClusters;
161
162   fHLTTracker->StartEvent();
163   if( fDoHLTPerformance ) AliHLTTPCCAPerformance::Instance().StartEvent();
164
165   if( !fParam ) return 1;
166
167   // load mc tracks
168   while( fDoHLTPerformance ){
169     if( !gAlice ) break;
170     AliRunLoader *rl = AliRunLoader::Instance();//gAlice->GetRunLoader(); 
171     if( !rl ) break;
172     rl->LoadKinematics();
173     AliStack *stack = rl->Stack();
174     if( !stack ) break;
175
176     AliHLTTPCCAPerformance::Instance().SetNMCTracks( stack->GetNtrack() );
177     
178     for( Int_t itr=0; itr<stack->GetNtrack(); itr++ ){
179       TParticle *part = stack->Particle(itr);
180       AliHLTTPCCAPerformance::Instance().ReadMCTrack( itr, part );
181     }
182
183     { // check for MC tracks at the TPC entrance
184
185       Bool_t *isTPC = 0;
186       isTPC = new Bool_t [stack->GetNtrack()];
187       for( Int_t i=0; i<stack->GetNtrack(); i++ ) isTPC[i] = 0;
188       rl->LoadTrackRefs();
189       TTree *mcTree = rl->TreeTR();
190       if( !mcTree ) break;
191       TBranch *branch=mcTree->GetBranch("TrackReferences");
192       if (!branch ) break;      
193       TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs=&tpcdummy;
194       branch->SetAddress(&tpcRefs);
195       Int_t nr=(Int_t)mcTree->GetEntries();
196       for (Int_t r=0; r<nr; r++) {
197         mcTree->GetEvent(r);
198         Int_t nref = tpcRefs->GetEntriesFast();
199         if (!nref) continue;
200         AliTrackReference *tpcRef= 0x0;  
201         for (Int_t iref=0; iref<nref; ++iref) {
202           tpcRef = (AliTrackReference*)tpcRefs->UncheckedAt(iref);
203           if (tpcRef->DetectorId() == AliTrackReference::kTPC) break;
204           tpcRef = 0x0;
205         }
206         if (!tpcRef) continue;
207
208         if( isTPC[tpcRef->Label()] ) continue;  
209
210         AliHLTTPCCAPerformance::Instance().ReadMCTPCTrack(tpcRef->Label(),
211                                         tpcRef->X(),tpcRef->Y(),tpcRef->Z(),
212                                         tpcRef->Px(),tpcRef->Py(),tpcRef->Pz() );
213         isTPC[tpcRef->Label()] = 1;
214         tpcRefs->Clear();
215       } 
216       if( isTPC ) delete[] isTPC;
217     }
218
219     while( fDoHLTPerformanceClusters ){
220       AliTPCLoader *tpcl = (AliTPCLoader*) rl->GetDetectorLoader("TPC");
221       if( !tpcl ) break;
222       if( tpcl->TreeH() == 0x0 ){
223         if( tpcl->LoadHits() ) break;
224       }
225       if( tpcl->TreeH() == 0x0 ) break;
226       
227       AliTPC *tpc = (AliTPC*) gAlice->GetDetector("TPC");
228       Int_t nEnt=(Int_t)tpcl->TreeH()->GetEntries();
229       Int_t nPoints = 0;
230       for (Int_t iEnt=0; iEnt<nEnt; iEnt++) {    
231         tpc->ResetHits();
232         tpcl->TreeH()->GetEvent(iEnt);
233         AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1);
234         for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ) nPoints++;
235       }
236       AliHLTTPCCAPerformance::Instance().SetNMCPoints( nPoints );
237
238       for (Int_t iEnt=0; iEnt<nEnt; iEnt++) {    
239         tpc->ResetHits();
240         tpcl->TreeH()->GetEvent(iEnt);
241         AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1);
242         for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ){
243           AliHLTTPCCAPerformance::Instance().ReadMCPoint( phit->GetTrack(),phit->X(), phit->Y(),phit->Z(),phit->Time(), phit->fSector%36);
244         }      
245       }
246       break;
247     }
248     break;
249   }
250   
251   TBranch * br = fromTree->GetBranch("Segment");
252   if( !br ) return 1;
253
254   AliTPCClustersRow *clrow = new AliTPCClustersRow;
255   clrow->SetClass("AliTPCclusterMI");
256   clrow->SetArray(0);
257   clrow->GetArray()->ExpandCreateFast(10000);
258   
259   br->SetAddress(&clrow);
260   
261   //
262   Int_t nEnt=Int_t(fromTree->GetEntries());
263
264   fNClusters = 0;
265   for (Int_t i=0; i<nEnt; i++) {
266     br->GetEntry(i);
267     Int_t sec,row;
268     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
269     fNClusters += clrow->GetArray()->GetEntriesFast();
270   }
271
272   fClusters = new AliTPCclusterMI [fNClusters];
273   fHLTTracker->SetNHits( fNClusters );
274   if( fDoHLTPerformance ) AliHLTTPCCAPerformance::Instance().SetNHits( fNClusters );
275   Int_t ind=0;
276   for (Int_t i=0; i<nEnt; i++) {
277     br->GetEntry(i);
278     Int_t sec,row;
279     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
280     Int_t nClu = clrow->GetArray()->GetEntriesFast();
281     Float_t x = fParam->GetPadRowRadii(sec,row);
282     for (Int_t icl=0; icl<nClu; icl++){
283       Int_t lab0 = -1;
284       Int_t lab1 = -1;
285       Int_t lab2 = -1;
286       AliTPCclusterMI* cluster = (AliTPCclusterMI*)(clrow->GetArray()->At(icl));
287       if( !cluster ) continue;
288       lab0 = cluster->GetLabel(0);
289       lab1 = cluster->GetLabel(1);
290       lab2 = cluster->GetLabel(2);
291
292       AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
293       if (!transform) {
294         AliFatal("Tranformations not in calibDB");
295       }
296       Double_t xx[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
297       Int_t id[1]={cluster->GetDetector()};
298       transform->Transform(xx,id,0,1);  
299       //if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
300       //if (cluster->GetDetector()%36>17){
301       //xx[1]*=-1;
302       //}
303       //}
304
305       cluster->SetX(xx[0]);
306       cluster->SetY(xx[1]);
307       cluster->SetZ(xx[2]);
308
309       TGeoHMatrix  *mat = fParam->GetClusterMatrix(cluster->GetDetector());
310       Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
311       Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
312       if (mat) mat->LocalToMaster(pos,posC);
313       else{
314         // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
315       }
316       cluster->SetX(posC[0]);
317       cluster->SetY(posC[1]);
318       cluster->SetZ(posC[2]);
319
320       Float_t y = cluster->GetY();
321       Float_t z = cluster->GetZ();        
322
323       if( sec>=36 ){
324         sec = sec - 36;
325         row = row + fParam->GetNRowLow(); 
326       }
327       
328       Int_t index = ind++;
329       fClusters[index] = *cluster;
330       fHLTTracker->ReadHit( x, y, z, 
331                             TMath::Sqrt(cluster->GetSigmaY2()), TMath::Sqrt(cluster->GetSigmaZ2()),                         
332                             cluster->GetMax(), index, sec, row );
333       if( fDoHLTPerformance ) AliHLTTPCCAPerformance::Instance().ReadHitLabel(index, lab0, lab1, lab2 );
334     }
335   }
336   delete clrow;
337   return 0;
338 }
339
340 AliCluster * AliTPCtrackerCA::GetCluster(Int_t index) const
341 {
342   return &(fClusters[index]);
343 }
344
345 Int_t AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event )
346 {
347   // reconstruction
348   //cout<<"Start of AliTPCtrackerCA"<<endl;
349   TStopwatch timer;
350
351   fHLTTracker->FindTracks();
352   //cout<<"Do performance.."<<endl;
353   if( fDoHLTPerformance ) AliHLTTPCCAPerformance::Instance().Performance();
354
355   if( 0 ) {// Write Event    
356     if( fStatNEvents == 0 ){
357       fstream geo;
358       geo.open("CAEvents/settings.dat", ios::out);
359       if( geo.is_open() ){
360         fHLTTracker->WriteSettings(geo);        
361       }
362       geo.close();
363     }
364
365     fstream hits;
366     char name[255];
367     sprintf( name,"CAEvents/%i.event.dat",fStatNEvents ); 
368     hits.open(name, ios::out);
369     if( hits.is_open() ){
370       fHLTTracker->WriteEvent(hits);    
371       fstream tracks;
372       sprintf( name,"CAEvents/%i.tracks.dat",fStatNEvents ); 
373       tracks.open(name, ios::out);
374       fHLTTracker->WriteTracks(tracks); 
375     }
376     hits.close();   
377     if( fDoHLTPerformance ){
378       fstream mcevent, mcpoints;
379       char mcname[255];
380       sprintf( mcname,"CAEvents/%i.mcevent.dat",fStatNEvents ); 
381       mcevent.open(mcname, ios::out);
382       if( mcevent.is_open() ){      
383         AliHLTTPCCAPerformance::Instance().WriteMCEvent(mcevent);
384       }
385       if(1 && fDoHLTPerformanceClusters ){
386         sprintf( mcname,"CAEvents/%i.mcpoints.dat",fStatNEvents ); 
387         mcpoints.open(mcname, ios::out);
388         if( mcpoints.is_open() ){      
389           AliHLTTPCCAPerformance::Instance().WriteMCPoints(mcpoints);
390         }
391         mcpoints.close();
392       }
393       mcevent.close();   
394     }
395   }
396   fStatNEvents++;
397
398   if( event ){
399    
400     Float_t bz = fHLTTracker->Slices()[0].Param().Bz();
401
402     for( Int_t itr=0; itr<fHLTTracker->NTracks(); itr++ ){
403       //AliTPCtrack tTPC;
404       AliTPCseed tTPC;
405       AliHLTTPCCAGBTrack &tCA = fHLTTracker->Tracks()[itr];
406       AliHLTTPCCATrackParam par = tCA.Param();  
407       AliHLTTPCCATrackConvertor::GetExtParam( par, tTPC, tCA.Alpha(), bz );
408       tTPC.SetMass(0.13957);
409       tTPC.SetdEdx( tCA.DeDx() );
410       if( TMath::Abs(tTPC.GetSigned1Pt())>1./0.02 ) continue;
411       Int_t nhits = tCA.NHits();
412       if( nhits>199 ) nhits=199;// kMaxRow ) nhits = kMaxRow;
413       tTPC.SetNumberOfClusters(nhits);
414  
415       Float_t alpha = tCA.Alpha();      
416       AliHLTTPCCATrackParam t0 = par;
417       for( Int_t ih=0; ih<nhits; ih++ ){
418         Int_t index = fHLTTracker->TrackHits()[tCA.FirstHitRef()+ih];
419         AliHLTTPCCAGBHit &h = fHLTTracker->Hits()[index];
420         Int_t extIndex = h.ID();
421         tTPC.SetClusterIndex(ih, extIndex);
422                 
423           AliTPCclusterMI *c = &(fClusters[extIndex]);
424           tTPC.SetClusterPointer(h.IRow(), c );
425           AliTPCTrackerPoint &point = *(tTPC.GetTrackPoint(h.IRow()));
426           {
427             Int_t iSlice = h.ISlice();
428             AliHLTTPCCATracker &slice = fHLTTracker->Slices()[iSlice];
429             if( slice.Param().Alpha()!=alpha ){
430               if( ! t0.Rotate(  slice.Param().Alpha() - alpha, .999 ) ) continue;
431               alpha = slice.Param().Alpha();
432             }
433             Float_t x = slice.Row(h.IRow()).X();
434             if( !t0.TransportToX( x, .999 ) ) continue;
435             Float_t sy2, sz2;
436             slice.GetErrors2( h.IRow(), t0, sy2, sz2 );
437             point.SetSigmaY(c->GetSigmaY2()/sy2);
438             point.SetSigmaZ(c->GetSigmaZ2()/sz2);
439             point.SetAngleY(TMath::Abs(t0.GetSinPhi()/t0.GetCosPhi()));
440             point.SetAngleZ(TMath::Abs(t0.GetDzDs()));  
441           }
442         
443         }
444       tTPC.CookdEdx(0.02,0.6);
445
446       CookLabel(&tTPC,0.1);           
447
448       AliESDtrack tESD;
449       tESD.UpdateTrackParams( &(tTPC),AliESDtrack::kTPCin);
450       //tESD.SetStatus( AliESDtrack::kTPCrefit );
451       //tESD.SetTPCPoints(tTPC.GetPoints());
452       Int_t   ndedx = tTPC.GetNCDEDX(0);
453       Float_t sdedx = tTPC.GetSDEDX(0);
454       Float_t dedx  = tTPC.GetdEdx();
455       tESD.SetTPCsignal(dedx, sdedx, ndedx); 
456       tESD.myTPC = tTPC;            
457
458       event->AddTrack(&tESD);
459     }
460   }
461   timer.Stop();
462   static double time=0, time1 = 0;
463   static Int_t ncalls = 0;
464   time+=timer.CpuTime();
465   time1+=timer.RealTime();
466   ncalls++;
467   //cout<<"\n\nCA tracker speed: cpu = "<<time/ncalls*1.e3<<" [ms/ev], real = "<<time1/ncalls*1.e3<<" [ms/ev], n calls = "<<ncalls<<endl<<endl;
468
469   //cout<<"End of AliTPCtrackerCA"<<endl;
470   return 0;
471 }
472
473
474 Int_t AliTPCtrackerCA::RefitInward (AliESDEvent *event)
475
476   //* back propagation of ESD tracks (not fully functional)
477
478   Float_t bz = fHLTTracker->Slices()[0].Param().Bz();
479   Float_t xTPC = fParam->GetInnerRadiusLow();
480   Float_t dAlpha = fParam->GetInnerAngle()/180.*TMath::Pi();
481   Float_t yMax = xTPC*TMath::Tan(dAlpha/2.); 
482
483   Int_t nentr=event->GetNumberOfTracks();
484      
485   for (Int_t i=0; i<nentr; i++) {
486     AliESDtrack *esd=event->GetTrack(i);
487     ULong_t status=esd->GetStatus(); 
488     if (!(status&AliESDtrack::kTPCin)) continue;
489     AliHLTTPCCATrackParam t0;
490     AliHLTTPCCATrackConvertor::SetExtParam(t0,*esd, bz );
491     Float_t alpha = esd->GetAlpha();
492     if( t0.TransportToXWithMaterial( xTPC, bz) ){
493       if (t0.GetY() > yMax) {
494         if (t0.Rotate(dAlpha)){ 
495           alpha+=dAlpha;  
496           t0.TransportToXWithMaterial( xTPC, bz);
497         }
498       } else if (t0.GetY() <-yMax) {
499         if (t0.Rotate(-dAlpha)){
500           alpha+=-dAlpha;
501           t0.TransportToXWithMaterial( xTPC, bz);
502         }
503       }    
504     }
505     AliTPCtrack tt(*esd);
506     AliHLTTPCCATrackConvertor::GetExtParam(t0,tt,alpha,bz);
507     esd->UpdateTrackParams( &tt,AliESDtrack::kTPCrefit); 
508   }
509   return 0;
510 }
511
512 Int_t AliTPCtrackerCA::PropagateBack(AliESDEvent *)
513
514   //* not implemented yet
515   return 0; 
516 }