cfab694cab512e2fe857334a97ebd018bc3f488b
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCMerger.cxx
1 //$Id$
2 // Original: AliL3Merger.cxx,v 1.16 2005/06/14 10:55:21 cvetan 
3
4 // Author: Uli Frankenfeld <mailto:franken@fi.uib.no>
5 //*-- Copyright &copy Uli 
6
7 #include "AliHLTTPCLogging.h"
8 #include "AliHLTTPCMerger.h"
9 #include "AliHLTTPCTrack.h"
10 #include "AliHLTTPCTrackSegmentData.h"
11 #include "AliHLTTPCTransform.h"
12 #include "AliHLTTPCTrackArray.h"
13
14 #ifdef use_root //use root ntuple for slow merge
15 #include <TNtuple.h>
16 #include <TTree.h>
17 #include <TFile.h>
18 #endif
19
20 /** \class AliHLTTPCMerger
21 <pre>
22 //_____________________________________________________________
23 // AliHLTTPCMerger
24 //
25 // The HLTTPC merger base class
26 //
27 </pre>
28 */
29
30 #if __GNUC__ >= 3
31 using namespace std;
32 #endif
33
34 ClassImp(AliHLTTPCMerger)
35
36 AliHLTTPCMerger::AliHLTTPCMerger()
37 {
38   //Default constructor
39   fInTrack=0;
40   fOutTrack=0;
41   fCurrentTracks=0;
42   fNIn=0;
43 }
44
45 AliHLTTPCMerger::~AliHLTTPCMerger()
46 {
47   //Destructor
48   DeleteArray();
49 }
50
51 void AliHLTTPCMerger::InitMerger(Int_t ntrackarrays,Char_t *tracktype)
52 {
53   //Used to setup all arrays
54   
55   if(strcmp(tracktype,"AliHLTTPCTrack")==0) fTrackType='t';
56   else if(strcmp(tracktype,"AliHLTTPCConfMapTrack")==0) fTrackType='c';
57 #ifdef INCLUDE_TPC_HOUGH
58   else if(strcmp(tracktype,"AliHLTTPCHoughTrack")==0) fTrackType='h';
59 #endif
60   else
61     LOG(AliHLTTPCLog::kError,"AliHLTTPCMerger::AliHLTTPCMerger","Track types")
62       <<"Unknown tracktype"<<ENDLOG;
63   SetArray(ntrackarrays);
64   fCurrentTracks=0;
65
66 }
67
68 void AliHLTTPCMerger::DeleteArray()
69 {
70   //delete arrays
71   for(Int_t i=0; i<fNIn;i++)
72     {
73       if(!fInTrack[i]) continue;
74       delete fInTrack[i];
75       fInTrack[i]=0;
76     }
77   if(fInTrack)
78     delete[] fInTrack;
79   if(fOutTrack)
80     delete fOutTrack;
81   fInTrack=0;
82   fOutTrack=0;
83 }
84
85 void AliHLTTPCMerger::SetArray(Int_t nin)
86
87   //set arrays
88   DeleteArray();//Make sure arrays are cleaned 
89   
90   fNIn = nin;
91   fInTrack = new AliHLTTPCTrackArray*[fNIn];
92   for(Int_t i=0; i<fNIn;i++)
93     {
94 #ifdef INCLUDE_TPC_HOUGH
95       if(fTrackType=='h')
96         fInTrack[i] = new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
97       else
98 #endif
99         fInTrack[i] = new AliHLTTPCTrackArray("AliHLTTPCTrack");
100       
101     }
102 #ifdef INCLUDE_TPC_HOUGH
103   if(fTrackType=='h')
104     fOutTrack= new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
105   else
106 #endif
107     fOutTrack= new AliHLTTPCTrackArray("AliHLTTPCTrack");
108 }
109
110 void AliHLTTPCMerger::Reset()
111
112   //reset
113   for(Int_t i=0; i<fNIn;i++)
114     {
115       fInTrack[i]->Reset();
116     }
117   fOutTrack->Reset();
118 }
119
120 void AliHLTTPCMerger::FillTracks(Int_t ntracks, AliHLTTPCTrackSegmentData* tr)
121 {
122   //Read tracks from shared memory (or memory)
123   AliHLTTPCTrackArray *destination = GetInTracks(fCurrentTracks);
124   if(Is2Global())
125     destination->FillTracks(ntracks, tr, fSlice);
126   else
127     destination->FillTracks(ntracks, tr);
128 }
129
130 void AliHLTTPCMerger::AddAllTracks()
131
132   //add all tracks
133   for(Int_t i=0; i<GetNIn();i++)
134     {
135       AliHLTTPCTrackArray *in = GetInTracks(i);
136       AliHLTTPCTrackArray *out = GetOutTracks();
137       out->AddTracks(in);
138     }
139 }
140
141 void AliHLTTPCMerger::SortGlobalTracks(AliHLTTPCTrack **tracks, Int_t ntrack)
142
143   //sort global tracks
144   AliHLTTPCTrack **tmp = new AliHLTTPCTrack*[ntrack]; 
145   for(Int_t i=0;i<ntrack;i++) tmp[i] = tracks[i];
146   Int_t *t = new Int_t[ntrack];
147   for(Int_t i=0;i<ntrack;i++) t[i]=-1;
148   
149   for(Int_t j=0;j<ntrack;j++)
150     {
151       Double_t minr=300;
152       Int_t    mini=0;
153       for(Int_t i=0;i<ntrack;i++)
154         {
155           if(!tracks[i]) continue;
156           Double_t rr=pow(tracks[i]->GetFirstPointX(),2)+pow(tracks[i]->GetFirstPointY(),2);
157           Double_t r=sqrt(rr);
158           if(r<minr){
159             minr=r;
160             mini=i;
161           }
162         }
163       t[j]=mini;
164       tracks[mini]=0;
165     }
166   for(Int_t i=0;i<ntrack;i++) tracks[i] = tmp[t[i]];
167   delete[] t;
168   delete[] tmp;
169 }
170
171 void AliHLTTPCMerger::SortTracks(AliHLTTPCTrack **tracks, Int_t ntrack) const
172
173   //sort tracks
174   AliHLTTPCTrack **tmp = new  AliHLTTPCTrack*[ntrack];
175   for(Int_t i=0;i<ntrack;i++) tmp[i] = tracks[i];
176   Int_t *t = new Int_t[ntrack];
177   for(Int_t i=0;i<ntrack;i++) t[i]=-1;
178   
179   for(Int_t j=0;j<ntrack;j++)
180     {
181       Double_t minx=300; 
182       Int_t    mini=0;
183       for(Int_t i=0;i<ntrack;i++)
184         {
185           if(!tracks[i]) continue;
186           if(tracks[i]->GetFirstPointX()<minx)
187             {
188               minx=tracks[i]->GetFirstPointX();
189               mini=i;
190             }     
191         }
192       t[j]=mini;  
193       tracks[mini]=0;
194     }
195   for(Int_t i=0;i<ntrack;i++) tracks[i] = tmp[t[i]];
196   delete[] t;
197   delete[] tmp;
198 }
199
200 void AliHLTTPCMerger::AddTrack(AliHLTTPCTrackArray *mergedtrack,AliHLTTPCTrack *track)
201
202   // add tracks
203   AliHLTTPCTrack *t[1];
204   t[0] = track;
205   MultiMerge(mergedtrack,t,1);
206 }
207
208 AliHLTTPCTrack * AliHLTTPCMerger::MergeTracks(AliHLTTPCTrackArray *mergedtrack,AliHLTTPCTrack *t0,AliHLTTPCTrack *t1)
209
210   //merge tracks
211   AliHLTTPCTrack *t[2];
212   t[0] = t0; 
213   t[1] = t1;
214   SortTracks(t,2);
215   return MultiMerge(mergedtrack,t,2);
216 }
217
218 AliHLTTPCTrack * AliHLTTPCMerger::MultiMerge(AliHLTTPCTrackArray *mergedtracks,AliHLTTPCTrack **tracks, Int_t ntrack)
219 {
220   //multi merge the tracks
221   //check npoints
222   Int_t nps = 0;
223   for(Int_t i=0;i<ntrack;i++)
224     {
225       nps+=tracks[i]->GetNHits();
226     }
227   if(nps>AliHLTTPCTransform::GetNRows())
228     {
229       LOG(AliHLTTPCLog::kWarning,"AliHLTTPCMerger::MultiMerge","Adding Points")
230         <<AliHLTTPCLog::kDec<<"Too many Points: "<<nps<<ENDLOG;
231       return 0;
232     }
233   
234   //create new track
235   AliHLTTPCTrack *newtrack = mergedtracks->NextTrack();
236   //copy points
237   UInt_t * nn = new UInt_t[AliHLTTPCTransform::GetNRows()];
238   nps = 0;
239   
240   //  for(Int_t i=0;i<ntrack;i++){
241   for(Int_t i=ntrack-1;i>=0;i--)
242     {
243       memcpy(&nn[nps],tracks[i]->GetHitNumbers(),tracks[i]->GetNHits()*sizeof(UInt_t));
244       nps+=tracks[i]->GetNHits();
245     }
246   AliHLTTPCTrack *tpf=tracks[0];
247   AliHLTTPCTrack *tpl=tracks[ntrack-1];
248   AliHLTTPCTrack *best = tpf;
249   if(tpf->GetNHits()<tpl->GetNHits() && Is2Global())
250     best = tpl;//Best means = most points and therefore best fit (in global case)
251   
252   newtrack->SetNHits(nps);
253   newtrack->SetHits(nps,nn);
254   newtrack->SetFirstPoint(tpf->GetFirstPointX(),tpf->GetFirstPointY(),tpf->GetFirstPointZ());
255   newtrack->SetLastPoint(tpl->GetLastPointX(),tpl->GetLastPointY(),tpl->GetLastPointZ());
256   newtrack->SetPt(best->GetPt());
257   newtrack->SetPsi(best->GetPsi());
258   newtrack->SetTgl(best->GetTgl());
259   newtrack->SetCharge(tpf->GetCharge());
260   delete [] nn;
261   return newtrack;
262 }
263
264 void* AliHLTTPCMerger::GetNtuple(char *varlist) const
265
266   //get ntuple
267 #ifdef use_root
268   TNtuple* nt = new TNtuple("ntuple","ntuple",varlist);
269   return (void*) nt;
270 #else
271   return 0;
272 #endif
273 }
274
275 void* AliHLTTPCMerger::GetNtuple() const
276
277   //get ntuple
278 #ifdef use_root
279   TNtuple* nt = new TNtuple("ntuple","ntuple",
280                             "dx:dy:dz:dk:dpsi:dtgl:dq:disx:disy:disz:dis:n0:n1:diff:drx:dry:drz");
281   return (void*) nt;
282 #else
283   return 0;
284 #endif
285 }
286
287 Bool_t AliHLTTPCMerger::WriteNtuple(char *filename, void* nt) const
288
289   //write ntuple
290 #ifdef use_root
291   TNtuple *ntuple=(TNtuple *) nt;
292   TFile *f = new TFile(filename,"RECREATE");
293   ntuple->Write();
294   f->Close();
295   delete ntuple; 
296   return kTRUE; 
297 #else
298   return kFALSE;
299 #endif
300 }
301
302 void AliHLTTPCMerger::FillNtuple(void *nt,AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
303
304   //fill ntuple
305   Float_t data[17];
306   if(outertrack->IsPoint()&&innertrack->IsPoint())
307     {
308       data[0] =Float_t(innertrack->GetPointX()-outertrack->GetPointX());
309       data[1] =Float_t(innertrack->GetPointY()-outertrack->GetPointY());
310       data[2] =Float_t(innertrack->GetPointZ()-outertrack->GetPointZ());
311       data[3] =Float_t(innertrack->GetKappa()-outertrack->GetKappa());
312       Double_t psi= innertrack->GetPointPsi() - outertrack->GetPointPsi();
313       if(psi>AliHLTTPCTransform::Pi()) psi-=AliHLTTPCTransform::TwoPi();
314       else if(psi<-AliHLTTPCTransform::Pi()) psi+=AliHLTTPCTransform::TwoPi();
315       data[4] =Float_t(psi);
316       data[5] =Float_t(innertrack->GetTgl()-outertrack->GetTgl());
317       data[6] =Float_t(innertrack->GetCharge()-outertrack->GetCharge());
318       data[7] =Float_t(innertrack->GetLastPointX()-outertrack->GetFirstPointX());
319       data[8] =Float_t(innertrack->GetLastPointY()-outertrack->GetFirstPointY());
320       data[9] =Float_t(innertrack->GetLastPointZ()-outertrack->GetFirstPointZ());
321       data[10] =sqrt(pow(data[7],2)+pow(data[8],2)+pow(data[9],2));
322       data[11]= outertrack->GetNHits();
323       data[12]= innertrack->GetNHits();
324       data[13] = Float_t(TrackDiff(innertrack,outertrack));
325       data[14]=0;
326       data[15]=0;
327       data[16]=0;
328 #ifdef use_root
329       TNtuple *ntuple = (TNtuple *) nt;
330       ntuple->Fill(data);
331 #endif
332     }
333 }
334
335 void AliHLTTPCMerger::FillNtuple(void *nt,Float_t *data) const
336
337   //fill ntuple
338 #ifdef use_root
339   TNtuple *ntuple = (TNtuple *) nt;
340   ntuple->Fill(data);
341 #endif
342 }
343
344 Double_t AliHLTTPCMerger::GetAngle(Double_t a1,Double_t a2)
345
346   //get angle
347   Double_t da = a1 - a2 + 4*AliHLTTPCTransform::Pi();
348   da = fmod(da,AliHLTTPCTransform::TwoPi());
349   if(da>AliHLTTPCTransform::Pi()) da = AliHLTTPCTransform::TwoPi()-da;
350   return da;
351 }
352
353 void AliHLTTPCMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxkappa, Double_t maxpsi, Double_t maxtgl)
354
355   //set parameters for merger
356   fMaxY = maxy;
357   fMaxZ = maxz;
358   fMaxKappa = maxkappa;
359   fMaxPsi = maxpsi;
360   fMaxTgl = maxtgl;
361 }
362
363 Bool_t AliHLTTPCMerger::IsTrack(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
364 {
365   //is track to be merged
366   if(innertrack->GetCharge()!=outertrack->GetCharge()) return kFALSE;
367   if( (!innertrack->IsPoint()) || (!outertrack->IsPoint()) )  return kFALSE; 
368   if(innertrack->GetNHits()+outertrack->GetNHits()>AliHLTTPCTransform::GetNRows()) return kFALSE;
369   
370   if(fabs(innertrack->GetPointY()-outertrack->GetPointY()) >fMaxY) return kFALSE;
371   if(fabs(innertrack->GetPointZ()-outertrack->GetPointZ()) >fMaxZ) return kFALSE;
372   if(fabs(innertrack->GetKappa()-outertrack->GetKappa())   >fMaxKappa) return kFALSE;
373   if(GetAngle(innertrack->GetPointPsi(),outertrack->GetPointPsi()) >fMaxPsi) return kFALSE;
374   if(fabs(innertrack->GetTgl()-outertrack->GetTgl()) >fMaxTgl) return kFALSE;
375   //if no rejection up to this point: merge!!
376   return kTRUE;
377 }
378
379 Bool_t AliHLTTPCMerger::IsRTrack(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
380
381   //same as IsTrack
382   return IsTrack(innertrack,outertrack);
383 }
384
385 Double_t AliHLTTPCMerger::TrackDiff(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
386
387   //return track difference
388   Double_t diff =-1;
389   Double_t x[4],y[4],z[4],dy[4],dz[4];
390   AliHLTTPCTrack *tracks[2]; 
391   
392   tracks[0] = innertrack;
393   tracks[1] = outertrack;
394   SortGlobalTracks(tracks,2);
395   innertrack = tracks[0]; 
396   outertrack = tracks[1];
397   
398   x[0] = innertrack->GetFirstPointX();
399   x[1] = innertrack->GetLastPointX();
400   x[2] = outertrack->GetFirstPointX();
401   x[3] = outertrack->GetLastPointX();
402   
403   y[0] = innertrack->GetFirstPointY();
404   y[1] = innertrack->GetLastPointY();
405   y[2] = outertrack->GetFirstPointY();
406   y[3] = outertrack->GetLastPointY();
407
408   z[0] = innertrack->GetFirstPointZ();
409   z[1] = innertrack->GetLastPointZ();
410   z[2] = outertrack->GetFirstPointZ();
411   z[3] = outertrack->GetLastPointZ();
412
413   
414   outertrack->CalculatePoint(x[0]);
415   if(!outertrack->IsPoint()) return diff;
416   dy[0] = fabs(y[0] - outertrack->GetPointY());
417   dz[0] = fabs(z[0] - outertrack->GetPointZ());
418   
419   outertrack->CalculatePoint(x[1]);
420   if(!outertrack->IsPoint()) return diff;
421   dy[1] = fabs(y[1] - outertrack->GetPointY());
422   dz[1] = fabs(z[1] - outertrack->GetPointZ());
423   
424   innertrack->CalculatePoint(x[2]);
425   if(!innertrack->IsPoint()) return diff;
426   dy[2] = fabs(y[2] - innertrack->GetPointY());
427   dz[2] = fabs(z[2] - innertrack->GetPointZ());
428   
429   innertrack->CalculatePoint(x[3]);
430   if(!innertrack->IsPoint()) return diff;
431   dy[3] = fabs(y[3] - innertrack->GetPointY());
432   dz[3] = fabs(z[3] - innertrack->GetPointZ());
433
434   diff=0;
435   for(Int_t i=0;i<4;i++)
436     diff+=sqrt(dy[i]*dy[i]+dz[i]*dz[i]);
437   return diff; 
438 }
439
440 void AliHLTTPCMerger::PrintDiff(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
441
442   // print difference
443   if(!innertrack->IsPoint()||!outertrack->IsPoint())
444     {
445       LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::PrintDiff","No Points")<<ENDLOG;
446       //cerr<<"AliHLTTPCMerger::PrintDiff: No Points"<<endl;
447       //cerr<<"---------------------------"<<endl;
448       return;
449     } 
450   
451   Double_t dx = innertrack->GetPointX()-outertrack->GetPointX();
452   Double_t dy = innertrack->GetPointY()-outertrack->GetPointY();
453   Double_t dz = innertrack->GetPointZ()-outertrack->GetPointZ();
454   Double_t dk = innertrack->GetKappa()-outertrack->GetKappa();
455   Double_t dpsi= innertrack->GetPointPsi() - outertrack->GetPointPsi();
456   if(dpsi>AliHLTTPCTransform::Pi()) dpsi-=AliHLTTPCTransform::TwoPi();
457   else if(dpsi<-AliHLTTPCTransform::Pi())dpsi+=AliHLTTPCTransform::TwoPi();
458   //Double_t dpsi = GetAngle(innertrack->GetPointPsi(),outertrack->GetPointPsi());
459   Double_t dtgl= innertrack->GetTgl()-outertrack->GetTgl();
460   Double_t dq =innertrack->GetCharge()-outertrack->GetCharge();
461   
462   LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::PrintDiff","Points") <<"dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<" dk: "<<dk<<" dpsi: "<<dpsi<<" dtgl: "<<dtgl<<" dq: "<<dq<<ENDLOG;
463   //fprintf(stderr,"dx: %4f dy: %4f dz: %4f dk: %4f dpsi: %4f dtgl: %4f dq: %4f\n",dx,dy,dz,dk,dpsi,dtgl,dq);
464   //cerr<<"---------------------------"<<endl;
465   
466 }
467
468 void AliHLTTPCMerger::Print()
469 {
470   // print some infos
471   for(Int_t i=0; i<fNIn; i++)
472     {
473       AliHLTTPCTrackArray *ttt= GetInTracks(i);
474       for(Int_t j =0;j<ttt->GetNTracks();j++)
475         {
476           AliHLTTPCTrack *track=ttt->GetCheckedTrack(j);
477           if(!track) continue;
478           track->CalculateHelix();
479           //      Double_t angle = atan2(track->GetLastPointY(),track->GetLastPointX());
480           //      if(angle<0) angle+=AliHLTTPCTransform::Pi();
481           if(track->CalculatePoint(135))
482             //      if(!track->CalculateEdgePoint(angle)) cerr<<"**************"<<endl;     
483             //      if(track->CalculatePoint(track->GetLastPointX()))
484             //      if(track->CalculatePoint(0))
485             {
486               //      PrintTrack(track);
487               //      track->CalculateReferencePoint(AliHLTTPCTransform::Pi()/180.);
488               track->CalculateReferencePoint(0.001);
489               Float_t dx=(float)track->GetPointX()-track->GetPointX();
490               Float_t dy=(float)track->GetPointY()-track->GetPointY();
491               Float_t dz=(float)track->GetPointZ()-track->GetPointZ();
492               LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::Print","RefPoint") <<"npt: "<<track->GetNHits()<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<ENDLOG;
493               
494               //fprintf(stderr,"npt: %3d dx: %8.5f dy: %8.5f dz: %8.5f\n",track->GetNHits(),dx,dy,dz);
495               //cerr<<"---------------------------"<<endl;
496             }
497         }  
498     }
499 }
500
501 void AliHLTTPCMerger::PrintTrack(AliHLTTPCTrack *track)
502
503   //print track info
504   fprintf(stderr,"npt: %3d pt: %.2f psi: %.2f tgl: %5.2f q: %2d\n",
505           track->GetNHits(),track->GetPt(),track->GetPsi(),
506           track->GetTgl(),track->GetCharge());
507   fprintf(stderr,
508           "x1: %6.2f y1: %6.2f z1: %6.2f xl: %6.2f yl: %6.2f zl: %6.2f\n",
509           track->GetFirstPointX(),track->GetFirstPointY(),track->GetFirstPointZ(),
510           track->GetLastPointX(),track->GetLastPointY(),track->GetLastPointZ());
511   if(track->IsPoint())
512     {
513       fprintf(stderr,
514               "R: %.2f Xc: %.2f Yc: %.2f Xp: %.2f Yp: %.2f Zp: %.2f Psip: %.2f\n",
515               track->GetRadius(),track->GetCenterX(),track->GetCenterY(),
516               track->GetPointX(),track->GetPointY(),track->GetPointZ(),
517               track->GetPointPsi());
518     }
519 }