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 combined 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 combined 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   fCombinedTrackfit = 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 = fCombinedTrackfit->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 = fCombinedTrackfit->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::fgkNclsMin
  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<fgkCutLeverArm
  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 > fgkMaxChi2
 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   , fLeverArm(-999)
 136   , fFitNcls(-999), fMissNcls(-999), fPreChi2(-999), fFitLeverArm(-999), fImpactD(-999), fImpactZ(-999)
 137 {
 138   //
 139   //Constructor
 140   //
 141   fInnerClusterUp.SetXYZ(-999,-999,-999);
 142   fInnerClusterLow.SetXYZ(-999,-999,-999);
 143
 144   if(fDebugLevel>0)
 145     fStreamer = new TTreeSRedirector(Form("CombinedTrackfit_%s.root", tag.Data()));
 146
 147   SetRow(0,1);
 148   SetX(-1e10, 1e10);
 149 }
 150
 151 AliTPCCosmicTrackfit::~AliTPCCosmicTrackfit()
 152 {
 153   //
 154   //Destructor
 155   //
 156
 157   delete fStreamer;
 158
 159   delete fTrackparUp;
 160   delete fTrackparLow;
 161   delete fIniTrackparUp;
 162   delete fIniTrackparLow;
 163 }
 164
 165 Bool_t AliTPCCosmicTrackfit::CombineESDtracks(AliESDtrack * &trk0, AliESDtrack *&trk1)
 166 {
 167   //
 168   //Get TPCseeds from the 2 ESDtracks, swap TPCseeds and ESDTracks (if necessary) according to y (0:upper 1:lower), perform trackfit using TPCseeds
 169   //if fStatus==0, i.e. combine is successful, swap of the ESDtracks is kept since pointer *& is used
 170   //
 171
 172   IniCombineESDtracks();
 173
 174   if(!GetTPCseeds(trk0, trk1)){
 175     return kFALSE;
 176   }
 177
 178   CombineTPCseeds();
 179
 180   if(fStatus == 0){
 181     if(fKswap){
 182       AliESDtrack * tmptrk = trk0;
 183       trk0 = trk1;
 184       trk1 = tmptrk;
 185     }
 186     return kTRUE;
 187   }
 188   else
 189     return kFALSE;
 190 }
 191
 192 Bool_t AliTPCCosmicTrackfit::CombineTPCseedsFast(AliTPCseed * tpcseeds[], const AliExternalTrackParam * trkpars[])
 193 {
 194   //
 195   //do combined track fit for given TPC seeds and initial trackpar, [0]: upper, [1]: lower
 196   //
 197
 198   IniCombineESDtracks();
 199
 200   fSeedUp = tpcseeds[0];
 201   fSeedLow = tpcseeds[1];
 202
 203   fTrackparUp =  new AliExternalTrackParam(*(trkpars[0]));
 204   fTrackparLow = new AliExternalTrackParam(*(trkpars[1]));
 205
 206   AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
 207   const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
 208   TTreeSRedirector * debugstreamer = 0x0;
 209   if(fDebugLevel&2){
 210     debugstreamer = fStreamer;
 211   }
 212
 213   AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
 214
 215   Update();
 216
 217   if(fStatus == 0){
 218     return kTRUE;
 219   }
 220   else
 221     return kFALSE;
 222 }
 223
 224 Bool_t AliTPCCosmicTrackfit::CombineTPCseeds(AliTPCseed * &seed0, AliTPCseed *&seed1)
 225 {
 226   //
 227   //same as AliTPCCosmicTrackfit::CombineESDtracks, except that the seeds are passed in from outside, which can be still unordered
 228   //if fStatus==0, i.e. combine is successful, swap of the TPCseeds is kept since pointer *& is used
 229   //
 230   IniCombineESDtracks();
 231
 232   fSeedUp  = seed0;
 233   fSeedLow = seed1;
 234
 235   CombineTPCseeds();
 236
 237   if(fStatus==0){
 238     if(fKswap){
 239       AliTPCseed * tmpseed = seed0;
 240       seed0 = seed1;
 241       seed1 = tmpseed;
 242     }
 243     return kTRUE;
 244   }
 245   else
 246     return kFALSE;
 247 }
 248
 249 void AliTPCCosmicTrackfit::Print() const
 250 {
 251   //
 252   //print out variable values
 253   //
 254   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);
 255 }
 256
 257 /*
 258 Double_t AliTPCCosmicTrackfit::ImpactParameter2D() const
 259 {
 260   //
 261   //calculate the 2D-impactparameter from (0,0)
 262   //
 263
 264   const TVector3 p0(0,0,0);
 265   const TVector3 l1(fInnerClusterUp.X(),  fInnerClusterUp.Y(),  0);
 266   const TVector3 l2(fInnerClusterLow.X(), fInnerClusterLow.Y(), 0);
 267
 268   return AliTPCCosmicUtils::Point2LineDist(p0, l1, l2);
 269 }
 270
 271 Double_t AliTPCCosmicTrackfit::ImpactParameter3D() const
 272 {
 273   //
 274   //calculate the 3D-impactparameter from (0,0,0)
 275   //
 276
 277   const TVector3 p0(0,0,0);
 278   const TVector3 l1(fInnerClusterUp);
 279   const TVector3 l2(fInnerClusterLow);
 280
 281   return AliTPCCosmicUtils::Point2LineDist(p0, l1, l2);
 282 }
 283 */
 284
 285 Double_t AliTPCCosmicTrackfit::MinPhi() const
 286 {
 287   //
 288   //the smaller phi of the two tracks w.r.t. horizon
 289   //
 290   Double_t fsp[] = {TMath::Abs(TMath::Sin(fTrackparUp->Phi())), TMath::Abs(TMath::Sin(fTrackparLow->Phi()))};;
 291   return TMath::ASin(TMath::Min(fsp[0], fsp[1])) * TMath::RadToDeg();
 292 }
 293 //===================================================================================================
 294 //===================================================================================================
 295
 296 void AliTPCCosmicTrackfit::IniCombineESDtracks()
 297 {
 298   //
 299   //initialization, for reuse of the same AliTPCCosmicTrackfit instance
 300   //
 301
 302   fSeedUp = 0x0;
 303   fSeedLow = 0x0;
 304   delete fTrackparUp;
 305   delete fTrackparLow;
 306   fTrackparUp = 0x0;
 307   fTrackparLow = 0x0;
 308
 309   delete fIniTrackparUp;
 310   delete fIniTrackparLow;
 311   fIniTrackparUp = 0x0;
 312   fIniTrackparLow = 0x0;
 313
 314   fStatus = 0;
 315   fKswap = kFALSE;
 316 }
 317
 318 void AliTPCCosmicTrackfit::CombineTPCseeds()
 319 {
 320   //
 321   //do combined trackfit using TPCseeds
 322   //
 323
 324   if(
 325      !CheckNcls()
 326      || !AnaSeeds() // fSeedUp/Low swapped here! check by runTest.C
 327      || !CheckLeverArm()
 328      )
 329     return;
 330
 331   //AliExternalTrackParam object created
 332   fTrackparUp  = AliTPCCosmicUtils::MakeSeed(fSeedUp);
 333   fTrackparLow = AliTPCCosmicUtils::MakeSeed(fSeedLow);
 334   if(!fTrackparUp || !fTrackparLow){
 335     fStatus = kFailMakeSeed;
 336     return;
 337   }
 338
 339   fIniTrackparUp = new AliExternalTrackParam(*fTrackparUp);
 340   fIniTrackparLow = new AliExternalTrackParam(*fTrackparLow);
 341
 342   AliExternalTrackParam * trackPars[]={fTrackparUp, fTrackparLow};
 343   const AliTPCseed *seeds[]={fSeedUp, fSeedLow};
 344   TTreeSRedirector * debugstreamer = 0x0;
 345   if(fDebugLevel&2){
 346     debugstreamer = fStreamer;
 347   }
 348
 349   AliTPCCosmicUtils::CombinedFit(trackPars, seeds, fRowStartShift, fRowStep, fXMin, fXMax, fFitNcls, fMissNcls, fPreChi2, fFitLeverArm, fImpactD, fImpactZ, debugstreamer);
 350
 351   Update();
 352
 353   return;
 354 }
 355
 356 void AliTPCCosmicTrackfit::Update()
 357 {
 358   //
 359   //Update variables depending on the fit result
 360   //
 361
 362   if(fMissNcls || fFitNcls==0){
 363     fStatus = kFailPropagation;
 364     return;
 365   }
 366
 367   fPreChi2 /= fFitNcls;
 368   if(fPreChi2>fgkMaxChi2){
 369     fStatus = kFailChi2;
 370     return;
 371   }
 372
 373   if(fImpactD<0){
 374     fStatus = kFailImpact;
 375     return;
 376   }
 377
 378   if( fStatus == 0 && (fDebugLevel&1) ){
 379     Double_t momup  = fTrackparUp->P();
 380     Double_t momlow = fTrackparLow->P();
 381     Double_t ptup   = fTrackparUp->Pt();
 382     Double_t ptlow  = fTrackparLow->Pt();
 383
 384     (*fStreamer)<<"TrackProp"<<
 385       "Tup.="<<fTrackparUp<<
 386       "Tlow.="<<fTrackparLow<<
 387       "icup.="<<&fInnerClusterUp<<
 388       "iclow.="<<&fInnerClusterLow<<
 389       "leverarm="<<fLeverArm<<
 390       "ncl="<<fFitNcls<<
 391       "nmiss="<<fMissNcls<<
 392       "chi2="<<fPreChi2<<
 393       "momup="<<  momup <<
 394       "momlow="<< momlow <<
 395       "ptup="<<   ptup <<
 396       "ptlow="<<  ptlow <<
 397       "\n";
 398   }
 399 }
 400
 401 Bool_t AliTPCCosmicTrackfit::CheckLeverArm()
 402 {
 403   //
 404   //if lever arm is too short, no need to use combined track fit.
 405   //On the other hand, short lever arm from two tracks mostly means they are fake pairs.
 406   //lever arm extents over one quadrant, e.g. (0,250)-(250,0): 250*sqrt(2)~350
 407   //
 408   if(fLeverArm<fgkCutLeverArm){
 409     fStatus = kFailLeverArm;
 410     return kFALSE;
 411   }
 412   else
 413     return kTRUE;
 414 }
 415
 416 Bool_t AliTPCCosmicTrackfit::AnaSeeds()
 417 {
 418   //
 419   //swap seeds (if necessary) so that (y of fSeedUp) > (y of fSeedLow)
 420   //
 421
 422   //---------------------------------- navigate through all clusters ----------------------------------
 423   AliTPCseed ** seeds[]={&fSeedUp, &fSeedLow};
 424
 425   //min, max according to y
 426   TVector3 singlemin[2], singlemax[2];
 427   for(Int_t ii=0; ii<2; ii++){
 428     singlemin[ii].SetXYZ( 1e10,  1e10,  1e10);
 429     singlemax[ii].SetXYZ(-1e10, -1e10, -1e10);
 430   }
 431
 432   for(Int_t itrk=0; itrk<2; itrk++){
 433     for(Int_t irow=0; irow<AliTPCCosmicUtils::fgkNRow; irow++){
 434       const AliTPCclusterMI * cls = (*seeds[itrk])->GetClusterPointer(irow);
 435       if(!cls)
 436         continue;
 437
 438       Float_t xyz[3]={-999,-999,-999};
 439       cls->GetGlobalXYZ(xyz);
 440       if(xyz[1]<singlemin[itrk].Y()){
 441         singlemin[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
 442       }
 443       if(xyz[1]>singlemax[itrk].Y()){
 444         singlemax[itrk].SetXYZ(xyz[0], xyz[1], xyz[2]);
 445       }
 446     }
 447   }
 448
 449   //--------------------------------
 450
 451   //kpass true if y of the two seeds clearly separate: min of one > max of the other
 452   Bool_t kpass = kFALSE;
 453
 454   fInnerClusterUp.SetXYZ(-999,-999,-999);
 455   fInnerClusterLow.SetXYZ(-999,-999,-999);
 456   TVector3 combinedmin, combinedmax;
 457   if(singlemin[0].Y()> singlemax[1].Y()){
 458     fInnerClusterUp  = singlemin[0];
 459     fInnerClusterLow = singlemax[1];
 460
 461     //no need to swap
 462     fKswap = kFALSE;
 463
 464     kpass = kTRUE;
 465
 466     combinedmax = singlemax[0];
 467     combinedmin = singlemin[1];
 468   }
 469   else if(singlemin[1].Y()> singlemax[0].Y()){
 470     fInnerClusterUp  = singlemin[1];
 471     fInnerClusterLow = singlemax[0];
 472
 473     //have to be swapped
 474     fKswap = kTRUE;
 475     AliTPCseed *tmp=*(seeds[0]);
 476     *(seeds[0])=*(seeds[1]);
 477     *(seeds[1])=tmp;
 478
 479     kpass = kTRUE;
 480
 481     combinedmax = singlemax[1];
 482     combinedmin = singlemin[0];
 483   }
 484   else
 485     kpass = kFALSE;
 486
 487   const TVector3 comdelta = combinedmax-combinedmin;
 488   fLeverArm = comdelta.Pt();
 489
 490   if(!kpass){
 491     fStatus = kFailSwapSeeds;
 492     return kFALSE;
 493   }
 494   else
 495     return kTRUE;
 496 }
 497
 498 Bool_t AliTPCCosmicTrackfit::CheckNcls()
 499 {
 500   //
 501   //check number of clusters in TPCseed, for too small number MakeSeed will fail
 502   //
 503   if( fSeedUp->GetNumberOfClusters()<AliTPCCosmicUtils::fgkNclsMin || fSeedLow->GetNumberOfClusters()<AliTPCCosmicUtils::fgkNclsMin ){
 504     fStatus = kFailNclsMin;
 505     return kFALSE;
 506   }
 507   else
 508     return kTRUE;
 509 }
 510
 511 Bool_t AliTPCCosmicTrackfit::GetTPCseeds(const AliESDtrack *trk0,  const AliESDtrack *trk1)
 512 {
 513   //
 514   //Get TPC seeds from ESDfriendTrack
 515   //
 516   fSeedUp  = AliTPCCosmicUtils::GetTPCseed(trk0);
 517   fSeedLow = AliTPCCosmicUtils::GetTPCseed(trk1);
 518
 519   if(!fSeedUp || !fSeedLow){
 520     fStatus = kFailGetTPCseeds;
 521     return kFALSE;
 522   }
 523
 524   return kTRUE;
 525 }
 526