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