TPC/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