When Pt is bad defined (ex. no field), the multiple scattering effect is calculated...
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCGlobalMerger.cxx
1 // @(#) $Id$
2 // Original: AliHLTGlobalMerger.cxx,v 1.13 2005/06/14 10:55:21 cvetan 
3
4 /**************************************************************************
5  * This file is property of and copyright by the ALICE HLT Project        * 
6  * ALICE Experiment at CERN, All rights reserved.                         *
7  *                                                                        *
8  * Primary Authors: Uli Frankenfeld, maintained by                          *
9  *                  Matthias Richter <Matthias.Richter@ift.uib.no>        *
10  *                  for The ALICE HLT Project.                            *
11  *                                                                        *
12  * Permission to use, copy, modify and distribute this software and its   *
13  * documentation strictly for non-commercial purposes is hereby granted   *
14  * without fee, provided that the above copyright notice appears in all   *
15  * copies and that both the copyright notice and this permission notice   *
16  * appear in the supporting documentation. The authors make no claims     *
17  * about the suitability of this software for any purpose. It is          *
18  * provided "as is" without express or implied warranty.                  *
19  **************************************************************************/
20
21 /** @file   AliHLTTPCGlobalMerger.cxx
22     @author Uli Frankenfeld, maintained by Matthias Richter
23     @date   
24     @brief  The HLT TPC slice merger
25 */
26
27 #include "AliHLTTPCLogging.h"
28 #include "AliHLTTPCGlobalMerger.h"
29 #include "AliHLTTPCTrack.h"
30 #include "AliHLTTPCTransform.h"
31 #include "AliHLTTPCTrackArray.h"
32
33 #if __GNUC__ >= 3
34 using namespace std;
35 #endif
36
37 ClassImp(AliHLTTPCGlobalMerger)
38
39 AliHLTTPCGlobalMerger::AliHLTTPCGlobalMerger()
40   :
41   fNSlices(0),
42   fFirst(0),
43   fLast(0)
44 {
45   //Default constructor. Use Setup to specify and setup the necessary parameters and arrays.
46   Is2Global(kTRUE);
47   SetParameter(0,0,0,0,0);
48 }
49
50
51 AliHLTTPCGlobalMerger::~AliHLTTPCGlobalMerger()
52 {
53   //Destructor
54 }
55
56 void AliHLTTPCGlobalMerger::Setup(Int_t first,Int_t last)
57 {
58   //Used to setup the arrays and everything
59   
60   fNSlices = last-first+1;
61   fFirst = first;
62   fLast = last;
63   InitMerger(last-first+1);
64 }
65
66 void AliHLTTPCGlobalMerger::InitSlice(Int_t slice)
67 {
68   // 
69   // Select Sector The following FillTracks call will 
70   // fill this Sector
71   //
72   fSlice = slice;
73   fCurrentTracks = fSlice - fFirst; 
74 }
75
76 Double_t AliHLTTPCGlobalMerger::CheckTracks(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack,Int_t slice)
77 {
78   //Compare the tracks by propagating the outermost track to the last and first point plane
79   //of the innermost track. This plane is defined by the padrow plane where these points
80   //are.
81   
82   if(innertrack->GetCharge()!=outertrack->GetCharge()) return -1;
83   
84   Float_t angle = 0;//perpendicular to the padrowplane (in local system)
85   AliHLTTPCTransform::Local2GlobalAngle(&angle,slice);
86   Double_t dx[2],dy[2],dz[2];
87   Double_t diff =-1;
88   AliHLTTPCTrack *tracks[2];
89   tracks[0] = innertrack;
90   tracks[1] = outertrack;
91   SortGlobalTracks(tracks,2);
92   innertrack = tracks[0]; 
93   outertrack = tracks[1];
94   
95   Float_t point[3];
96   
97   point[0]=innertrack->GetLastPointX();
98   point[1]=innertrack->GetLastPointY();
99   point[2]=innertrack->GetLastPointZ();
100   AliHLTTPCTransform::Global2LocHLT(point,slice);
101   
102   outertrack->CalculateReferencePoint(angle,point[0]);//local x = global distance to padrowplane
103   if(!outertrack->IsPoint()) return diff;
104   dx[0] = fabs(outertrack->GetPointX()-innertrack->GetLastPointX());
105   dy[0] = fabs(outertrack->GetPointY()-innertrack->GetLastPointY());
106   dz[0] = fabs(outertrack->GetPointZ()-innertrack->GetLastPointZ());
107   
108   point[0]=innertrack->GetFirstPointX();
109   point[1]=innertrack->GetFirstPointY();
110   point[2]=innertrack->GetFirstPointZ();
111   AliHLTTPCTransform::Global2LocHLT(point,slice);
112   
113   outertrack->CalculateReferencePoint(angle,point[0]);//local x = global distance to padrowplane
114   if(!outertrack->IsPoint()) return diff;
115   dx[1] = fabs(outertrack->GetPointX()-innertrack->GetFirstPointX());
116   dy[1] = fabs(outertrack->GetPointY()-innertrack->GetFirstPointY());
117   dz[1] = fabs(outertrack->GetPointZ()-innertrack->GetFirstPointZ());
118   
119   diff=0;//This was a tough bug to find....
120   for(Int_t i=0; i<2; i++)
121     diff += sqrt(dx[i]*dx[i] + dy[i]*dy[i] + dz[i]*dz[i]);
122   return diff;
123 }
124
125 void AliHLTTPCGlobalMerger::SlowMerge(const Char_t *path)
126 {
127   //Tuning of parameters. This matches _all_ tracks between two neighbouring
128   //slices, and merges the ones which are closest in space. The difference
129   //is written to a ntuppel, which can be used as a input for the SetParameters
130   //when using normal Merge function.
131   
132   
133   void* ntuple=GetNtuple();
134   AliHLTTPCTrack *track[2];
135   AliHLTTPCTrackArray *tout = GetOutTracks();
136   if(fNSlices<2)
137     {
138       LOG(AliHLTTPCLog::kWarning,"AliHLTTPCGlobalMerger::SlowMerge","Slice Number")
139         <<"Need more than one Slice!"<<ENDLOG;
140       return;
141     }
142   
143   for(Int_t i=0; i<fNSlices; i++)
144     {
145       //if(fNSlices!=18 && i+1 == fNSlices) continue; 
146       Int_t slice = fFirst + i;
147       AliHLTTPCTrackArray *ttt0=GetInTracks(i);
148       Int_t slice2 = i+1;
149       //if(slice2==fNSlices) slice2 =0; 
150       
151       //Make sure slices are on the same side of the TPC
152       if(slice2 == 18) slice2=0;
153       else if(slice2 == 36) slice2=18;
154       AliHLTTPCTrackArray *ttt1=GetInTracks(slice2);
155       //10 degrees -> the border of the slices in local coordinates
156       Float_t angle = AliHLTTPCTransform::Pi()/18; 
157       AliHLTTPCTransform::Local2GlobalAngle(&angle,slice);
158       
159       //In the two following cases, the angle is 2*pi, so set it back to 0 in order for
160       //the calculation of crossing points to be correct.
161       if(slice==17 || slice==35) 
162         angle=0;
163       if(i==0)
164         ttt0->QSort();
165       ttt1->QSort();
166       for(Int_t s0=0;s0<ttt0->GetNTracks();s0++)
167         {
168           AliHLTTPCTrack *track0=ttt0->GetCheckedTrack(s0);
169           if(!track0) continue;
170           track0->CalculateHelix();
171           track0->CalculateEdgePoint(angle);
172           //      if(track0->IsPoint()) AddTrack(tout,track0);
173         }
174       for(Int_t s1=0;s1<ttt1->GetNTracks();s1++)
175         {
176           AliHLTTPCTrack *track1=ttt1->GetCheckedTrack(s1);
177           if(!track1) continue;
178           track1->CalculateHelix();
179           track1->CalculateEdgePoint(angle);
180           //      if(track1->IsPoint())  AddTrack(tout,track1); 
181         }
182       Bool_t merge = kTRUE;
183       while(merge)
184         {
185           Int_t min0=-1,min1=-1;
186           Double_t min=5;
187           Int_t n0=ttt0->GetNTracks(),n1=ttt1->GetNTracks();
188           for(Int_t s0=0;s0<n0;s0++)
189             {
190               AliHLTTPCTrack *track0=ttt0->GetCheckedTrack(s0);
191               if(!track0) continue;
192               if(!track0->IsPoint()) continue;
193               for(Int_t s1=0;s1<n1;s1++)
194                 {
195                   AliHLTTPCTrack *track1=ttt1->GetCheckedTrack(s1);
196                   if(!track1) continue;
197                   if(!track1->IsPoint()) continue;
198                   
199                   //Double_t diff = TrackDiff(track0,track1,angle);
200                   Double_t diff = CheckTracks(track0,track1,slice);
201                   //PrintDiff(track0,track1);
202                   if(diff>=0&&diff<min)
203                     {
204                       min=diff;
205                       min0=s0;
206                       min1=s1;
207                     }
208                 }
209             }
210           if(min0>=0&&min1>=0)
211             {
212               AliHLTTPCTrack *track0=ttt0->GetTrack(min0);
213               AliHLTTPCTrack *track1=ttt1->GetTrack(min1);
214               track[0] = track0;
215               track[1] = track1;
216               SortGlobalTracks(track,2);
217               track1->CalculateEdgePoint((angle+AliHLTTPCTransform::Pi()/9));
218               if(track1->IsPoint())//Check if the track will cross the boundary of yet another slice.
219                 MultiMerge(ttt1,track,2);
220               else
221                 MultiMerge(tout,track,2); 
222               track0->CalculateReferencePoint(angle);
223               track1->CalculateReferencePoint(angle);
224               //PrintDiff(track0,track1);
225               FillNtuple(ntuple,track0,track1);
226               ttt0->Remove(min0);
227               ttt1->Remove(min1);
228               
229             }
230           else merge = kFALSE;
231         }
232       ttt0->Compress();
233       ttt1->Compress();
234       LOG(AliHLTTPCLog::kInformational,"AliHLTTPCGlobalMerger::SlowMerge","Result")
235         <<AliHLTTPCLog::kDec<<"Merged Tracks: "<<tout->GetNTracks()<<" at:"
236         <<angle<<ENDLOG;
237     }
238   Char_t fname[1024];
239   sprintf(fname,"%s/merge_parameters.root",path);
240   WriteNtuple(fname,ntuple);
241 }
242
243 void AliHLTTPCGlobalMerger::Merge()
244 {
245   //Normal merging procedure. Matches tracks which are within limits
246   //set by SetParameters. Parameters can be tuned by SlowMerge.
247   
248   AliHLTTPCTrack *track[2];
249   AliHLTTPCTrackArray *tout = GetOutTracks();
250   if(fNSlices<2)
251     {
252       LOG(AliHLTTPCLog::kWarning,"AliHLTTPCGlobalMerger::Merge","Slice Number")
253         <<"Need more than one Slice!"<<ENDLOG;
254       return;
255     }
256   for(Int_t i=0; i<fNSlices; i++)
257     {
258       //if(fNSlices!=18 && i+1 == fNSlices) continue; 
259       Int_t slice = fFirst + i;
260       AliHLTTPCTrackArray *ttt0=GetInTracks(i);
261       Int_t slice2 = i+1;
262       //if(slice2==fNSlices) slice2 =0;
263       
264       //Make sure slices are on the same side of the TPC
265       if(slice2 == 18) slice2=0;
266       else if(slice2 == 36) slice2=18;
267       AliHLTTPCTrackArray *ttt1=GetInTracks(slice2);
268       //10 degrees -> the border of the slices in local coordinates
269       Float_t angle = AliHLTTPCTransform::Pi()/18; 
270       AliHLTTPCTransform::Local2GlobalAngle(&angle,slice);
271       
272       //In the two following cases, the angle is 2*pi, so set it back to 0 in order for
273       //the calculation of crossing points to be correct.
274       if(slice==17 || slice==35)
275         angle=0;
276       if(i==0)
277         ttt0->QSort();
278       ttt1->QSort();
279       Bool_t *ismatched0  = new Bool_t[ttt0->GetNTracks()];
280       Bool_t *ismatched1  = new Bool_t[ttt1->GetNTracks()];
281       Int_t n0=0,n1=0;
282       for(Int_t s0=0;s0<ttt0->GetNTracks();s0++)
283         {
284           ismatched0[s0]=kFALSE;
285           AliHLTTPCTrack *track0=ttt0->GetCheckedTrack(s0);
286           if(!track0) continue;
287           //track0->CalculateHelix();    //This it done in TrackArray
288           track0->CalculateEdgePoint(angle);
289           if(track0->IsPoint()) 
290             {
291               n0++;
292               track0->CalculateReferencePoint(angle);
293             }
294         }
295       for(Int_t s1=0;s1<ttt1->GetNTracks();s1++)
296         {
297           ismatched1[s1]=kFALSE;
298           AliHLTTPCTrack *track1=ttt1->GetCheckedTrack(s1);
299           if(!track1) continue;
300           //track1->CalculateHelix();   //This is done in TrackArray
301           track1->CalculateEdgePoint(angle);
302           if(track1->IsPoint()) 
303             {
304               n1++;
305               track1->CalculateReferencePoint(angle);
306             }
307         }
308       for(Int_t s0=0;s0<ttt0->GetNTracks();s0++)
309         {
310           if(ismatched0[s0]) continue;
311           AliHLTTPCTrack *track0=ttt0->GetCheckedTrack(s0);
312           if(!track0) continue;
313           if(!track0->IsPoint()) continue;
314           for(Int_t s1=0;s1<ttt1->GetNTracks();s1++)
315             {
316               if(ismatched1[s1]) continue;
317               AliHLTTPCTrack *track1=ttt1->GetCheckedTrack(s1);
318               if(!track1) continue;
319               if(!track1->IsPoint()) continue;
320               if(IsRTrack(track0,track1))
321                 {
322                   track[0] = track0;
323                   track[1] = track1;
324                   SortGlobalTracks(track,2);
325                   Double_t r0 = pow(track[0]->GetLastPointX(),2)+
326                     pow(track[0]->GetLastPointY(),2);
327                   Double_t r1 = pow(track[1]->GetFirstPointX(),2)+
328                     pow(track[1]->GetFirstPointY(),2);
329                   if(r0<r1)
330                     {
331                       MultiMerge(tout,track,2); 
332                       ismatched0[s0]=kTRUE;
333                       ismatched1[s1]=kTRUE;
334                       ttt0->Remove(s0);
335                       ttt1->Remove(s1);
336                       break;
337                       /*
338                         The track is merged, so we will _not_ look for more matches.
339                         Because there could easily be more matches, if a track is being
340                         split within the sector....
341                       */
342                     }
343                 }
344             }
345         }
346       LOG(AliHLTTPCLog::kInformational,"AliHLTTPCGlobalMerger::Merge","Result")
347         <<AliHLTTPCLog::kDec<<"slice0: "<<n0<<" slice1: "<<n1
348         <<" Merged Tracks: "<<tout->GetNTracks()<<ENDLOG;
349       delete [] ismatched0;
350       delete [] ismatched1;
351     }
352 }