TPC module
[u/mrichter/AliRoot.git] / TPC / TPCrec / AliTPCCosmicTrackfit.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 //-------------------------------------------------------
17 //          Implementation of the TPC cosmic trackfit
18 //
19 //   Origin:  Xianguo Lu lu@physi.uni-heidelberg.de  Xianguo.Lu@cern.ch
20 //
21 //=========================================================================================
22 // Motivation:
23 // 
24 // In the default reconstruction in the ALICE the cosmic tracks are found  as two independent particles.
25 //
26 // In general any of  subtracks can be used for the physics studies. In order to avoid the double counting,
27 // the track from the upper hemisphere can be used. 
28 //
29 // The momentum resolution is determined by the lever arm (1/L^2) and by the number of clusters
30 // used for the track fitting (1/sqrt(Ncl)). 
31 // Combining/refitting  the two segments together significantly better momentum resolution can be obtained.
32 //    sigma_{1/pt} ~ 8x10^-3  - defaul tracking (e.g only upper track)
33 //    sigma_{1/pt} ~ 8x10^-4  - combined tracking
34 //===========================================================================================
35 // 
36 // Interface/Implementation:
37 // This class provide the functionality to combine two ESD tracks and perform the trackfit using all related track points.
38 // Input: AliESDtrack's (or AliTPCseed's, expert use)
39 // Output: the trackparams from the combined fit are stored as data member of this class, and they can be accessed via getter or copy.
40 //
41 //===========================================================================================      
42 //
43 // Algorithm:
44 // 1. Read in two ESDtracks (see function CombineESDtracks())
45 //      TPCseeds are got from ESDfriends
46 //      Cuts are applied on the following quantities: #TPC cluster, relation between global-y coordinate of the two tracks, lever arm.
47 // 2. AliTPCCosmicUtils::CombinedFit is called and refit is performed using the TPCseeds from both tracks
48 // 3. The fitting quality is analyzed and stored and can be accessed via function GetStatus(). The detailed judgement of the quality is recorded as "enum CombineStatus"
49 //
50 //
51 //===========================================================================================
52 // Algorithm numerical debugging:
53 // The class can be used in the different debug/verbose level (see fDebugLevel) 
54 // Several intermediate variables can be stored in the trees, printout, or draw.
55 // Given functionality (dumping of variables to the tree) was also used for the tuning of the pair
56 // selection criterias, and for validation of the fit functionality. 
57 //
58 // ----- Debug output:
59 // for (debuglevel & 1)==1 && good fit, the following info saved:
60 /*
61   (*fStreamer)<<"TrackProp"<<
62   "Tup.="<<fTrackparUp<<          //AliExternalTrackParam at uppermost cluster obtained by upward propagation
63   "Tlow.="<<fTrackparLow<<        //AliExternalTrackParam at lowermost cluster obtained by downward propagation
64   "icup.="<<&fInnerClusterUp<<    //TVector3 position of the innermost cluster of the upper track
65   "iclow.="<<&fInnerClusterLow<<
66   "leverarm="<<fLeverArm<<
67   "ncl="<<fFitNcls<<              //number of clusters used in successful propagation  
68   "nmiss="<<fMissNcls<<           //number of clusters failed in propagation, should always be 0 in this case.
69   "chi2="<<fPreChi2<<             //chi2/nfit  
70   "momup="<<  momup <<            //momentum at uppermost cluster with upward propagation
71   "momlow="<< momlow <<           //momentum at lowermost cluster with downward propagation
72   "ptup="<<   ptup <<
73   "ptlow="<<  ptlow <<
74       "\n";
75 */
76 // for (debuglevel & 2)==1, debug info in AliTPCCosmicUtils::FitKernel saved
77 //
78 //===========================================================================================
79 // Usage
80 /*
81   fCosmicTrackfit = new AliTPCCosmicTrackfit(debuglevel, "anystring");
82
83   //order not important; will be internally ordered (potinters modified due to &) such that track0 is the upper one
84   //kfit = kTRUE: good fit, kFALSE: bad fit
85   const Bool_t kfit = fCosmicTrackfit->CombineESDtracks(esdtrack0, esdtrack1);
86
87   //status = 0 for good fit (i.e. kfit=kTRUE), non-0 for bad fit (i.e. kfit=kFALSE), see "enum CombineStatus" definition in header file
88   const Int_t status = fCosmicTrackfit->GetStatus(); 
89 */
90 //===========================================================================================
91 // Efficiency
92 //
93 // for 2011 Feb. cosmic data nch=2 events, the kfit and status look like:
94 /*
95 kfit,status (  0,   1):   16337 /  443901 =   3.680%          //kFailGetTPCseeds
96 kfit,status (  0,   2):    3692 /  443901 =   0.832%          //not both tracks have ncl > AliTPCCosmicUtils::NclsMin()
97 kfit,status (  0,   3):    6527 /  443901 =   1.470%          //clusters in two tracks should be clearly separated in y, i.e. lowest cluster of upper track higher than highest cluster of lower track; otherwise fail 
98 kfit,status (  0,   4):    7033 /  443901 =   1.584%          //fLeverArm<CutLeverArm()
99 kfit,status (  0,   6):    4398 /  443901 =   0.991%          //fail in propagation of at least one cluster
100 kfit,status (  0,   7):     508 /  443901 =   0.114%          //chi2/nfit > MaxChi2()
101 kfit,status (  0,   8):      52 /  443901 =   0.012%          //fail in geting impact parameters
102 kfit,status (  1,   0):  405354 /  443901 =  91.316%          //i.e. 91% of nch=2 events are successfully fitted.
103 */
104 //
105 //===========================================================================================
106 // Resolution
107 //
108 // for muon momentum small than 20 GeV, energy loss in muon filter is visable when compaing fTrackparUp and fTrackparLow; energy loss estimated as 5 MeV/cm.
109 // particle traversing muon filter can be rejected by requiring "fInnerClusterUp.fZ > -40 && fInnerClusterLow.fZ > -40"
110 // momentum resolution is estimated by comparing the trackfit result by upward propagation through odd pad rows and that by downward propagation through even pad rows. Number of clusters used in this case is only half of that in normal usage.
111 // RMS of log10 p = 0.01 at 10 GeV/c, 0.1 at 100 GeV/c, 0.5 at 1 TeV/c.
112 // muon filter deteriorates momentum resolution by about +0.01 (absolute value).
113 //
114
115 #include <TAxis.h>
116 #include <TCanvas.h>
117 #include <TFile.h>
118 #include <TGraph.h>
119 #include <TTreeStream.h>
120 #include <TVector3.h>
121
122 #include "AliESDtrack.h"
123 #include "AliESDfriendTrack.h"
124 #include "AliTPCseed.h"
125 #include "AliTrackerBase.h"
126 #include "AliTrackPointArray.h"
127
128 #include "AliTPCCosmicUtils.h"
129 #include "AliTPCCosmicTrackfit.h"
130
131 AliTPCCosmicTrackfit::AliTPCCosmicTrackfit(const Int_t dlev, const TString tag):
132   fStreamer(0x0), fDebugLevel(dlev)
133   , fSeedUp(0x0), fSeedLow(0x0), fTrackparUp(0x0), fTrackparLow(0x0), fIniTrackparUp(0x0), fIniTrackparLow(0x0)
134   , fStatus(-999)
135   , fKswap(0)
136   , fInnerClusterUp(-999,-999,-999)
137   , fInnerClusterLow(-999,-999,-999)
138   , fLeverArm(-999)
139   , fFitNcls(-999), fMissNcls(-999), fPreChi2(-999), fFitLeverArm(-999), fImpactD(-999), fImpactZ(-999)
140   , fRowStartShift(-999)
141   , fRowStep(-999)
142   , fXMin(-999)
143   , fXMax(-999)
144 {
145   //
146   //Constructor
147   //
148
149   if(fDebugLevel>0)
150     fStreamer = new TTreeSRedirector(Form("CosmicTrackfit_%s.root", tag.Data()));
151
152   SetRow(0,1);
153   SetX(-1e10, 1e10);
154 }
155
156 AliTPCCosmicTrackfit::~AliTPCCosmicTrackfit()
157 {
158   //
159   //Destructor
160   //
161
162   delete fStreamer;
163   
164   delete fTrackparUp;
165   delete fTrackparLow;
166   delete fIniTrackparUp;
167   delete fIniTrackparLow;
168 }
169
170 Bool_t AliTPCCosmicTrackfit::CombineESDtracks(AliESDtrack * &trk0, AliESDtrack *&trk1)
171 {
172   //
173   //Get TPCseeds from the 2 ESDtracks, swap TPCseeds and ESDTracks (if necessary) according to y (0:upper 1:lower), perform trackfit using TPCseeds
174   //if fStatus==0, i.e. combine is successful, swap of the ESDtracks is kept since pointer *& is used
175   //
176
177   IniCombineESDtracks();
178
179   if(!GetTPCseeds(trk0, trk1)){
180     return kFALSE; 
181   }
182
183   CombineTPCseeds();
184
185   if(fStatus == 0){
186     if(fKswap){
187       AliESDtrack * tmptrk = trk0;
188       trk0 = trk1;
189       trk1 = tmptrk;
190     }
191     return kTRUE;
192   }
193   else
194     return kFALSE;
195 }
196
197 Bool_t AliTPCCosmicTrackfit::CombineTPCseedsFast(AliTPCseed * tpcseeds[], const AliExternalTrackParam * trkpars[])
198 {
199   //
200   //do combined track fit for given TPC seeds and initial trackpar, [0]: upper, [1]: lower
201   //
202
203   IniCombineESDtracks();
204
205   fSeedUp = tpcseeds[0];
206   fSeedLow = tpcseeds[1];
207
208   fTrackparUp =  new AliExternalTrackParam(*(trkpars[0]));
209   fTrackparLow = new AliExternalTrackParam(*(trkpars[1]));
210
211   AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
212   const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
213   TTreeSRedirector * debugstreamer = 0x0;
214   if(fDebugLevel&2){
215     debugstreamer = fStreamer;
216   }
217
218   AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
219
220   Update();
221
222   if(fStatus == 0){
223     return kTRUE;
224   }
225   else
226     return kFALSE;
227 }
228
229 Bool_t AliTPCCosmicTrackfit::CombineTPCseeds(AliTPCseed * &seed0, AliTPCseed *&seed1)
230 {
231   //
232   //same as AliTPCCosmicTrackfit::CombineESDtracks, except that the seeds are passed in from outside, which can be still unordered
233   //if fStatus==0, i.e. combine is successful, swap of the TPCseeds is kept since pointer *& is used
234   //
235   IniCombineESDtracks();
236
237   fSeedUp  = seed0;
238   fSeedLow = seed1;
239   
240   CombineTPCseeds();
241
242   if(fStatus==0){
243     if(fKswap){
244       AliTPCseed * tmpseed = seed0;
245       seed0 = seed1;
246       seed1 = tmpseed;
247     }
248     return kTRUE;
249   }
250   else 
251     return kFALSE;
252 }
253
254 void AliTPCCosmicTrackfit::Print() const
255 {
256   //
257   //print out variable values
258   //
259   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);
260 }
261
262 TVector3 AliTPCCosmicTrackfit::ImpactParameter2D() const
263 {
264   //
265   //calculate the 2D-impactparameter from (0,0)
266   //
267
268   const TVector3 p0(0,0,0);
269   const TVector3 l1(fInnerClusterUp.X(),  fInnerClusterUp.Y(),  0);
270   const TVector3 l2(fInnerClusterLow.X(), fInnerClusterLow.Y(), 0);
271
272   TVector3 vtx;
273   AliTPCCosmicUtils::Point2LineDist(p0, l1, l2, &vtx);
274   return vtx;
275 }
276
277 TVector3 AliTPCCosmicTrackfit::ImpactParameter3D() const
278 {
279   //
280   //calculate the 3D-impactparameter from (0,0,0)
281   //
282
283   const TVector3 p0(0,0,0);
284   const TVector3 l1(fInnerClusterUp);
285   const TVector3 l2(fInnerClusterLow);
286
287   TVector3 vtx;
288   AliTPCCosmicUtils::Point2LineDist(p0, l1, l2, &vtx);
289
290   //==========================
291   /*
292   printf("testing... \n");
293   TVector3 tmpv = ImpactParameter2D();
294   if(fabs(tmpv.Mag()-vtx.Pt())>1e-6){
295     printf("strange!!! %e %e %e\n", tmpv.Mag(), vtx.Pt(), fabs(tmpv.Mag()-vtx.Pt()));
296     vtx.Print();
297     tmpv.Print();
298     fInnerClusterUp.Print();
299     fInnerClusterLow.Print();
300   }
301   */
302   //==========================
303   return vtx;
304 }
305
306 Double_t AliTPCCosmicTrackfit::MinPhi() const
307 {
308   //
309   //the smaller phi of the two tracks w.r.t. horizon
310   //
311   Double_t fsp[] = {TMath::Abs(TMath::Sin(fTrackparUp->Phi())), TMath::Abs(TMath::Sin(fTrackparLow->Phi()))};;
312   return TMath::ASin(TMath::Min(fsp[0], fsp[1])) * TMath::RadToDeg();
313 }
314 //===================================================================================================
315 //===================================================================================================
316
317 void AliTPCCosmicTrackfit::IniCombineESDtracks()
318 {
319   //
320   //initialization, for reuse of the same AliTPCCosmicTrackfit instance
321   //
322
323   fSeedUp = 0x0;
324   fSeedLow = 0x0;
325   delete fTrackparUp;
326   delete fTrackparLow;
327   fTrackparUp = 0x0;
328   fTrackparLow = 0x0;
329
330   delete fIniTrackparUp;
331   delete fIniTrackparLow;
332   fIniTrackparUp = 0x0;
333   fIniTrackparLow = 0x0;
334
335   fStatus = 0;
336   fKswap = kFALSE;
337 }
338
339 void AliTPCCosmicTrackfit::CombineTPCseeds()
340 {
341   //
342   //do combined trackfit using TPCseeds
343   //
344
345   if(
346      !CheckNcls()
347      || !AnaSeeds() // fSeedUp/Low swapped here! check by runTest.C 
348      || !CheckLeverArm()
349      )
350     return;
351
352   //AliExternalTrackParam object created
353   fTrackparUp  = AliTPCCosmicUtils::MakeSeed(fSeedUp);
354   fTrackparLow = AliTPCCosmicUtils::MakeSeed(fSeedLow);
355   if(!fTrackparUp || !fTrackparLow){
356     fStatus = kFailMakeSeed;
357     return;
358   }
359
360   fIniTrackparUp = new AliExternalTrackParam(*fTrackparUp);
361   fIniTrackparLow = new AliExternalTrackParam(*fTrackparLow);
362
363   AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
364   const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
365   TTreeSRedirector * debugstreamer = 0x0;
366   if(fDebugLevel&2){
367     debugstreamer = fStreamer;
368   }
369
370   AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
371
372   Update();
373
374   return;
375 }
376
377 void AliTPCCosmicTrackfit::Update()
378 {
379   //
380   //Update variables depending on the fit result
381   //
382
383   if(fMissNcls || fFitNcls==0){
384     fStatus = kFailPropagation;
385     return;
386   }
387
388   fPreChi2 /= fFitNcls;
389   if(fPreChi2>MaxChi2()){
390     fStatus = kFailChi2;
391     return;
392   }
393
394   if(fImpactD<0){
395     fStatus = kFailImpact;
396     return;
397   }
398
399   if( fStatus == 0 && (fDebugLevel&1) ){
400     Double_t momup  = fTrackparUp->P();
401     Double_t momlow = fTrackparLow->P();
402     Double_t ptup   = fTrackparUp->Pt();
403     Double_t ptlow  = fTrackparLow->Pt();
404
405     (*fStreamer)<<"TrackProp"<<
406       "Tup.="<<fTrackparUp<<
407       "Tlow.="<<fTrackparLow<<
408       "icup.="<<&fInnerClusterUp<<
409       "iclow.="<<&fInnerClusterLow<<
410       "leverarm="<<fLeverArm<<
411       "ncl="<<fFitNcls<<
412       "nmiss="<<fMissNcls<<
413       "chi2="<<fPreChi2<<
414       "momup="<<  momup <<
415       "momlow="<< momlow <<
416       "ptup="<<   ptup <<
417       "ptlow="<<  ptlow <<
418       "\n";
419   }
420 }
421
422 Bool_t AliTPCCosmicTrackfit::CheckLeverArm()
423 {
424   //
425   //if lever arm is too short, no need to use combined track fit. 
426   //On the other hand, short lever arm from two tracks mostly means they are fake pairs.
427   //lever arm extents over one quadrant, e.g. (0,250)-(250,0): 250*sqrt(2)~350
428   //
429   if(fLeverArm<CutLeverArm()){
430     fStatus = kFailLeverArm;
431     return kFALSE;
432   }
433   else 
434     return kTRUE;
435 }
436
437 Bool_t AliTPCCosmicTrackfit::AnaSeeds()
438 {
439   //
440   //swap seeds (if necessary) so that (y of fSeedUp) > (y of fSeedLow)
441   //
442
443   //---------------------------------- navigate through all clusters ----------------------------------
444   AliTPCseed ** seeds[]={&fSeedUp, &fSeedLow};
445
446   //min, max according to y
447   TVector3 singlemin[2], singlemax[2];
448   for(Int_t ii=0; ii<2; ii++){
449     singlemin[ii].SetXYZ( 1e10,  1e10,  1e10);
450     singlemax[ii].SetXYZ(-1e10, -1e10, -1e10);
451   }
452   
453   for(Int_t itrk=0; itrk<2; itrk++){
454     for(Int_t irow=0; irow<AliTPCCosmicUtils::NRow(); irow++){
455       const AliTPCclusterMI * cls = (*seeds[itrk])->GetClusterPointer(irow);
456       if(!cls)
457         continue;
458       
459       Float_t xyz[3]={-999,-999,-999};
460       cls->GetGlobalXYZ(xyz);
461       if(xyz[1]<singlemin[itrk].Y()){
462         singlemin[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
463       }
464       if(xyz[1]>singlemax[itrk].Y()){
465         singlemax[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
466       }
467     }
468   }
469
470   //--------------------------------
471
472   //kpass true if y of the two seeds clearly separate: min of one > max of the other
473   Bool_t kpass = kFALSE;
474
475   fInnerClusterUp.SetXYZ(-999,-999,-999);
476   fInnerClusterLow.SetXYZ(-999,-999,-999);
477   TVector3 combinedmin, combinedmax;
478   if(singlemin[0].Y()> singlemax[1].Y()){
479     fInnerClusterUp  = singlemin[0];
480     fInnerClusterLow = singlemax[1];
481
482     //no need to swap
483     fKswap = kFALSE;
484
485     kpass = kTRUE;
486
487     combinedmax = singlemax[0];
488     combinedmin = singlemin[1];
489   }
490   else if(singlemin[1].Y()> singlemax[0].Y()){
491     fInnerClusterUp  = singlemin[1];
492     fInnerClusterLow = singlemax[0];
493   
494     //have to be swapped
495     fKswap = kTRUE;
496     AliTPCseed *tmp=*(seeds[0]);
497     *(seeds[0])=*(seeds[1]);
498     *(seeds[1])=tmp;
499     
500     kpass = kTRUE;
501
502     combinedmax = singlemax[1];
503     combinedmin = singlemin[0];
504   }           
505   else
506     kpass = kFALSE;
507
508   const TVector3 comdelta = combinedmax-combinedmin;
509   fLeverArm = comdelta.Pt();
510
511   if(!kpass){
512     fStatus = kFailSwapSeeds;
513     return kFALSE;
514   }
515   else
516     return kTRUE;
517 }
518
519 Bool_t AliTPCCosmicTrackfit::CheckNcls()
520 {
521   //
522   //check number of clusters in TPCseed, for too small number MakeSeed will fail
523   //
524   if( fSeedUp->GetNumberOfClusters()<AliTPCCosmicUtils::NclsMin() || fSeedLow->GetNumberOfClusters()<AliTPCCosmicUtils::NclsMin() ){
525     fStatus = kFailNclsMin;
526     return kFALSE;
527   }
528   else
529     return kTRUE;
530 }
531
532 Bool_t AliTPCCosmicTrackfit::GetTPCseeds(const AliESDtrack *trk0,  const AliESDtrack *trk1)
533 {
534   //
535   //Get TPC seeds from ESDfriendTrack
536   //
537   fSeedUp  = AliTPCCosmicUtils::GetTPCseed(trk0);
538   fSeedLow = AliTPCCosmicUtils::GetTPCseed(trk1);
539
540   if(!fSeedUp || !fSeedLow){
541     fStatus = kFailGetTPCseeds;
542     return kFALSE;
543   }
544
545   return kTRUE;
546 }
547