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