unnecessary header file
[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 "AliRun.h"
27 #include "AliRunLoader.h"
28 #include "AliStack.h"
29
30 #include "AliHLTTPCCATracker.h"
31 #include "AliHLTTPCCAGBHit.h"
32 #include "AliHLTTPCCAGBTracker.h"
33 #include "AliHLTTPCCAGBTrack.h"
34 #include "AliHLTTPCCAMCTrack.h"
35 #include "AliHLTTPCCAOutTrack.h"
36 #include "AliHLTTPCCAPerformance.h"
37 #include "AliHLTTPCCAParam.h"
38 #include "AliHLTTPCCATrackConvertor.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 "AliTPCReconstructor.h"
47 #include "AliTPCtrack.h"
48 #include "AliESDtrack.h"
49 #include "AliESDEvent.h"
50 #include "AliTrackReference.h"
51
52 //#include <fstream.h>
53
54 ClassImp(AliTPCtrackerCA)
55
56 AliTPCtrackerCA::AliTPCtrackerCA()
57   :AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
58 {
59   //* default constructor
60 }
61
62 AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCtrackerCA &):
63   AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
64 {
65   //* dummy
66 }
67
68 AliTPCtrackerCA & AliTPCtrackerCA::operator=(const AliTPCtrackerCA& )
69 {
70   //* dummy 
71   return *this;
72 }
73
74
75 AliTPCtrackerCA::~AliTPCtrackerCA() 
76 {
77   //* destructor
78   if( fClusters ) delete[] fClusters;
79   if( fHLTTracker ) delete fHLTTracker;
80   if( fHLTPerformance ) delete fHLTPerformance;
81 }
82
83 AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCParam *par): 
84   AliTracker(),fParam(par), fClusters(0), fNClusters(0), fHLTTracker(0), fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
85 {
86   //* constructor
87
88   DoHLTPerformance() = 1;
89   DoHLTPerformanceClusters() = 1;
90
91   fHLTTracker = new AliHLTTPCCAGBTracker;
92   fHLTTracker->SetNSlices( fParam->GetNSector()/2 );
93
94   if( fDoHLTPerformance ){
95     fHLTPerformance = new AliHLTTPCCAPerformance;
96     fHLTPerformance->SetTracker( fHLTTracker );
97   }
98
99   for( int iSlice=0; iSlice<fHLTTracker->NSlices(); iSlice++ ){
100   
101     Float_t bz = AliTracker::GetBz();
102
103     Float_t inRmin = fParam->GetInnerRadiusLow();
104     //Float_t inRmax = fParam->GetInnerRadiusUp();
105     //Float_t outRmin = fParam->GetOuterRadiusLow(); 
106     Float_t outRmax = fParam->GetOuterRadiusUp();
107     Float_t plusZmin = 0.0529937; 
108     Float_t plusZmax = 249.778; 
109     Float_t minusZmin = -249.645; 
110     Float_t minusZmax = -0.0799937; 
111     Float_t dalpha = 0.349066;
112     Float_t alpha = 0.174533 + dalpha*iSlice;
113     
114     Bool_t zPlus = (iSlice<18 );
115     Float_t zMin =  zPlus ?plusZmin :minusZmin;
116     Float_t zMax =  zPlus ?plusZmax :minusZmax;
117     //TPCZmin = -249.645, ZMax = 249.778    
118     //Float_t rMin =  inRmin;
119     //Float_t rMax =  outRmax;
120         
121     Float_t padPitch = 0.4;
122     Float_t sigmaZ = 0.228808;
123
124     Int_t NRows = fParam->GetNRowLow()+fParam->GetNRowUp();
125
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.YErrorCorrection() = 1.;//.33;
137     param.ZErrorCorrection() = 1.;//.33;
138     param.MaxTrackMatchDRow() = 5;
139     param.TrackConnectionFactor() = 3.5;
140     fHLTTracker->Slices()[iSlice].Initialize( param ); 
141   }
142 }
143
144
145
146 Int_t AliTPCtrackerCA::LoadClusters (TTree * tree)
147
148   fNClusters = 0;
149   if( fClusters ) delete[] fClusters;
150
151   fHLTTracker->StartEvent();
152   if( fDoHLTPerformance ) fHLTPerformance->StartEvent();
153
154   if( !fParam ) return 1;
155
156   // load mc tracks
157   while( fDoHLTPerformance ){
158     if( !gAlice ) break;
159     AliRunLoader *rl = gAlice->GetRunLoader(); 
160     if( !rl ) break;
161     rl->LoadKinematics();
162     AliStack *stack = rl->Stack();
163     if( !stack ) break;
164
165     fHLTPerformance->SetNMCTracks( stack->GetNtrack() );
166     
167     for( Int_t itr=0; itr<stack->GetNtrack(); itr++ ){
168       TParticle *part = stack->Particle(itr);
169       fHLTPerformance->ReadMCTrack( itr, part );
170     }
171
172     { // check for MC tracks at the TPC entrance
173       Bool_t *isTPC = 0;
174       isTPC = new Bool_t [stack->GetNtrack()];
175       for( Int_t i=0; i<stack->GetNtrack(); i++ ) isTPC[i] = 0;
176       rl->LoadTrackRefs();
177       TTree *TR = rl->TreeTR();
178       if( !TR ) break;
179       TBranch *branch=TR->GetBranch("TrackReferences");
180       if (!branch ) break;      
181       TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs=&tpcdummy;
182       branch->SetAddress(&tpcRefs);
183       Int_t nr=(Int_t)TR->GetEntries();
184       for (Int_t r=0; r<nr; r++) {
185         TR->GetEvent(r);
186         Int_t nref = tpcRefs->GetEntriesFast();
187         if (!nref) continue;
188         AliTrackReference *tpcRef= 0x0;  
189         for (Int_t iref=0; iref<nref; ++iref) {
190           tpcRef = (AliTrackReference*)tpcRefs->UncheckedAt(iref);
191           if (tpcRef->DetectorId() == AliTrackReference::kTPC) break;
192           tpcRef = 0x0;
193         }
194         if (!tpcRef) continue;
195
196         if( isTPC[tpcRef->Label()] ) continue;  
197
198         fHLTPerformance->ReadMCTPCTrack(tpcRef->Label(),
199                                         tpcRef->X(),tpcRef->Y(),tpcRef->Z(),
200                                         tpcRef->Px(),tpcRef->Py(),tpcRef->Pz() );
201         isTPC[tpcRef->Label()] = 1;
202         tpcRefs->Clear();
203       } 
204       if( isTPC ) delete[] isTPC;
205     }
206
207     while( fDoHLTPerformanceClusters ){
208       AliTPCLoader *tpcl = (AliTPCLoader*) rl->GetDetectorLoader("TPC");
209       if( !tpcl ) break;
210       if( tpcl->TreeH() == 0x0 ){
211         if( tpcl->LoadHits() ) break;
212       }
213       if( tpcl->TreeH() == 0x0 ) break;
214       
215       AliTPC *tpc = (AliTPC*) gAlice->GetDetector("TPC");
216       Int_t nEnt=(Int_t)tpcl->TreeH()->GetEntries();
217       Int_t nPoints = 0;
218       for (Int_t iEnt=0; iEnt<nEnt; iEnt++) {    
219         tpc->ResetHits();
220         tpcl->TreeH()->GetEvent(iEnt);
221         AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1);
222         for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ) nPoints++;
223       }
224       fHLTPerformance->SetNMCPoints( nPoints );
225
226       for (Int_t iEnt=0; iEnt<nEnt; iEnt++) {    
227         tpc->ResetHits();
228         tpcl->TreeH()->GetEvent(iEnt);
229         AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1);
230         for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ){
231           fHLTPerformance->ReadMCPoint( phit->GetTrack(),phit->X(), phit->Y(),phit->Z(),phit->Time(), phit->fSector%36);
232         }      
233       }
234       break;
235     }
236     break;
237   }
238   
239   TBranch * br = tree->GetBranch("Segment");
240   if( !br ) return 1;
241
242   AliTPCClustersRow *clrow = new AliTPCClustersRow;
243   clrow->SetClass("AliTPCclusterMI");
244   clrow->SetArray(0);
245   clrow->GetArray()->ExpandCreateFast(10000);
246   
247   br->SetAddress(&clrow);
248   
249   //
250   Int_t NEnt=Int_t(tree->GetEntries());
251
252   fNClusters = 0;
253   for (Int_t i=0; i<NEnt; i++) {
254     br->GetEntry(i);
255     Int_t sec,row;
256     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
257     fNClusters += clrow->GetArray()->GetEntriesFast();
258   }
259
260   fClusters = new AliTPCclusterMI [fNClusters];
261   fHLTTracker->SetNHits( fNClusters );
262   if( fDoHLTPerformance ) fHLTPerformance->SetNHits( fNClusters );
263   int ind=0;
264   for (Int_t i=0; i<NEnt; i++) {
265     br->GetEntry(i);
266     Int_t sec,row;
267     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
268     int NClu = clrow->GetArray()->GetEntriesFast();
269     Float_t x = fParam->GetPadRowRadii(sec,row);
270     for (Int_t icl=0; icl<NClu; icl++){
271       Int_t lab0 = -1;
272       Int_t lab1 = -1;
273       Int_t lab2 = -1;
274       AliTPCclusterMI* cluster = (AliTPCclusterMI*)(clrow->GetArray()->At(icl));
275       if( !cluster ) continue;
276       lab0 = cluster->GetLabel(0);
277       lab1 = cluster->GetLabel(1);
278       lab2 = cluster->GetLabel(2);
279
280       AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
281       if (!transform) {
282         AliFatal("Tranformations not in calibDB");
283       }
284       Double_t xx[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
285       Int_t id[1]={cluster->GetDetector()};
286       transform->Transform(xx,id,0,1);  
287       //if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
288       if (cluster->GetDetector()%36>17){
289         xx[1]*=-1;
290       }
291       //}
292
293       cluster->SetX(xx[0]);
294       cluster->SetY(xx[1]);
295       cluster->SetZ(xx[2]);
296
297       TGeoHMatrix  *mat = fParam->GetClusterMatrix(cluster->GetDetector());
298       Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
299       Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
300       if (mat) mat->LocalToMaster(pos,posC);
301       else{
302         // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
303       }
304       cluster->SetX(posC[0]);
305       cluster->SetY(posC[1]);
306       cluster->SetZ(posC[2]);
307
308       Float_t y = cluster->GetY();
309       Float_t z = cluster->GetZ();        
310
311       if( sec>=36 ){
312         sec = sec - 36;
313         row = row + fParam->GetNRowLow(); 
314       }
315       
316       Int_t index = ind++;
317       fClusters[index] = *cluster;
318       fHLTTracker->ReadHit( x, y, z, 
319                             TMath::Sqrt(cluster->GetSigmaY2()), TMath::Sqrt(cluster->GetSigmaZ2()),                         
320                             cluster->GetMax(), index, sec, row );
321       if( fDoHLTPerformance ) fHLTPerformance->ReadHitLabel(index, lab0, lab1, lab2 );
322     }
323   }
324   delete clrow;
325   return 0;
326 }
327
328 AliCluster * AliTPCtrackerCA::GetCluster(Int_t index) const
329 {
330   return &(fClusters[index]);
331 }
332
333 Int_t AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event )
334 {
335   //cout<<"Start of AliTPCtrackerCA"<<endl;
336
337   fHLTTracker->FindTracks();
338   if( fDoHLTPerformance ) fHLTPerformance->Performance();
339
340   if( 0 ) {// Write Event    
341     if( fStatNEvents == 0 ){
342       fstream geo;
343       geo.open("CAEvents/settings.dat", ios::out);
344       if( geo.is_open() ){
345         fHLTTracker->WriteSettings(geo);        
346       }
347       geo.close();
348     }
349
350     fstream event;
351     char name[255];
352     sprintf( name,"CAEvents/%i.event.dat",fStatNEvents ); 
353     event.open(name, ios::out);
354     if( event.is_open() ){
355       fHLTTracker->WriteEvent(event);   
356       fstream tracks;
357       sprintf( name,"CAEvents/%i.tracks.dat",fStatNEvents ); 
358       tracks.open(name, ios::out);
359       fHLTTracker->WriteTracks(tracks); 
360     }
361     event.close();   
362     if( fDoHLTPerformance ){
363       fstream mcevent, mcpoints;
364       char mcname[255];
365       sprintf( mcname,"CAEvents/%i.mcevent.dat",fStatNEvents ); 
366       mcevent.open(mcname, ios::out);
367       if( mcevent.is_open() ){      
368         fHLTPerformance->WriteMCEvent(mcevent);
369       }
370       if(fDoHLTPerformanceClusters ){
371         sprintf( mcname,"CAEvents/%i.mcpoints.dat",fStatNEvents ); 
372         mcpoints.open(mcname, ios::out);
373         if( mcpoints.is_open() ){      
374           fHLTPerformance->WriteMCPoints(mcpoints);
375         }
376         mcpoints.close();
377       }
378       mcevent.close();   
379     }
380   }
381   fStatNEvents++;
382
383   if( event ){
384    
385     Float_t bz = fHLTTracker->Slices()[0].Param().Bz();
386
387     for( Int_t itr=0; itr<fHLTTracker->NTracks(); itr++ ){
388       AliTPCtrack tTPC;
389       AliHLTTPCCAGBTrack &tCA = fHLTTracker->Tracks()[itr];
390       AliHLTTPCCATrackParam &par = tCA.Param(); 
391       AliHLTTPCCATrackConvertor::GetExtParam( par, tTPC, tCA.Alpha(), bz );
392       tTPC.SetMass(0.13957);
393       tTPC.SetdEdx( tCA.DeDx() );
394       if( TMath::Abs(tTPC.GetSigned1Pt())>1./0.1 ) continue;
395       int nhits = tCA.NHits();
396       if( nhits>kMaxRow ) nhits = kMaxRow;
397       tTPC.SetNumberOfClusters(nhits);      
398       for( Int_t ih=0; ih<nhits; ih++ ){
399         Int_t index = fHLTTracker->TrackHits()[tCA.FirstHitRef()+ih];
400         Int_t ext_index = fHLTTracker->Hits()[index].ID();
401         tTPC.SetClusterIndex(ih, ext_index);
402       }
403       CookLabel(&tTPC,0.1);           
404       {
405         Double_t xTPC=fParam->GetInnerRadiusLow();
406         Double_t dAlpha = fParam->GetInnerAngle()/180.*TMath::Pi();
407         if (tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz)) {   
408           Double_t y=tTPC.GetY();
409           Double_t ymax=xTPC*TMath::Tan(dAlpha/2.); 
410           if (y > ymax) {
411             if (tTPC.Rotate(dAlpha)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz);
412           } else if (y <-ymax) {
413             if (tTPC.Rotate(-dAlpha)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz);
414           }         
415         }
416       }
417
418       AliESDtrack tESD;
419       tESD.UpdateTrackParams( &(tTPC),AliESDtrack::kTPCin);
420       //tESD.SetStatus( AliESDtrack::kTPCrefit );
421       //tESD.SetTPCPoints(tTPC.GetPoints());
422       //tESD.myTPC = tTPC;            
423       event->AddTrack(&tESD);
424     }
425   }
426
427   //cout<<"End of AliTPCtrackerCA"<<endl;
428   return 0;
429 }
430
431
432 Int_t AliTPCtrackerCA::RefitInward (AliESDEvent *event)
433
434   //* back propagation of ESD tracks (not fully functional)
435
436   Float_t bz = fHLTTracker->Slices()[0].Param().Bz();
437   Float_t xTPC = fParam->GetInnerRadiusLow();
438   Float_t dAlpha = fParam->GetInnerAngle()/180.*TMath::Pi();
439   Float_t yMax = xTPC*TMath::Tan(dAlpha/2.); 
440
441   Int_t nentr=event->GetNumberOfTracks();
442      
443   for (Int_t i=0; i<nentr; i++) {
444     AliESDtrack *esd=event->GetTrack(i);
445     ULong_t status=esd->GetStatus();
446     if (!(status&AliESDtrack::kTPCin)) continue;
447     AliHLTTPCCATrackParam t0;
448     AliHLTTPCCATrackConvertor::SetExtParam(t0,*esd, bz );
449     Float_t alpha = esd->GetAlpha();
450     if( t0.TransportToXWithMaterial( xTPC, bz) ){
451       if (t0.GetY() > yMax) {
452         if (t0.Rotate(dAlpha)){ 
453           alpha+=dAlpha;  
454           t0.TransportToXWithMaterial( xTPC, bz);
455         }
456       } else if (t0.GetY() <-yMax) {
457         if (t0.Rotate(-dAlpha)){
458           alpha+=-dAlpha;
459           t0.TransportToXWithMaterial( xTPC, bz);
460         }
461       }    
462     }
463     AliTPCtrack tt(*esd);
464     AliHLTTPCCATrackConvertor::GetExtParam(t0,tt,alpha,bz);
465     esd->UpdateTrackParams( &tt,AliESDtrack::kTPCrefit); 
466   }
467   return 0;
468 }
469
470 Int_t AliTPCtrackerCA::PropagateBack(AliESDEvent *)
471
472   //* not implemented yet
473   return 0; 
474 }