HBTAN/AliHBTTrackPoints.cxx

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