ATO-17 - reset cach matrix and current matrix after dumping to the tree (not before)
[u/mrichter/AliRoot.git] / TPC / AliTPCCombinedTrackfit.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 //
16 // Combine cosmic track pairs (upper, lower) and do track fitting
17 // oooooOOOOOooooo
18 // oooooOOOOOooooo
19 // oooooOOOOOooooo
20 //
21 //  Xianguo Lu <lu@physi.uni-heidelberg.de>
22
23 #include <TAxis.h>
24 #include <TCanvas.h>
25 #include <TFile.h>
26 #include <TGraph.h>
27 #include <TTreeStream.h>
28 #include <TVector3.h>
29
30 #include "AliESDtrack.h"
31 #include "AliESDfriendTrack.h"
32 #include "AliTPCseed.h"
33 #include "AliTrackerBase.h"
34 #include "AliTrackPointArray.h"
35
36 #include "AliTPCCosmicUtils.h"
37 #include "AliTPCCombinedTrackfit.h"
38
39 AliTPCCombinedTrackfit::AliTPCCombinedTrackfit(const Int_t dlev, const TString tag):
40   fStreamer(0x0), fDebugLevel(dlev)
41   , fSeedUp(0x0), fSeedLow(0x0), fTrackparUp(0x0), fTrackparLow(0x0)
42   , fStatus(-999)
43   , fLeverArm(-999)
44   , fFitNcls(-999), fMissNcls(-999), fPreChi2(-999)
45 {
46   //
47   //Constructor
48   //
49   fInnerClusterUp.SetXYZ(-999,-999,-999);
50   fInnerClusterLow.SetXYZ(-999,-999,-999);
51
52   if(fDebugLevel>0)
53     fStreamer = new TTreeSRedirector(Form("CombinedTrackfit_%s.root", tag.Data()));
54 }
55
56 AliTPCCombinedTrackfit::~AliTPCCombinedTrackfit()
57 {
58   //
59   //Destructor
60   //
61
62   delete fStreamer;
63   
64   delete fTrackparUp;
65   delete fTrackparLow;
66 }
67
68 Bool_t AliTPCCombinedTrackfit::CombineESDtracks(AliESDtrack * &trk0, AliESDtrack *&trk1)
69 {
70   //
71   //Get TPCseeds from the 2 ESDtracks, swap TPCseeds and ESDTracks (if necessary) according to y (0:upper 1:lower), perform trackfit using TPCseeds
72   //if fStatus==0, i.e. combine is successful, swap of the ESDtracks is kept since pointer *& is used
73   //
74
75   IniCombineESDtracks();
76
77   if(!GetTPCseeds(trk0, trk1)){
78     return kFALSE; 
79   }
80
81   Bool_t kswap = kFALSE;
82   CombineTPCseeds(kswap);
83
84   if(fStatus == 0){
85     if(kswap){
86       AliESDtrack * tmptrk = trk0;
87       trk0 = trk1;
88       trk1 = tmptrk;
89     }
90     return kTRUE;
91   }
92   else
93     return kFALSE;
94 }
95
96 Bool_t AliTPCCombinedTrackfit::CombineTPCseeds(AliTPCseed * &seed0, AliTPCseed *&seed1)
97 {
98   //
99   //same as AliTPCCombinedTrackfit::CombineESDtracks, except that the seeds are passed in from outside, which can be still unordered
100   //if fStatus==0, i.e. combine is successful, swap of the TPCseeds is kept since pointer *& is used
101   //
102   IniCombineESDtracks();
103
104   fSeedUp  = seed0;
105   fSeedLow = seed1;
106   
107   Bool_t kswap = kFALSE;
108   CombineTPCseeds(kswap);
109
110   if(fStatus==0){
111     if(kswap){
112       AliTPCseed * tmpseed = seed0;
113       seed0 = seed1;
114       seed1 = tmpseed;
115     }
116     return kTRUE;
117   }
118   else 
119     return kFALSE;
120 }
121
122 void AliTPCCombinedTrackfit::Print() const
123 {
124   //
125   //print out variable values
126   //
127   printf("Status %2d NclsU %3d NclsD %3d ZinnerU %7.2f ZinnerD %7.2f LeverArm %7.2f\n", fStatus, fSeedUp->GetNumberOfClusters(), fSeedLow->GetNumberOfClusters(), fInnerClusterUp.Z(), fInnerClusterLow.Z(), fLeverArm);
128 }
129
130 Double_t AliTPCCombinedTrackfit::ImpactParameter() const
131 {
132   //
133   //calculate the impactparameter from (0,0,0)
134   //
135   const TVector3 p0(0,0,0);
136   const TVector3 va = p0 - fInnerClusterUp;
137   const TVector3 vb = fInnerClusterLow - fInnerClusterUp;
138
139   const TVector3 dd = va.Cross(vb);
140
141   return dd.Mag()/vb.Mag();
142 }
143
144 Double_t AliTPCCombinedTrackfit::MinPhi() const
145 {
146   //
147   //the smaller phi of the two tracks w.r.t. horizon
148   //
149   Double_t fsp[] = {fabs(sin(fTrackparUp->Phi())), fabs(sin(fTrackparLow->Phi()))};;
150   return asin(TMath::Min(fsp[0], fsp[1])) * TMath::RadToDeg();
151 }
152 //===================================================================================================
153 //===================================================================================================
154
155 void AliTPCCombinedTrackfit::IniCombineESDtracks()
156 {
157   //
158   //initialization, for reuse of the same AliTPCCombinedTrackfit instance
159   //
160
161   fSeedUp = 0x0;
162   fSeedLow = 0x0;
163   delete fTrackparUp;
164   delete fTrackparLow;
165   fTrackparUp = 0x0;
166   fTrackparLow = 0x0;
167
168   fStatus = 0;
169 }
170
171 void AliTPCCombinedTrackfit::CombineTPCseeds(Bool_t &kswap)
172 {
173   //
174   //do combined trackfit using TPCseeds
175   //
176
177   if(
178      !CheckNcls()
179      || !AnaSeeds(kswap) 
180      || !CheckLeverArm()
181      )
182     return;
183
184   //AliExternalTrackParam object created
185   fTrackparUp  = AliTPCCosmicUtils::MakeSeed(fSeedUp);
186   fTrackparLow = AliTPCCosmicUtils::MakeSeed(fSeedLow);
187   if(!fTrackparUp || !fTrackparLow){
188     fStatus = kFailMakeSeed;
189     return;
190   }
191
192   AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
193   const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
194   TTreeSRedirector * debugstreamer = 0x0;
195   if(fDebugLevel&2){
196     debugstreamer = fStreamer;
197   }
198
199   AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fFitNcls, fMissNcls, fPreChi2, debugstreamer);
200
201   Update();
202
203   return;
204 }
205
206 void AliTPCCombinedTrackfit::Update()
207 {
208   //
209   //Update variables depending on the fit result
210   //
211
212   if(fMissNcls || fFitNcls==0){
213     fStatus = kFailPropagation;
214     return;
215   }
216
217   fPreChi2 /= fFitNcls;
218   if(fPreChi2>fgkMaxChi2){
219     fStatus = kFailChi2;
220     return;
221   }
222
223   if( fStatus == 0 && (fDebugLevel&1) ){
224     Double_t momup  = fTrackparUp->P();
225     Double_t momlow = fTrackparLow->P();
226     Double_t ptup   = fTrackparUp->Pt();
227     Double_t ptlow  = fTrackparLow->Pt();
228
229     (*fStreamer)<<"TrackProp"<<
230       "Tup.="<<fTrackparUp<<
231       "Tlow.="<<fTrackparLow<<
232       "icup.="<<&fInnerClusterUp<<
233       "iclow.="<<&fInnerClusterLow<<
234       "leverarm="<<fLeverArm<<
235       "ncl="<<fFitNcls<<
236       "nmiss="<<fMissNcls<<
237       "chi2="<<fPreChi2<<
238       "momup="<<  momup <<
239       "momlow="<< momlow <<
240       "ptup="<<   ptup <<
241       "ptlow="<<  ptlow <<
242       "\n";
243   }
244 }
245
246 Bool_t AliTPCCombinedTrackfit::CheckLeverArm()
247 {
248   //
249   //if lever arm is too short, no need to use combined track fit. 
250   //On the other hand, short lever arm from two tracks mostly means they are fake pairs.
251   //lever arm extents over one quadrant, e.g. (0,250)-(250,0): 250*sqrt(2)~350
252   //
253   if(fLeverArm<fgkCutLeverArm){
254     fStatus = kFailLeverArm;
255     return kFALSE;
256   }
257   else 
258     return kTRUE;
259 }
260
261 Bool_t AliTPCCombinedTrackfit::AnaSeeds(Bool_t &kswap)
262 {
263   //
264   //swap seeds (if necessary) so that (y of fSeedUp) > (y of fSeedLow)
265   //
266
267   //---------------------------------- navigate through all clusters ----------------------------------
268   AliTPCseed ** seeds[]={&fSeedUp, &fSeedLow};
269
270   //min, max according to y
271   TVector3 singlemin[2], singlemax[2];
272   for(Int_t ii=0; ii<2; ii++){
273     singlemin[ii].SetXYZ( 1e10,  1e10,  1e10);
274     singlemax[ii].SetXYZ(-1e10, -1e10, -1e10);
275   }
276   
277   for(Int_t itrk=0; itrk<2; itrk++){
278     for(Int_t irow=0; irow<AliTPCCosmicUtils::fgkNRow; irow++){
279       const AliTPCclusterMI * cls = (*seeds[itrk])->GetClusterPointer(irow);
280       if(!cls)
281         continue;
282       
283       Float_t xyz[3]={-999,-999,-999};
284       cls->GetGlobalXYZ(xyz);
285       if(xyz[1]<singlemin[itrk].Y()){
286         singlemin[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
287       }
288       if(xyz[1]>singlemax[itrk].Y()){
289         singlemax[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
290       }
291     }
292   }
293
294   //--------------------------------
295
296   //kpass true if y of the two seeds clearly separate: min of one > max of the other
297   Bool_t kpass = kFALSE;
298
299   fInnerClusterUp.SetXYZ(-999,-999,-999);
300   fInnerClusterLow.SetXYZ(-999,-999,-999);
301   TVector3 combinedmin, combinedmax;
302   if(singlemin[0].Y()> singlemax[1].Y()){
303     fInnerClusterUp  = singlemin[0];
304     fInnerClusterLow = singlemax[1];
305
306     //no need to swap
307     kswap = kFALSE;
308
309     kpass = kTRUE;
310
311     combinedmax = singlemax[0];
312     combinedmin = singlemin[1];
313   }
314   else if(singlemin[1].Y()> singlemax[0].Y()){
315     fInnerClusterUp  = singlemin[1];
316     fInnerClusterLow = singlemax[0];
317   
318     //have to be swapped
319     kswap = kTRUE;
320     AliTPCseed *tmp=*(seeds[0]);
321     *(seeds[0])=*(seeds[1]);
322     *(seeds[1])=tmp;
323     
324     kpass = kTRUE;
325
326     combinedmax = singlemax[1];
327     combinedmin = singlemin[0];
328   }           
329   else
330     kpass = kFALSE;
331
332   const TVector3 comdelta = combinedmax-combinedmin;
333   fLeverArm = comdelta.Pt();
334
335   if(!kpass){
336     fStatus = kFailSwapSeeds;
337     return kFALSE;
338   }
339   else
340     return kTRUE;
341 }
342
343 Bool_t AliTPCCombinedTrackfit::CheckNcls()
344 {
345   //
346   //check number of clusters in TPCseed, for too small number MakeSeed will fail
347   //
348   if( fSeedUp->GetNumberOfClusters()<AliTPCCosmicUtils::fgkNclsMin || fSeedLow->GetNumberOfClusters()<AliTPCCosmicUtils::fgkNclsMin ){
349     fStatus = kFailNclsMin;
350     return kFALSE;
351   }
352   else
353     return kTRUE;
354 }
355
356 Bool_t AliTPCCombinedTrackfit::GetTPCseeds(const AliESDtrack *trk0,  const AliESDtrack *trk1)
357 {
358   //
359   //Get TPC seeds from ESDfriendTrack
360   //
361   fSeedUp  = AliTPCCosmicUtils::GetTPCseed(trk0);
362   fSeedLow = AliTPCCosmicUtils::GetTPCseed(trk1);
363
364   if(!fSeedUp || !fSeedLow){
365     fStatus = kFailGetTPCseeds;
366     return kFALSE;
367   }
368
369   return kTRUE;
370 }
371