During simulation: fill STU region w/ non null time sums
[u/mrichter/AliRoot.git] / TPC / AliCosmicTracker.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 //
18 //
19 //-------------------------------------------------------
20 //          Implementation of the Cosmic tracker
21 //
22 //   Origin:  Xianguo Lu lu@physi.uni-heidelberg.de  Xianguo.Lu@cern.ch
23 //
24 //=========================================================================================
25 // Motivation:
26 // 
27 // In the default reconstruction in the ALICE the cosmic tracks are found  as two independent particles.
28 //
29 // In general any of  subtracks can be used for the physics studies. In order to avoid the double counting,
30 // the track from the upper hemisphere can be used. 
31 //
32 // The momentum resolution is determined by the lever arm (1/L^2) and by the number of clusters
33 // used for the track fitting (1/sqrt(Ncl)). 
34 // Combining/refitting  the two segments together significantly better momentum resolution can be obtained.
35 //    sigma_{1/pt} ~ 8x10^-3  - defaul tracking (e.g only upper track)
36 //    sigma_{1/pt} ~ 8x10^-4  - combined tracking
37 //===========================================================================================
38 // 
39 // Interface/Implementation:
40 // The class AliCosmicTracker provides functionality to find and refit the cosmic tracks. As a starting point, the events reconstruccted using standard tracking are used. 
41 // Input:  AliESDEvent
42 // Output: TClonesArray of the AliESDCosmicTrack
43 //         The array is stored as a data member of the tracker. 
44 //
45 // The cosmic tracker can be called in the user analysis code (standard analisys train, using the AliAnalysisTask, 
46 //    see e.g. AliAnalysisTaskCosmicTracker.h). In oreder to make an analysis simpler and faster it is planned to use the tracker already in the standard reconstruction (To be done).
47 //   
48 //===========================================================================================
49 // Algorithm:
50 // 1. Reads an ESD event      -  SetESDEvent() function
51 // 2. Loop over single tracks -  Process() function 
52 //    cuts are applied for individual ESD tracks (see function ESDtrckCut()). Only ESD tracks with TPCrefit, 
53 //    no kink and with ESDfriends will be selected. 
54 //    User defined cuts (as a pointer to the static function) can be used. (Expert usage)  
55
56 // 3. Double loop over tracks -  Process() function
57 //      a.) if not pair ( see function IsPair() ) continue; 
58 //         To accept the pair the tracks should be close together in the parameter space (AliExternalTrackParam - fP[0]-[4], also cut on ESD-phi and -theta)
59 //         Absolute, and relative (pull) cut are used
60 //         The cuts can be modified beyond default values via SetCut***().
61 //                  
62 //      b.) Each pair is fit via AliTPCCosmicTrackfit::CombineESDtracks
63 //      c.) For each pair one AliESDCosmicTrack object is stored in the fTrackStack, which can be passed out via TClonesArray *arr = fCosmicTracker->GetTrackStack();
64 //
65 //
66 //===========================================================================================
67 // Algorithm numerical debugging:
68 // The AliCosmicTracker can be used in the different debug/verbose level (see fDebugLevel) 
69 // Several intermediate variables can be stored in the trees, printout, or draw.
70 // Given functionality (dumping of variables to the tree) was also used for the tuning of the pair
71 // selection criterias, and for validation of the fit functionality. 
72 //
73 //===========================================================================================
74 // Usage:
75 // AliCosmicTracker *fCosmicTracker = new AliCosmicTracker(debuglevel, tag);
76 // fCosmicTracker->SetESDEvent(fESDEvent); //fTrackStack will be automatically cleared/emptied
77 // Int_t npair = fCosmicTracker->Process(processtag, kprint); //processtag only relavant if (debuglevel & 4) to draw the tracks in png; number of cosmic candidates are returned; if kprint the event is draw to png
78 //
79 //
80 // Advanced usage:
81 // fUserCut can be assigned externally so that additional ESDtrack cut can be applied in the very beginning together with those in ESDtrackCut()
82 //
83 // Example:
84 /*
85 //define static (important!!)  cut function in analysis task, e.g. AliAnalysisTaskCosmicTracker
86 //1) in AliAnalysisTaskCosmicTracker.h
87 static Bool_t TrackCut(AliESDtrack *trk);
88
89 //2) in AliAnalysisTaskCosmicTracker.cxx
90 Bool_t AliAnalysisTaskCosmicTracker::TrackCut(AliESDtrack *trk)
91 {
92   //
93   //external track cut in addition to the one in AliCosmicTracker (example)
94   //
95   if(!trk->GetTRDncls())
96     return kFALSE;
97
98   return kTRUE;
99 }
100 //set user cut function
101 fCosmicTracker = new AliCosmicTracker;
102 fCosmicTracker->SetUserESDtrackCut(AliAnalysisTaskCosmicTracker::TrackCut);
103 */ 
104
105 #include <TTreeStream.h>
106
107 #include "AliESDEvent.h"
108 #include "AliTPCseed.h"
109 #include "AliTrackerBase.h"
110
111 #include "AliESDCosmicTrack.h"
112 #include "AliCosmicTracker.h"
113 #include "AliTPCCosmicUtils.h"
114 #include "AliTPCCosmicTrackfit.h"
115
116 AliCosmicTracker::AliCosmicTracker(const Int_t dlev, const TString tag): 
117   fUserCut(0x0)
118   , fStreamer(0x0), fDebugLevel(dlev)
119   , fESDEvent(0x0)
120   , fCombinedTrackfit(0x0)
121   , fTrackStack(0x0)
122   , fCombRefPt(-999)
123   , fChi2PerCluster(-999)
124   , fImpactD(-999)
125   , fImpactZ(-999)
126   , fIsReuse(kFALSE)
127   , fFindableRatio(-999)
128   , fdPhi(-999)
129   , fdTheta(-999)
130 {
131   //
132   //constructor
133   //
134
135   if(fDebugLevel & 1)
136     fStreamer = new TTreeSRedirector(Form("CosmicTracker_%s.root", tag.Data()));
137
138   fCombinedTrackfit = new AliTPCCosmicTrackfit(0, "AliCosmicTracker");
139   fTrackStack = new TClonesArray("AliESDCosmicTrack",100);
140  
141   for(Int_t ii=0; ii<5; ii++){
142     fDelta[ii] = -999;
143     fPull[ii] = -999;
144   }
145
146   fCutdPhi   = 19e-3*5;
147   fCutdTheta = 10e-3*5;
148   
149   fCutPull[0] = 1.9 *10;
150   fCutPull[1] = 1.5 *1e10;
151   fCutPull[2] = 1.9 *10;//bug-fixed!
152   fCutPull[3] = 0.4 *1e10;
153   fCutPull[4] = 3.6 *10;
154
155   fCutDelta[0] = 0.8   * 10;
156   fCutDelta[1] = 2.7   * 10;
157   fCutDelta[2] = 0.012 * 10;//bug-fixed!
158   fCutDelta[3] = 0.007 * 10;
159   fCutDelta[4] = 0.05  * 10;
160 }
161
162 AliCosmicTracker::~AliCosmicTracker()
163 {
164   //
165   //destructor
166   //
167   delete fStreamer;
168   delete fCombinedTrackfit;
169   delete fTrackStack;
170 }
171
172 void AliCosmicTracker::SetESDEvent(AliESDEvent *esd)
173 {
174   //
175   //set esd event
176   //
177
178   fESDEvent = esd;
179   fTrackStack->Clear();
180 }
181
182 Int_t AliCosmicTracker::Process(const TString tag, const Bool_t kprint)
183 {
184   //
185   //double loop over combinations of esd tracks, cosmic event candidates sotred in fTrackStack
186   //
187
188   Int_t npair=0;
189   const Int_t ntrk = fESDEvent->GetNumberOfTracks();
190   Int_t trkcounter[ntrk];
191   for(Int_t ii=0; ii<ntrk; ii++){
192     trkcounter[ii]=0;
193   }
194
195   Double_t findabler0 = -999;
196   Double_t findabler1 = -999;
197
198   for(Int_t itrk=0; itrk<ntrk; itrk++){
199     if(!ESDtrackCut(fESDEvent->GetTrack(itrk), findabler0))
200       continue;
201
202     for(Int_t jtrk=itrk+1; jtrk<ntrk; jtrk++){
203       if(!ESDtrackCut(fESDEvent->GetTrack(jtrk), findabler1))
204         continue;
205
206       AliESDtrack * trk0 = fESDEvent->GetTrack(itrk);
207       AliESDtrack * trk1 = fESDEvent->GetTrack(jtrk);
208       if( IsPair(trk0, trk1) ){
209         if( fCombinedTrackfit->CombineESDtracks(trk0, trk1)  ){
210           fCombTrackUp = *(fCombinedTrackfit->GetTrackParamUp());
211           fCombTrackLow = *(fCombinedTrackfit->GetTrackParamLow());
212
213           fCombRefPt = fCombinedTrackfit->GetTrackParamLow()->Pt();
214           fChi2PerCluster = fCombinedTrackfit->GetChi2PerCluster();
215
216           fImpactD = fCombinedTrackfit->GetImpactD();
217           fImpactZ = fCombinedTrackfit->GetImpactZ();
218
219           fFindableRatio = TMath::Min(findabler0, findabler1);
220
221           trkcounter[itrk]++;
222           trkcounter[jtrk]++;
223
224           fIsReuse = (trkcounter[itrk]>1 || trkcounter[jtrk]>1);
225           if(fDebugLevel & 1)
226             WriteStreamer(ntrk);
227
228           Int_t idup = itrk;
229           Int_t idlow = jtrk;
230           if(fCombinedTrackfit->IsSwap()){
231             const Int_t idtmp = idup;
232             idup = idlow;
233             idlow = idtmp;
234
235             AliExternalTrackParam tptmp = fTrack0;
236             fTrack0 = fTrack1;
237             fTrack1 = tptmp;
238           }
239
240           if(
241              (fDebugLevel & 4) && 
242              ( (fIsReuse && ntrk<=4) || kprint )
243              ){
244             AliESDtrack * trks[]={fESDEvent->GetTrack(idup), fESDEvent->GetTrack(idlow)};
245             AliTPCCosmicUtils::DrawTracks(trks, Form("reuse_%03d_%03d_%03d_%s", ntrk, itrk, jtrk, tag.Data()));
246           }
247           if(
248              (fDebugLevel & 8) &&
249              fImpactD > 200
250              ){
251             AliESDtrack * trks[]={fESDEvent->GetTrack(idup), fESDEvent->GetTrack(idlow)};
252             AliTPCCosmicUtils::DrawTracks(trks, Form("largevtd_%.f_%03d_%03d_%03d_%s", fImpactD, ntrk, itrk, jtrk, tag.Data()));
253           }
254
255           AliESDCosmicTrack costrk(idup, idlow, fCombinedTrackfit->GetTrackParamUp(), fCombinedTrackfit->GetTrackParamLow(), &fTrack0, &fTrack1, fCombinedTrackfit->GetFitNcls(), fCombinedTrackfit->GetFitLeverArm(), fCombinedTrackfit->GetChi2PerCluster(), fCombinedTrackfit->GetImpactD(), fCombinedTrackfit->GetImpactZ(), fIsReuse, fFindableRatio);
256           new((*fTrackStack)[npair]) AliESDCosmicTrack(costrk);
257           npair++;
258         }
259         else{
260           if(fDebugLevel & 16){
261             if(ntrk==2){
262               AliESDtrack * trks[]={trk0, trk1};
263               AliTPCCosmicUtils::DrawTracks(trks, Form("failCombinedFit_%02d_%03d_%03d_%03d_%s", fCombinedTrackfit->GetStatus(), ntrk, itrk, jtrk, tag.Data()));
264             }
265           }
266         }
267       }
268       else{
269         if(fDebugLevel & 32){
270           if(ntrk==2){
271             AliESDtrack * trks[]={trk0, trk1};
272             AliTPCCosmicUtils::DrawTracks(trks, Form("failIsPair_%02d_%03d_%03d_%03d_%s", fErrFlagIsPair, ntrk, itrk, jtrk, tag.Data()));
273           }
274         }
275       }
276     }
277   }
278   
279   return npair;
280 }
281
282 Bool_t AliCosmicTracker::IsPair(AliESDtrack * trk0, AliESDtrack * trk1)
283 {
284   //
285   //check whether the two tracks come from one cosmic ray
286   //
287
288   fErrFlagIsPair = 0;
289
290   //dphi + pi should = 0
291   fdPhi   = AliTPCCosmicUtils::AngleInRange(trk0->Phi()   - trk1->Phi()   + TMath::Pi());
292   if( TMath::Abs(fdPhi) > fCutdPhi ){
293     fErrFlagIsPair = 1;
294     return kFALSE;
295   }
296
297   fdTheta = AliTPCCosmicUtils::AngleInRange(trk0->Theta() + trk1->Theta() + TMath::Pi());
298   if( TMath::Abs(fdTheta) > fCutdTheta ){
299     fErrFlagIsPair = 2;
300     return kFALSE;
301   }
302
303   //use fIp, the comments on the web is wrong (M. Ivanov)
304   if(!trk0->GetInnerParam()){
305     fErrFlagIsPair = 3;
306     return kFALSE;
307   }
308   if(!trk1->GetInnerParam()){
309     fErrFlagIsPair = 4;
310     return kFALSE;
311   }
312
313   AliExternalTrackParam tmptrk[]={*(trk0->GetInnerParam()), *(trk1->GetInnerParam())};
314
315   if(fDebugLevel & 2){
316     printf("\n************************ raw ESD:\n");
317     AliTPCCosmicUtils::PrintTrackParam(100, &(tmptrk[0]));
318     AliTPCCosmicUtils::PrintTrackParam(101, &(tmptrk[1]));
319     tmptrk[0].Print();
320     tmptrk[1].Print();
321   }
322
323   Double_t xyz0[3], xyz1[3];
324   tmptrk[0].GetXYZ(xyz0);
325   tmptrk[1].GetXYZ(xyz1);
326   const TVector3 gpos0(xyz0), gpos1(xyz1);
327
328   //============================== rotate to common angle (M. Ivanov), since it is not possible to rotate from alpha1 to alpha0 via AliExternalTrackParam::Rotate
329   const Double_t meanalpha = (gpos0-gpos1).Phi();
330   const Double_t alpha0 = tmptrk[0].GetAlpha();
331
332   //track0 closer to mean alpha
333   if( TMath::Abs(AliTPCCosmicUtils::AngleInRange(meanalpha-alpha0)) <TMath::PiOver2() ){
334     if( !AliTPCCosmicUtils::RotateSafe(&(tmptrk[0]), meanalpha) || 
335         !AliTPCCosmicUtils::RotateSafe(&(tmptrk[1]), meanalpha+TMath::Pi()) ){
336       fErrFlagIsPair = 5;
337       return kFALSE;
338     }
339   }
340   //track1 closer to mean alpha
341   else{
342     if( !AliTPCCosmicUtils::RotateSafe(&(tmptrk[1]), meanalpha) || 
343         !AliTPCCosmicUtils::RotateSafe(&(tmptrk[0]), meanalpha+TMath::Pi()) ){
344       fErrFlagIsPair = 6;
345       return kFALSE;
346     }
347   }
348
349   if(fDebugLevel & 2){
350     printf("\n************************ after rotation!!\n");
351     AliTPCCosmicUtils::PrintTrackParam(300, &(tmptrk[0]));
352     AliTPCCosmicUtils::PrintTrackParam(301, &(tmptrk[1]));
353     tmptrk[0].Print();
354     tmptrk[1].Print();
355   }
356
357   //============================== propagate from TPC inner wall to x=0 with correct dedx
358   const Double_t xTogo = 0.0;
359   const Double_t maxStep = 1;
360   const Bool_t rotateTo = kFALSE;
361   const Double_t maxSnp = 0.8;
362   Double_t eloss[2]={-1, 1};
363   //default [0]=upper [1]=lower
364   //tmptrk[0].phi<0 ==> [0]=lower
365   if(AliTPCCosmicUtils::AngleInRange(tmptrk[0].Phi())<0){
366     eloss[0]= 1;
367     eloss[1]=-1;
368   }
369
370   for(Int_t ii=0; ii<2; ii++){
371     if(!AliTrackerBase::PropagateTrackToBxByBz(&(tmptrk[ii]), xTogo, AliTPCCosmicUtils::fgkMass, maxStep, rotateTo, maxSnp, eloss[ii])){
372       fErrFlagIsPair = 7;
373       return kFALSE;
374     }
375   }
376   
377   if(fDebugLevel & 2){
378     printf("\n************************ after dedx corr:\n");
379     AliTPCCosmicUtils::PrintTrackParam(200, &(tmptrk[0]));
380     AliTPCCosmicUtils::PrintTrackParam(201, &(tmptrk[1]));
381     tmptrk[0].Print();
382     tmptrk[1].Print();
383   }
384
385   //____________________________________________________________________________________
386   //____________________________________________________________________________________
387
388   //ESD tracks after reconstruction all have x=0 and
389   //TMath::Abs(alpha0 - alpha1)~ pi ==> back-to-back with angular resolution
390   //[0]: local Y-coordinate of a track (cm);  9.945702e+01 -9.961257e+01 ==> opposite
391   //[1]: local Z-coordinate of a track (cm); 2.677805e+01 2.711143e+01 ==> same
392   //[2]: local sine of the track momentum azimuthal angle; should be the same!! bug-fixed
393   //[3]: tangent of the track momentum dip angle; 1.563563e-01 -1.542005e-01 ==> opposite
394   //[4]: 1/pt (1/(GeV/c)); -1.774935e-01 1.705480e-01 ==> opposite
395
396   fDelta[0] = tmptrk[0].GetParameter()[0] + tmptrk[1].GetParameter()[0];
397   fDelta[1] = tmptrk[0].GetParameter()[1] - tmptrk[1].GetParameter()[1];
398   fDelta[2] = tmptrk[0].GetParameter()[2] - tmptrk[1].GetParameter()[2];//bug-fixed!! should use "-"
399   fDelta[3] = tmptrk[0].GetParameter()[3] + tmptrk[1].GetParameter()[3];
400   fDelta[4] = tmptrk[0].GetParameter()[4] + tmptrk[1].GetParameter()[4];
401
402   fPull[0] = fDelta[0]/TMath::Sqrt(tmptrk[0].GetCovariance()[0] +tmptrk[1].GetCovariance()[0]);
403   fPull[1] = fDelta[1]/TMath::Sqrt(tmptrk[0].GetCovariance()[2] +tmptrk[1].GetCovariance()[2]);
404   fPull[2] = fDelta[2]/TMath::Sqrt(tmptrk[0].GetCovariance()[5] +tmptrk[1].GetCovariance()[5]);
405   fPull[3] = fDelta[3]/TMath::Sqrt(tmptrk[0].GetCovariance()[9] +tmptrk[1].GetCovariance()[9]);
406   fPull[4] = fDelta[4]/TMath::Sqrt(tmptrk[0].GetCovariance()[14]+tmptrk[1].GetCovariance()[14]);
407
408   if(fDebugLevel & 2){
409     for(Int_t ii=0; ii<5; ii++){
410       printf("test %d %e %e -- %e\n", ii, tmptrk[0].GetParameter()[ii], tmptrk[1].GetParameter()[ii], fPull[ii]);
411     }
412   }
413
414   for(Int_t ii=0; ii<5; ii++){
415     if( TMath::Abs(fPull[ii])  > fCutPull[ii] ){
416       fErrFlagIsPair = 10+ii;
417       return kFALSE;
418     }
419     if( TMath::Abs(fDelta[ii]) > fCutDelta[ii] ){
420       fErrFlagIsPair = 20+ii;
421       return kFALSE;
422     }
423   }
424
425   fTrack0 = tmptrk[0];
426   fTrack1 = tmptrk[1];
427
428   return kTRUE;
429 }
430
431 Bool_t AliCosmicTracker::ESDtrackCut(AliESDtrack * trk, Double_t &findabler) 
432 {
433   //
434   //cut on track quality (kink, TPCrefit, findable ratio) and require TPC seed
435   //
436
437   if(fUserCut){
438     if(!fUserCut(trk))
439       return kFALSE;
440   }
441
442   //reject kink
443   if(trk->GetKinkIndex(0)>0)
444     return kFALSE;
445
446   //require refit
447   if(!trk->IsOn(AliESDtrack::kTPCrefit))
448     return kFALSE;
449
450   // due to drift velocity calibration, a track crossing Z=0 may be reconstructed as 2 ESD tracks, so two pairs are formed, each with one part of this track. Solution: cut on findable ratio (require > 0.5) to remove split tracks due to drift velocity calibration systematics on different sides
451   //there is some remaining with isreuse = true, user should cut on findable ratio according to the fraction of isreuse
452   findabler = -999;
453   if(!trk->GetTPCNclsF())
454     return kFALSE;
455
456   findabler = (Double_t)trk->GetTPCNcls()/(Double_t) trk->GetTPCNclsF();
457
458   if(findabler < fgkCutFindable ){
459     return kFALSE;
460   }
461
462   //require ESDfriends
463   if(!AliTPCCosmicUtils::GetTPCseed(trk))
464     return kFALSE;
465
466   return kTRUE;
467 }
468
469 void AliCosmicTracker::WriteStreamer(Int_t ntrk)
470 {
471   //
472   //output to streamer
473   //
474
475   (*fStreamer)<<"CosmicTracker_Streamer"<<
476     "ntrk="<<ntrk<<
477     "combup="<<&fCombTrackUp<<
478     "comblow="<<&fCombTrackLow<<
479     "track0="<<&fTrack0<<
480     "track1="<<&fTrack1<<
481     "vtxd="<<fImpactD<<
482     "vtxz="<<fImpactZ<<
483     "isreuse="<<fIsReuse<<
484     "findabler="<<fFindableRatio<<
485     "pt="<<fCombRefPt<<
486     "chi2="<<fChi2PerCluster<<
487     "dphi="<<fdPhi<<
488     "dtheta="<<fdTheta<<
489     "pull0="<<fPull[0]<<
490     "pull1="<<fPull[1]<<
491     "pull2="<<fPull[2]<<
492     "pull3="<<fPull[3]<<
493     "pull4="<<fPull[4]<<
494     "delta0="<<fDelta[0]<<
495     "delta1="<<fDelta[1]<<
496     "delta2="<<fDelta[2]<<
497     "delta3="<<fDelta[3]<<
498     "delta4="<<fDelta[4]<<
499     "\n";
500 }