HBTAN/AliHBTTrackPoints.cxx

   1 #include "AliHBTTrackPoints.h"
   2 //_________________________________
   3 ////////////////////////////////////////////////////////////
   4 //                                                        //
   5 // class AliHBTTrackPoints                                //
   6 //                                                        //
   7 // used by Anti-Merging cut                               //
   8 // contains set of poits the lay on track trajectory      //
   9 // according to reconstructed track parameters -          //
  10 // NOT CLUSTERS POSITIONS!!!                              //
  11 // Anti-Merging cut is applied only on tracks coming from //
  12 // different events (that are use to fill deniminators)   //
  13 //                                                        //
  14 ////////////////////////////////////////////////////////////
  15 #include "AliTrackReference.h"
  16
  17 #include "AliTPCtrack.h"
  18 #include "AliESDtrack.h"
  19
  20 #include <TMath.h>
  21 #include <TClonesArray.h>
  22
  23 ClassImp(AliHBTTrackPoints)
  24
  25 Int_t AliHBTTrackPoints::fgDebug = 0;
  26 AliHBTTrackPoints::AliHBTTrackPoints():
  27  fN(0),
  28  fX(0x0),
  29  fY(0x0),
  30  fZ(0x0)
  31 {
  32   //constructor
  33 }
  34 /***************************************************************/
  35
  36 AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
  37  fN(n),
  38  fX(new Float_t[fN]),
  39  fY(new Float_t[fN]),
  40  fZ(new Float_t[fN])
  41 {
  42   //constructor
  43   //mf - magnetic field in kG - needed to calculated curvature out of Pt
  44   //r0 - starting radius
  45   //dr - calculate points every dr cm, default every 30cm
  46   if (track == 0x0)
  47    {
  48      Error("AliHBTTrackPoints","ESD track is null");
  49      fN = 0;
  50      delete [] fX;
  51      delete [] fY;
  52      delete [] fZ;
  53      fX = fY = fZ = 0x0;
  54      return;
  55    }
  56
  57   Double_t x;
  58   Double_t par[5];
  59   track->GetInnerExternalParameters(x,par);     //get properties of the track
  60   if (par[4] == 0)
  61    {
  62      Error("AliHBTTrackPoints","This ESD track does not contain TPC information");
  63      return;
  64    }
  65
  66   if (mf == 0.0)
  67    {
  68      Error("AliHBTTrackPoints","Zero Magnetic field passed as parameter.");
  69      return;
  70    }
  71
  72   Double_t alpha = track->GetInnerAlpha();
  73   Double_t cc = 1000./0.299792458/mf;//conversion constant
  74   Double_t c=par[4]/cc;
  75
  76   MakePoints(dr,r0,x,par,c,alpha);
  77
  78 }
  79 /***************************************************************/
  80
  81 AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
  82  fN(n),
  83  fX(new Float_t[fN]),
  84  fY(new Float_t[fN]),
  85  fZ(new Float_t[fN])
  86 {
  87   //constructor
  88   //r0 starting radius
  89   //dr - calculate points every dr cm, default every 30cm
  90   if (track == 0x0)
  91    {
  92      Error("AliHBTTrackPoints","TPC track is null");
  93      fN = 0;
  94      delete [] fX;
  95      delete [] fY;
  96      delete [] fZ;
  97      fX = fY = fZ = 0x0;
  98      return;
  99    }
 100   track->PropagateTo(r0);
 101
 102  //* This formation is now fixed in the following way:                         *
 103  //*      external param0:   local Y-coordinate of a track (cm)                *
 104  //*      external param1:   local Z-coordinate of a track (cm)                *
 105  //*      external param2:   local sine of the track momentum azimuth angle    *
 106  //*      external param3:   tangent of the track momentum dip angle           *
 107  //*      external param4:   1/pt (1/(GeV/c))                                  *
 108
 109   Double_t x = 0;
 110   Double_t par[5];
 111   track->GetExternalParameters(x,par);     //get properties of the track
 112
 113   Double_t alpha = track->GetAlpha();
 114   Double_t c=track->GetC();
 115   MakePoints(dr,r0,x,par,c,alpha);
 116 }
 117
 118 void AliHBTTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
 119 {
 120   //Calculates points starting at radius r0
 121   //spacing every dr (in radial direction)
 122   // according to track parameters
 123   // x - position in sector local reference frame. x i parallel to R and sector is symmetric with respect to x
 124   // par - track external parameters; array with 5 elements;  look at AliTPCtrack.h or AliESDtrack.h for their meaning
 125   // c - track curvature
 126   // alpha - sector's rotation angle (phi) == angle needed for local to global transformation
 127
 128   Double_t y = par[0];
 129   Double_t z0 = par[1];
 130
 131   Double_t phi0local = TMath::ATan2(y,x);
 132   Double_t phi0global = phi0local + alpha;
 133
 134   if (phi0local<0) phi0local+=2*TMath::Pi();
 135   if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
 136
 137   if (phi0global<0) phi0global+=2*TMath::Pi();
 138   if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
 139
 140   Double_t r = TMath::Hypot(x,y);
 141
 142
 143   if (GetDebug() > 9)
 144     Info("AliHBTTrackPoints","Radius0 %f, Real Radius %f",r0,r);
 145
 146   if (GetDebug() > 5)
 147     Info("AliHBTTrackPoints","Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local);
 148
 149   Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
 150
 151   //this calculattions are assuming new (current) model
 152   Double_t tmp = par[2];
 153   tmp = 1. - tmp*tmp;
 154   tmp = c*y + TMath::Sqrt(tmp);
 155   Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
 156
 157   //Here is old model Cold=Cnew/2.
 158   Double_t dcasq = dca*dca;
 159   Double_t c2 = c/2.;
 160   Double_t cst1 = (1.+c2*dca)*dca;//first constant
 161   Double_t cst2 = 1. + 2.*c2*dca;//second constant
 162
 163   Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
 164   Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
 165
 166   for(Int_t i = 0; i<fN; i++)
 167    {
 168      Double_t rc = r0 + i*dr;
 169      Double_t ftmp = (c2*rc + cst1/rc)/cst2;
 170      if (ftmp > 1.0)
 171       {
 172         Warning("AliHBTTrackPoints","ASin argument > 1 %f:",ftmp);
 173         ftmp=1.0;
 174       }
 175      else if (ftmp < -1.0)
 176       {
 177         Warning("AliHBTTrackPoints","ASin argument < -1 %f:",ftmp);
 178         ftmp=-1.0;
 179       }
 180
 181      Double_t factorPhi = TMath::ASin( ftmp );//factor phi od rc
 182      Double_t phi = phi0global + factorPhi - factorPhi0;
 183
 184      ftmp = (rc*rc-dcasq)/cst2;
 185      if (ftmp < 0.0)
 186       {
 187         Warning("AliHBTTrackPoints","Sqrt argument < 0: %f",ftmp);
 188         ftmp=1.0;
 189       }
 190
 191      ftmp = c2*TMath::Sqrt(ftmp);
 192      if (ftmp > 1.0)
 193       {
 194         Warning("AliHBTTrackPoints","ASin argument > 1: %f",ftmp);
 195         ftmp=1.0;
 196       }
 197      else if (ftmp < -1.0)
 198       {
 199         Warning("AliHBTTrackPoints","ASin argument < -1: %f",ftmp);
 200         ftmp=-1.0;
 201       }
 202      Double_t factorZ = TMath::ASin(ftmp)*par[3]/c2;
 203      fZ[i] = z0 + factorZ - factorZ0;
 204      fX[i] = rc*TMath::Cos(phi);
 205      fY[i] = rc*TMath::Sin(phi);
 206
 207      if ( GetDebug() > 2 )
 208       {
 209         Info("AliHBTTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
 210              fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i]));
 211       }
 212    }
 213 }
 214 /***************************************************************/
 215
 216 AliHBTTrackPoints::~AliHBTTrackPoints()
 217 {
 218   //destructor
 219   delete [] fX;
 220   delete [] fY;
 221   delete [] fZ;
 222 }
 223 /***************************************************************/
 224
 225 void AliHBTTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
 226 {
 227   //returns position at point n
 228   if ((n<0) || (n>fN))
 229    {
 230      Error("PositionAt","Point %d out of range",n);
 231      return;
 232    }
 233
 234   x = fX[n];
 235   y = fY[n];
 236   z = fZ[n];
 237   if ( GetDebug() > 1 )
 238     {
 239       Info("AliHBTTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
 240     }
 241 }
 242 /***************************************************************/
 243
 244 Double_t AliHBTTrackPoints::AvarageDistance(const AliHBTTrackPoints& tr)
 245 {
 246   //returns the aritmethic avarage distance between two tracks
 247   if (fN != tr.fN)
 248    {
 249      Error("AvarageDistance","Number of points is not equal");
 250      return -1;
 251    }
 252   if ( (fN <= 0) || (tr.fN <=0) )
 253    {
 254      Error("AvarageDistance","One of tracks is empty");
 255      return -1;
 256    }
 257
 258   Double_t sum = 0;
 259   for (Int_t i = 0; i<fN; i++)
 260    {
 261      if (GetDebug()>9)
 262       {
 263 //       Float_t r1sq = fX[i]*fX[i]+fY[i]*fY[i];
 264 //       Float_t r2sq = tr.fX[i]*tr.fX[i]+tr.fY[i]*tr.fY[i];
 265        Float_t r1sq = TMath::Hypot(fX[i],fY[i]);
 266        Float_t r2sq = TMath::Hypot(tr.fX[i],tr.fY[i]);
 267        Info("AvarageDistance","radii: %f %f",r1sq,r2sq);
 268       }
 269
 270
 271      Double_t dx = fX[i]-tr.fX[i];
 272      Double_t dy = fY[i]-tr.fY[i];
 273      Double_t dz = fZ[i]-tr.fZ[i];
 274      sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 275
 276      if (GetDebug()>0)
 277       {
 278        Info("AvarageDistance","Diff: x ,y z: %f , %f, %f",dx,dy,dz);
 279        Info("AvarageDistance","xxyyzz %f %f %f %f %f %f",
 280             fX[i],tr.fX[i],fY[i],tr.fY[i],fZ[i],tr.fZ[i]);
 281       }
 282    }
 283
 284   Double_t retval = sum/((Double_t)fN);
 285   if ( GetDebug() )
 286     {
 287       Info("AvarageDistance","Avarage distance is %f.",retval);
 288     }
 289   return retval;
 290 }
 291 /***************************************************************/
 292 /***************************************************************/
 293 /***************************************************************/
 294 /***************************************************************/
 295 /***************************************************************/
 296 /***************************************************************/
 297
 298 #include "AliRun.h"
 299 #include "AliESD.h"
 300 #include "AliRunLoader.h"
 301 #include "AliTPCtrack.h"
 302 #include "TTree.h"
 303 #include "TBranch.h"
 304 #include "TH2D.h"
 305 #include "TCanvas.h"
 306 #include "AliMagF.h"
 307
 308
 309
 310
 311 void AliHBTTrackPoints::testesd(Int_t entr,const char* fname )
 312 {
 313   delete gAlice;
 314   gAlice = 0x0;
 315   AliRunLoader* rl = AliRunLoader::Open();
 316   rl->LoadgAlice();
 317
 318   Float_t mf = rl->GetAliRun()->Field()->SolenoidField();
 319
 320
 321   TFile* fFile = TFile::Open(fname);
 322
 323   if (fFile == 0x0)
 324    {
 325      printf("testesd: There is no suche a ESD file\n");
 326      return;
 327    }
 328   AliESD* esd = dynamic_cast<AliESD*>(fFile->Get("0"));
 329   AliESDtrack *t = esd->GetTrack(entr);
 330   if (t == 0x0)
 331    {
 332      ::Error("testesd","Can not get track %d",entr);
 333      return;
 334    }
 335
 336
 337   Int_t N = 170;
 338   AliHBTTrackPoints* tp = new AliHBTTrackPoints(N,t,mf,1.);
 339
 340   Float_t xmin = -250;
 341   Float_t xmax =  250;
 342
 343   Float_t ymin = -250;
 344   Float_t ymax =  250;
 345
 346   Float_t zmin = -250;
 347   Float_t zmax =  250;
 348
 349   TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
 350   TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
 351   TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
 352
 353   TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
 354   TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
 355   TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
 356
 357   hxyt->SetDirectory(0x0);
 358   hxy->SetDirectory(0x0);
 359   hxyTR->SetDirectory(0x0);
 360
 361   hxzt->SetDirectory(0x0);
 362   hxz->SetDirectory(0x0);
 363   hxzTR->SetDirectory(0x0);
 364
 365   Float_t x,y,z;
 366
 367   for (Int_t i = 0;i<N;i++)
 368    {
 369       Double_t r = 84.1+i;
 370       tp->PositionAt(i,x,y,z);
 371       hxy->Fill(x,y);
 372       hxz->Fill(x,z);
 373       printf("Rdemanded %f\n",r);
 374       printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
 375
 376    }
 377
 378   rl->LoadTrackRefs();
 379   TTree* treeTR = rl->TreeTR();
 380   TBranch* b = treeTR->GetBranch("TPC");
 381
 382   TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
 383   AliTrackReference* tref;
 384   b->SetAddress(&trackrefs);
 385
 386   Int_t tlab = TMath::Abs(t->GetLabel());
 387
 388   Int_t netr = (Int_t)treeTR->GetEntries();
 389   printf("Found %d entries in TR tree\n",netr);
 390
 391   for (Int_t e = 0; e < netr; e++)
 392    {
 393      treeTR->GetEntry(e);
 394      tref = (AliTrackReference*)trackrefs->At(0);
 395      if (tref == 0x0) continue;
 396      if (tref->GetTrack() != tlab) continue;
 397
 398      printf("Found %d entries in TR array\n",trackrefs->GetEntries());
 399
 400      for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
 401       {
 402         tref = (AliTrackReference*)trackrefs->At(i);
 403         if (tref->GetTrack() != tlab) continue;
 404         x = tref->X();
 405         y = tref->Y();
 406         z = tref->Z();
 407         printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
 408
 409         hxzTR->Fill(x,z);
 410         hxyTR->Fill(x,y);
 411         for (Int_t j = 1; j < 10; j++)
 412          {
 413             hxyTR->Fill(x, y+j*0.1);
 414             hxyTR->Fill(x, y-j*0.1);
 415             hxyTR->Fill(x+j*0.1,y);
 416             hxyTR->Fill(x-j*0.1,y);
 417
 418             hxzTR->Fill(x,z-j*0.1);
 419             hxzTR->Fill(x,z+j*0.1);
 420             hxzTR->Fill(x-j*0.1,z);
 421             hxzTR->Fill(x+j*0.1,z);
 422          }
 423       }
 424       break;
 425     }
 426   hxy->Draw("");
 427 //  hxzt->Draw("same");
 428   hxyTR->Draw("same");
 429
 430   delete rl;
 431 }
 432
 433 /***************************************************************/
 434 /***************************************************************/
 435 /***************************************************************/
 436
 437 void AliHBTTrackPoints::testtpc(Int_t entr)
 438 {
 439   delete gAlice;
 440   gAlice = 0x0;
 441   AliRunLoader* rl = AliRunLoader::Open();
 442   AliLoader* l = rl->GetLoader("TPCLoader");
 443   rl->LoadgAlice();
 444   AliKalmanTrack::SetConvConst(100/0.299792458/0.2/rl->GetAliRun()->Field()->Factor());
 445   l->LoadTracks();
 446   AliTPCtrack* t = new AliTPCtrack();
 447   TBranch* b=l->TreeT()->GetBranch("tracks");
 448   b->SetAddress(&t);
 449   l->TreeT()->GetEntry(entr);
 450   Int_t N = 160;
 451   AliHBTTrackPoints* tp = new AliHBTTrackPoints(N,t,1.);
 452
 453   Float_t xmin = -250;
 454   Float_t xmax =  250;
 455
 456   Float_t ymin = -250;
 457   Float_t ymax =  250;
 458
 459   Float_t zmin = -250;
 460   Float_t zmax =  250;
 461
 462   TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
 463   TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
 464   TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
 465
 466   TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
 467   TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
 468   TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
 469
 470   hxyt->SetDirectory(0x0);
 471   hxy->SetDirectory(0x0);
 472   hxyTR->SetDirectory(0x0);
 473
 474   hxzt->SetDirectory(0x0);
 475   hxz->SetDirectory(0x0);
 476   hxzTR->SetDirectory(0x0);
 477
 478   Float_t x,y,z;
 479
 480   for (Int_t i = 0;i<N;i++)
 481    {
 482       Double_t r = 84.1+i;
 483       tp->PositionAt(i,x,y,z);
 484       hxy->Fill(x,y);
 485       hxz->Fill(x,z);
 486       printf("Rdemanded %f\n",r);
 487       printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
 488
 489       //BUT they are local!!!!
 490       t->PropagateTo(r);
 491 //      Double_t phi = t->Phi();
 492       Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius
 493
 494       Double_t alpha = t->GetAlpha();
 495       Double_t salpha = TMath::Sin(alpha);
 496       Double_t calpha = TMath::Cos(alpha);
 497       x = t->GetX()*calpha - t->GetY()*salpha;
 498       y = t->GetX()*salpha + t->GetY()*calpha;
 499       z = t->GetZ();
 500
 501       printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
 502
 503       printf("tpz - tz %f\n",z-t->GetZ());
 504       printf("\n");
 505       hxyt->Fill(x,y);
 506       hxzt->Fill(x,z);
 507
 508    }
 509
 510   rl->LoadTrackRefs();
 511   TTree* treeTR = rl->TreeTR();
 512   b = treeTR->GetBranch("TPC");
 513
 514   TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
 515   AliTrackReference* tref;
 516   b->SetAddress(&trackrefs);
 517
 518   Int_t tlab = TMath::Abs(t->GetLabel());
 519
 520   Int_t netr = (Int_t)treeTR->GetEntries();
 521   printf("Found %d entries in TR tree\n",netr);
 522
 523   for (Int_t e = 0; e < netr; e++)
 524    {
 525      treeTR->GetEntry(e);
 526      tref = (AliTrackReference*)trackrefs->At(0);
 527      if (tref == 0x0) continue;
 528      if (tref->GetTrack() != tlab) continue;
 529
 530      printf("Found %d entries in TR array\n",trackrefs->GetEntries());
 531
 532      for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
 533       {
 534         tref = (AliTrackReference*)trackrefs->At(i);
 535         if (tref->GetTrack() != tlab) continue;
 536         x = tref->X();
 537         y = tref->Y();
 538         z = tref->Z();
 539         printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
 540
 541         hxzTR->Fill(x,z);
 542         hxyTR->Fill(x,y);
 543         for (Int_t j = 1; j < 10; j++)
 544          {
 545             hxyTR->Fill(x, y+j*0.1);
 546             hxyTR->Fill(x, y-j*0.1);
 547             hxyTR->Fill(x+j*0.1,y);
 548             hxyTR->Fill(x-j*0.1,y);
 549
 550             hxzTR->Fill(x,z-j*0.1);
 551             hxzTR->Fill(x,z+j*0.1);
 552             hxzTR->Fill(x-j*0.1,z);
 553             hxzTR->Fill(x+j*0.1,z);
 554          }
 555       }
 556       break;
 557     }
 558   hxz->Draw("");
 559 //  hxzt->Draw("same");
 560   hxzTR->Draw("same");
 561
 562   delete rl;
 563 }