1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 //-------------------------------------------------------
20 // Implementation of the Cosmic tracker
22 // Origin: Xianguo Lu lu@physi.uni-heidelberg.de Xianguo.Lu@cern.ch
24 //=========================================================================================
27 // In the default reconstruction in the ALICE the cosmic tracks are found as two independent particles.
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.
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 //===========================================================================================
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.
42 // Output: TClonesArray of the AliESDCosmicTrack
43 // The array is stored as a data member of the tracker.
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).
48 //===========================================================================================
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)
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***().
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();
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.
73 //===========================================================================================
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
81 // fUserCut can be assigned externally so that additional ESDtrack cut can be applied in the very beginning together with those in ESDtrackCut()
85 //define static (important!!) cut function in analysis task, e.g. AliAnalysisTaskCosmicTracker
86 //1) in AliAnalysisTaskCosmicTracker.h
87 static Bool_t TrackCut(AliESDtrack *trk);
89 //2) in AliAnalysisTaskCosmicTracker.cxx
90 Bool_t AliAnalysisTaskCosmicTracker::TrackCut(AliESDtrack *trk)
93 //external track cut in addition to the one in AliCosmicTracker (example)
95 if(!trk->GetTRDncls())
100 //set user cut function
101 fCosmicTracker = new AliCosmicTracker;
102 fCosmicTracker->SetUserESDtrackCut(AliAnalysisTaskCosmicTracker::TrackCut);
105 #include <TTreeStream.h>
107 #include "AliESDEvent.h"
108 #include "AliTPCseed.h"
109 #include "AliTrackerBase.h"
111 #include "AliESDCosmicTrack.h"
112 #include "AliCosmicTracker.h"
113 #include "AliTPCCosmicUtils.h"
114 #include "AliTPCCosmicTrackfit.h"
116 AliCosmicTracker::AliCosmicTracker(const Int_t dlev, const TString tag):
118 , fStreamer(0x0), fDebugLevel(dlev)
120 , fCombinedTrackfit(0x0)
123 , fChi2PerCluster(-999)
127 , fFindableRatio(-999)
136 fStreamer = new TTreeSRedirector(Form("CosmicTracker_%s.root", tag.Data()));
138 fCombinedTrackfit = new AliTPCCosmicTrackfit(0, "AliCosmicTracker");
139 fTrackStack = new TClonesArray("AliESDCosmicTrack",100);
141 for(Int_t ii=0; ii<5; ii++){
147 fCutdTheta = 10e-3*5;
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;
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;
162 AliCosmicTracker::~AliCosmicTracker()
168 delete fCombinedTrackfit;
172 void AliCosmicTracker::SetESDEvent(AliESDEvent *esd)
179 fTrackStack->Clear();
182 Int_t AliCosmicTracker::Process(const TString tag, const Bool_t kprint)
185 //double loop over combinations of esd tracks, cosmic event candidates sotred in fTrackStack
189 const Int_t ntrk = fESDEvent->GetNumberOfTracks();
190 Int_t trkcounter[ntrk];
191 for(Int_t ii=0; ii<ntrk; ii++){
195 Double_t findabler0 = -999;
196 Double_t findabler1 = -999;
198 for(Int_t itrk=0; itrk<ntrk; itrk++){
199 if(!ESDtrackCut(fESDEvent->GetTrack(itrk), findabler0))
202 for(Int_t jtrk=itrk+1; jtrk<ntrk; jtrk++){
203 if(!ESDtrackCut(fESDEvent->GetTrack(jtrk), findabler1))
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());
213 fCombRefPt = fCombinedTrackfit->GetTrackParamLow()->Pt();
214 fChi2PerCluster = fCombinedTrackfit->GetChi2PerCluster();
216 fImpactD = fCombinedTrackfit->GetImpactD();
217 fImpactZ = fCombinedTrackfit->GetImpactZ();
219 fFindableRatio = TMath::Min(findabler0, findabler1);
224 fIsReuse = (trkcounter[itrk]>1 || trkcounter[jtrk]>1);
230 if(fCombinedTrackfit->IsSwap()){
231 const Int_t idtmp = idup;
235 AliExternalTrackParam tptmp = fTrack0;
242 ( (fIsReuse && ntrk<=4) || kprint )
244 AliESDtrack * trks[]={fESDEvent->GetTrack(idup), fESDEvent->GetTrack(idlow)};
245 AliTPCCosmicUtils::DrawTracks(trks, Form("reuse_%03d_%03d_%03d_%s", ntrk, itrk, jtrk, tag.Data()));
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()));
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);
260 if(fDebugLevel & 16){
262 AliESDtrack * trks[]={trk0, trk1};
263 AliTPCCosmicUtils::DrawTracks(trks, Form("failCombinedFit_%02d_%03d_%03d_%03d_%s", fCombinedTrackfit->GetStatus(), ntrk, itrk, jtrk, tag.Data()));
269 if(fDebugLevel & 32){
271 AliESDtrack * trks[]={trk0, trk1};
272 AliTPCCosmicUtils::DrawTracks(trks, Form("failIsPair_%02d_%03d_%03d_%03d_%s", fErrFlagIsPair, ntrk, itrk, jtrk, tag.Data()));
282 Bool_t AliCosmicTracker::IsPair(AliESDtrack * trk0, AliESDtrack * trk1)
285 //check whether the two tracks come from one cosmic ray
290 //dphi + pi should = 0
291 fdPhi = AliTPCCosmicUtils::AngleInRange(trk0->Phi() - trk1->Phi() + TMath::Pi());
292 if( TMath::Abs(fdPhi) > fCutdPhi ){
297 fdTheta = AliTPCCosmicUtils::AngleInRange(trk0->Theta() + trk1->Theta() + TMath::Pi());
298 if( TMath::Abs(fdTheta) > fCutdTheta ){
303 //use fIp, the comments on the web is wrong (M. Ivanov)
304 if(!trk0->GetInnerParam()){
308 if(!trk1->GetInnerParam()){
313 AliExternalTrackParam tmptrk[]={*(trk0->GetInnerParam()), *(trk1->GetInnerParam())};
316 printf("\n************************ raw ESD:\n");
317 AliTPCCosmicUtils::PrintTrackParam(100, &(tmptrk[0]));
318 AliTPCCosmicUtils::PrintTrackParam(101, &(tmptrk[1]));
323 Double_t xyz0[3], xyz1[3];
324 tmptrk[0].GetXYZ(xyz0);
325 tmptrk[1].GetXYZ(xyz1);
326 const TVector3 gpos0(xyz0), gpos1(xyz1);
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();
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()) ){
340 //track1 closer to mean alpha
342 if( !AliTPCCosmicUtils::RotateSafe(&(tmptrk[1]), meanalpha) ||
343 !AliTPCCosmicUtils::RotateSafe(&(tmptrk[0]), meanalpha+TMath::Pi()) ){
350 printf("\n************************ after rotation!!\n");
351 AliTPCCosmicUtils::PrintTrackParam(300, &(tmptrk[0]));
352 AliTPCCosmicUtils::PrintTrackParam(301, &(tmptrk[1]));
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){
370 for(Int_t ii=0; ii<2; ii++){
371 if(!AliTrackerBase::PropagateTrackToBxByBz(&(tmptrk[ii]), xTogo, AliTPCCosmicUtils::fgkMass, maxStep, rotateTo, maxSnp, eloss[ii])){
378 printf("\n************************ after dedx corr:\n");
379 AliTPCCosmicUtils::PrintTrackParam(200, &(tmptrk[0]));
380 AliTPCCosmicUtils::PrintTrackParam(201, &(tmptrk[1]));
385 //____________________________________________________________________________________
386 //____________________________________________________________________________________
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
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];
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]);
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]);
414 for(Int_t ii=0; ii<5; ii++){
415 if( TMath::Abs(fPull[ii]) > fCutPull[ii] ){
416 fErrFlagIsPair = 10+ii;
419 if( TMath::Abs(fDelta[ii]) > fCutDelta[ii] ){
420 fErrFlagIsPair = 20+ii;
431 Bool_t AliCosmicTracker::ESDtrackCut(AliESDtrack * trk, Double_t &findabler)
434 //cut on track quality (kink, TPCrefit, findable ratio) and require TPC seed
443 if(trk->GetKinkIndex(0)>0)
447 if(!trk->IsOn(AliESDtrack::kTPCrefit))
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
453 if(!trk->GetTPCNclsF())
456 findabler = (Double_t)trk->GetTPCNcls()/(Double_t) trk->GetTPCNclsF();
458 if(findabler < fgkCutFindable ){
463 if(!AliTPCCosmicUtils::GetTPCseed(trk))
469 void AliCosmicTracker::WriteStreamer(Int_t ntrk)
475 (*fStreamer)<<"CosmicTracker_Streamer"<<
477 "combup="<<&fCombTrackUp<<
478 "comblow="<<&fCombTrackLow<<
479 "track0="<<&fTrack0<<
480 "track1="<<&fTrack1<<
483 "isreuse="<<fIsReuse<<
484 "findabler="<<fFindableRatio<<
486 "chi2="<<fChi2PerCluster<<
494 "delta0="<<fDelta[0]<<
495 "delta1="<<fDelta[1]<<
496 "delta2="<<fDelta[2]<<
497 "delta3="<<fDelta[3]<<
498 "delta4="<<fDelta[4]<<