PWG1/AliESDRecInfo.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 //  Time Projection Chamber                                                  //
  20 //  Comparison macro for reconstructed tracks                                  //
  21 //  responsible:
  22 //  marian.ivanov@cern.ch                                                    //
  23 //
  24 //
  25
  26
  27
  28
  29
  30 //ROOT includes
  31 #include "Rtypes.h"
  32 //
  33 //ALIROOT includes
  34 //
  35 #include "AliESDtrack.h"
  36 #include "AliTracker.h"
  37 #include "AliTPCParam.h"
  38 #include "AliTrackReference.h"
  39 #include "AliTPCParamSR.h"
  40 #include "AliESDfriend.h"
  41 #include "AliESDtrack.h"
  42 #include "AliTPCseed.h"
  43 #include "AliITStrackMI.h"
  44 #include "AliTRDtrackV1.h"
  45 #include "AliMCInfo.h"
  46 #include "AliESDRecInfo.h"
  47
  48
  49
  50 ClassImp(AliESDRecInfo)
  51
  52
  53
  54
  55 AliTPCParam * GetTPCParam(){
  56   AliTPCParamSR * par = new AliTPCParamSR;
  57   par->Update();
  58   return par;
  59 }
  60
  61
  62
  63
  64 AliESDRecInfo::AliESDRecInfo():
  65   fITSOn(0),           // ITS refitted inward
  66   fTRDOn(0),           // ITS refitted inward
  67   fDeltaP(0),          //delta of momenta
  68   fSign(0),           // sign
  69   fReconstructed(0),         //flag if track was reconstructed
  70   fFake(0),             // fake track
  71   fMultiple(0),         // number of reconstructions
  72   fTPCOn(0),           // TPC refitted inward
  73   fBestTOFmatch(0),        //best matching between times
  74   fESDtrack(0),        // esd track
  75   fTrackF(0),      // friend track
  76   fTPCtrack(0),        // tpc track
  77   fITStrack(0),        // its track
  78   fTRDtrack(0)        // trd track
  79 {
  80   //
  81   //  default constructor
  82   //
  83 }
  84
  85
  86 AliESDRecInfo::AliESDRecInfo(const AliESDRecInfo& recinfo):
  87   TObject(),
  88   fITSOn(0),           // ITS refitted inward
  89   fTRDOn(0),           // ITS refitted inward
  90   fDeltaP(0),          //delta of momenta
  91   fSign(0),           // sign
  92   fReconstructed(0),         //flag if track was reconstructed
  93   fFake(0),             // fake track
  94   fMultiple(0),         // number of reconstructions
  95   fTPCOn(0),           // TPC refitted inward
  96   fBestTOFmatch(0),        //best matching between times
  97   fESDtrack(0),        // esd track
  98   fTrackF(0),      // friend track
  99   fTPCtrack(0),        // tpc track
 100   fITStrack(0),        // its track
 101   fTRDtrack(0)        // trd track
 102 {
 103   //
 104   //
 105   //
 106   memcpy(this,&recinfo, sizeof(recinfo));
 107   fESDtrack=0; fTrackF=0; fTPCtrack=0;fITStrack=0;fTRDtrack=0;
 108   SetESDtrack(recinfo.GetESDtrack());
 109 }
 110
 111
 112 AliESDRecInfo& AliESDRecInfo::operator=(const AliESDRecInfo& info) {
 113   //
 114   // Assignment operator
 115   //
 116   this->~AliESDRecInfo();
 117   new (this) AliESDRecInfo(info);
 118   return *this;
 119 }
 120
 121
 122
 123 AliESDRecInfo::~AliESDRecInfo()
 124
 125 {
 126   //
 127   //  destructor
 128   //
 129   if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
 130   if (fTrackF)   { delete fTrackF;   fTrackF=0;}
 131   if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
 132   if (fITStrack) { delete fITStrack; fITStrack=0;}
 133   if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
 134
 135 }
 136
 137
 138
 139 void AliESDRecInfo::Reset()
 140 {
 141   //
 142   // reset info
 143   //
 144   fMultiple =0;
 145   fFake     =0;
 146   fReconstructed=0;
 147   if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
 148   if (fTrackF)   { delete fTrackF;   fTrackF=0;}
 149   if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
 150   if (fITStrack) { delete fITStrack; fITStrack=0;}
 151   if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
 152 }
 153
 154 void AliESDRecInfo::SetESDtrack(const AliESDtrack *track){
 155   //
 156   //
 157   //
 158   if (fESDtrack) delete fESDtrack;
 159   fESDtrack = (AliESDtrack*)track->Clone();
 160   if (track->GetFriendTrack()){
 161     if (fTrackF) delete fTrackF;
 162     fTrackF = (AliESDfriendTrack*)track->GetFriendTrack()->Clone();
 163     Int_t icalib=0;
 164     TObject *cobject=0;
 165     //
 166     while (fTrackF->GetCalibObject(icalib)){
 167       cobject=fTrackF->GetCalibObject(icalib);
 168       if (dynamic_cast<AliTPCseed*>(cobject)){
 169                                 if (fTPCtrack) delete fTPCtrack;
 170                                 fTPCtrack = (AliTPCseed*)(dynamic_cast<AliTPCseed*>(cobject))->Clone();
 171       } else if (dynamic_cast<AliTRDtrackV1*>(cobject)){
 172                                 if (fTRDtrack) delete fTRDtrack;
 173                                 fTRDtrack = (AliTRDtrackV1*)(dynamic_cast<AliTRDtrackV1*>(cobject))->Clone();
 174       }
 175       icalib++;
 176     }
 177   }
 178
 179 }
 180
 181 void  AliESDRecInfo::UpdatePoints(AliESDtrack*track)
 182 {
 183   //
 184   //
 185   Int_t iclusters[200];
 186   Float_t density[160];
 187   for (Int_t i=0;i<160;i++) density[i]=-1.;
 188   fTPCPoints[0]= 160;
 189   fTPCPoints[1] = -1;
 190   //
 191   if (fTPCPoints[0]<fTPCPoints[1]) return;
 192   //  Int_t nclusters=track->GetTPCclusters(iclusters);
 193
 194   Int_t ngood=0;
 195   Int_t undeff=0;
 196   Int_t nall =0;
 197   Int_t range=20;
 198   for (Int_t i=0;i<160;i++){
 199     Int_t last = i-range;
 200     if (nall<range) nall++;
 201     if (last>=0){
 202       if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
 203       if (iclusters[last]==-1) undeff--;
 204     }
 205     if (iclusters[i]>0&& (iclusters[i]&0x8000)==0)   ngood++;
 206     if (iclusters[i]==-1) undeff++;
 207     if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
 208   }
 209   Float_t maxdens=0;
 210   Int_t indexmax =0;
 211   for (Int_t i=0;i<160;i++){
 212     if (density[i]<0) continue;
 213     if (density[i]>maxdens){
 214       maxdens=density[i];
 215       indexmax=i;
 216     }
 217   }
 218   //
 219   //max dens point
 220   fTPCPoints[3] = maxdens;
 221   fTPCPoints[1] = indexmax;
 222   //
 223   // last point
 224   for (Int_t i=indexmax;i<160;i++){
 225     if (density[i]<0) continue;
 226     if (density[i]<maxdens/2.) {
 227       break;
 228     }
 229     fTPCPoints[2]=i;
 230   }
 231   //
 232   // first point
 233   for (Int_t i=indexmax;i>0;i--){
 234     if (density[i]<0) continue;
 235     if (density[i]<maxdens/2.) {
 236       break;
 237     }
 238     fTPCPoints[0]=i;
 239   }
 240   //
 241   // Density at the last 30 padrows
 242   //
 243   //
 244   nall  = 0;
 245   ngood = 0;
 246   for (Int_t i=159;i>0;i--){
 247     if (iclusters[i]==-1) continue; //dead zone
 248     nall++;
 249     if (iclusters[i]>0)   ngood++;
 250     if (nall>20) break;
 251   }
 252   fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
 253   //
 254   if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
 255
 256
 257 }
 258
 259 //
 260 //
 261 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed)
 262 {
 263   //
 264   //calculates derived variables
 265   //
 266   UpdatePoints(fESDtrack);
 267   UpdateStatus(info,reconstructed);
 268   UpdateITS(info);
 269   UpdateTPC(info);
 270   UpdateTOF(info);
 271 }
 272
 273
 274 void  AliESDRecInfo::UpdateStatus(AliMCInfo* info, Bool_t reconstructed){
 275   //
 276   // Interpret bit mask flags
 277   //
 278   for (Int_t i=0;i<4;i++) fStatus[i] =0;
 279   fReconstructed = kFALSE;
 280   fTPCOn = kFALSE;
 281   fITSOn = kFALSE;
 282   fTRDOn = kFALSE;
 283   if (reconstructed==kFALSE) return;
 284
 285   fLabels[0] = info->fLabel;
 286   fLabels[1] = info->fPrimPart;
 287   fReconstructed = kTRUE;
 288   fTPCOn = ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
 289   fITSOn = ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
 290   fTRDOn = ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
 291   //
 292   //
 293   if ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0){
 294     fStatus[1] =3;
 295   }
 296   else{
 297     if ((fESDtrack->GetStatus()&AliESDtrack::kTPCout)>0){
 298       fStatus[1] =2;
 299     }
 300     else{
 301       if ((fESDtrack->GetStatus()&AliESDtrack::kTPCin)>0)
 302         fStatus[1]=1;
 303     }
 304   }
 305   //
 306   if ((fESDtrack->GetStatus()&AliESDtrack::kITSout)>0){
 307     fStatus[0] =2;
 308   }
 309   else{
 310     if ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0){
 311       fStatus[0] =1;
 312     }
 313     else{
 314       fStatus[0]=0;
 315     }
 316   }
 317   //
 318   //
 319   if ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0){
 320     fStatus[2] =2;
 321   }
 322   else{
 323     if ((fESDtrack->GetStatus()&AliESDtrack::kTRDout)>0){
 324       fStatus[2] =1;
 325     }
 326   }
 327   if ((fESDtrack->GetStatus()&AliESDtrack::kTRDStop)>0){
 328     fStatus[2] =10;
 329   }
 330 }
 331
 332 void AliESDRecInfo::UpdateTOF(AliMCInfo* info){
 333   //
 334   // Update TOF related comparison information
 335   //
 336   fBestTOFmatch=1000;
 337   if (((fESDtrack->GetStatus()&AliESDtrack::kTOFout)>0)){
 338     //
 339     // best tof match
 340     Double_t times[5];
 341     fESDtrack->GetIntegratedTimes(times);
 342     for (Int_t i=0;i<5;i++){
 343       if ( TMath::Abs(fESDtrack->GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
 344         fBestTOFmatch = fESDtrack->GetTOFsignal()-times[i];
 345       }
 346     }
 347     Int_t toflabel[3];
 348     fESDtrack->GetTOFLabel(toflabel);
 349     Bool_t toffake=kTRUE;
 350     Bool_t tofdaughter=kFALSE;
 351     for (Int_t i=0;i<3;i++){
 352       if (toflabel[i]<0) continue;
 353       if (toflabel[i]== TMath::Abs(fESDtrack->GetLabel()))  toffake=kFALSE;
 354       if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE;  // decay product of original particle
 355       fStatus[3]=1;
 356     }
 357     if (toffake) fStatus[3] =3;       //total fake
 358     if (tofdaughter) fStatus[3]=2;    //fake because of decay
 359   }else{
 360     fStatus[3]=0;
 361   }
 362 }
 363
 364
 365
 366 void AliESDRecInfo::UpdateITS(AliMCInfo* info){
 367   //
 368   // Update ITS related comparison information
 369   //
 370   fITSinP0[0]=info->fParticle.Px();
 371   fITSinP0[1]=info->fParticle.Py();
 372   fITSinP0[2]=info->fParticle.Pz();
 373   fITSinP0[3]=info->fParticle.Pt();
 374   //
 375   fITSinR0[0]=info->fParticle.Vx();
 376   fITSinR0[1]=info->fParticle.Vy();
 377   fITSinR0[2]=info->fParticle.Vz();
 378   fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
 379   fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
 380   //
 381   //
 382   if (fITSinP0[3]>0.0000001){
 383     fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
 384     fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
 385   }
 386   if (fITSOn){
 387     // ITS
 388     Double_t param[5],x;
 389     fESDtrack->GetExternalParameters(x,param);
 390     //    fESDtrack->GetConstrainedExternalParameters(x,param);
 391     Double_t cov[15];
 392     fESDtrack->GetExternalCovariance(cov);
 393     //fESDtrack->GetConstrainedExternalCovariance(cov);
 394     if (TMath::Abs(param[4])<0.0000000001) return;
 395
 396     fESDtrack->GetXYZ(fITSinR1);
 397     fESDtrack->GetPxPyPz(fITSinP1);
 398     fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
 399     //
 400     fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
 401     fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
 402     //
 403     //
 404     if (fITSinP1[3]>0.0000001){
 405       fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
 406       fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);
 407     }
 408     //
 409     //
 410     fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3];  //delta rfi
 411     fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
 412     fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]);              //delta z
 413     fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
 414     fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
 415     fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
 416     fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
 417     fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
 418     fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);
 419     Double_t sign = (param[4]>0) ? 1:-1;
 420     fSign = sign;
 421     fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
 422   }
 423 }
 424
 425 void AliESDRecInfo::UpdateTPC(AliMCInfo* info){
 426   //
 427   //
 428   // Update TPC related information
 429   //
 430   AliTrackReference * ref = &(info->fTrackRef);
 431   fTPCinR0[0] = info->fTrackRef.X();
 432   fTPCinR0[1] = info->fTrackRef.Y();
 433   fTPCinR0[2] = info->fTrackRef.Z();
 434   fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
 435   fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
 436   //
 437   fTPCinP0[0] = ref->Px();
 438   fTPCinP0[1] = ref->Py();
 439   fTPCinP0[2] = ref->Pz();
 440   fTPCinP0[3] = ref->Pt();
 441   fTPCinP0[4] = ref->P();
 442   fDeltaP     = (ref->P()-info->fParticle.P())/info->fParticle.P();
 443   //
 444   //
 445   if (fTPCinP0[3]>0.0000001){
 446     //
 447     fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
 448     fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
 449   }
 450   //
 451   if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
 452     //TPC
 453     fESDtrack->GetInnerXYZ(fTPCinR1);
 454     fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
 455     fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
 456     fESDtrack->GetInnerPxPyPz(fTPCinP1);
 457     fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
 458     fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
 459     //
 460     //
 461     if (fTPCinP1[3]>0.000000000000001){
 462       fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
 463       fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
 464     }
 465     Double_t cov[15], param[5],x, alpha;
 466     fESDtrack->GetInnerExternalCovariance(cov);
 467     fESDtrack->GetInnerExternalParameters(alpha, x,param);
 468     if (x<50) return ;
 469     //
 470     fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
 471     fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
 472     fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
 473     fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
 474     fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
 475     fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
 476     fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
 477     fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
 478     fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
 479     Double_t sign = (param[4]>0)? 1.:-1;
 480     fSign =sign;
 481     fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14]));
 482   }
 483 }
 484
 485
 486
 487
 488