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