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