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 }