HBTAN/AliHBTTrackPoints.cxx

   1 #include "AliHBTTrackPoints.h"
   2 #include "AliTPCtrack.h"
   3 #include <TMath.h>
   4 #include "AliTrackReference.h"
   5 #include <TClonesArray.h>
   6
   7 ClassImp(AliHBTTrackPoints)
   8
   9 AliHBTTrackPoints::AliHBTTrackPoints():
  10  fN(0),
  11  fX(0x0),
  12  fY(0x0),
  13  fZ(0x0)
  14 {
  15   //constructor
  16 }
  17
  18 AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
  19  fN(n),
  20  fX(new Float_t[fN]),
  21  fY(new Float_t[fN]),
  22  fZ(new Float_t[fN])
  23 {
  24   //constructor
  25   //r0 starting radius
  26   //dr - calculate points every dr cm, default every 30cm
  27   if (track == 0x0)
  28    {
  29      Error("AliHBTTrackPoints","TPC track is null");
  30      fN = 0;
  31      delete [] fX;
  32      delete [] fY;
  33      delete [] fZ;
  34      fX = fY = fZ = 0x0;
  35      return;
  36    }
  37   track->PropagateTo(r0);
  38
  39  //* This formation is now fixed in the following way:                         *
  40  //*      external param0:   local Y-coordinate of a track (cm)                *
  41  //*      external param1:   local Z-coordinate of a track (cm)                *
  42  //*      external param2:   local sine of the track momentum azimuth angle    *
  43  //*      external param3:   tangent of the track momentum dip angle           *
  44  //*      external param4:   1/pt (1/(GeV/c))                                  *
  45
  46   Double_t xk;
  47   Double_t par[5];
  48   track->GetExternalParameters(xk,par);     //get properties of the track
  49
  50   Double_t x = track->GetX();
  51   Double_t y = track->GetY();
  52   Double_t z0 = track->GetZ();
  53
  54   Double_t phi0local = TMath::ATan2(y,x);
  55   Double_t phi0global = phi0local + track->GetAlpha();
  56
  57   if (phi0local<0) phi0local+=2*TMath::Pi();
  58   if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
  59
  60   if (phi0global<0) phi0global+=2*TMath::Pi();
  61   if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
  62
  63   Double_t r = TMath::Hypot(x,y);
  64
  65
  66   if (GetDebug())
  67     Info("AliHBTTrackPoints","Radius0 %f, Real Radius %f",r0,r);
  68
  69   if (GetDebug())
  70     Info("AliHBTTrackPoints","Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local);
  71   Double_t c=track->GetC();
  72   Double_t eta = track->GetEta();
  73
  74   //this calculattions are assuming new (current) model
  75   Double_t tmp = c*x - eta;
  76   tmp = 1. - tmp*tmp;
  77   tmp = c*y + TMath::Sqrt(tmp);
  78   Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
  79
  80   //Here is old model Cold=Cnew/2.
  81   Double_t dcasq = dca*dca;
  82   Double_t c2 = c/2.;
  83   Double_t cst1 = (1.+c2*dca)*dca;//first constant
  84   Double_t cst2 = 1. + 2.*c2*dca;//second constant
  85
  86   Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
  87   Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
  88
  89   for(Int_t i = 0; i<fN; i++)
  90    {
  91      Double_t rc = r0 + i*dr;
  92      Double_t factorPhi = TMath::ASin( (c2*rc + cst1/rc)/cst2 );//factor phi od rc
  93      Double_t phi = phi0global + factorPhi - factorPhi0;
  94
  95      Double_t factorZ = TMath::ASin(c2*TMath::Sqrt((rc*rc-dcasq)/cst2))*par[3]/c2;
  96      fZ[i] = z0 + factorZ - factorZ0;
  97      fX[i] = rc*TMath::Cos(phi);
  98      fY[i] = rc*TMath::Sin(phi);
  99
 100      if ( GetDebug() )
 101       {
 102         Info("AliHBTTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
 103              fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i]));
 104       }
 105    }
 106 }
 107
 108
 109 AliHBTTrackPoints::~AliHBTTrackPoints()
 110 {
 111   //destructor
 112   delete [] fX;
 113   delete [] fY;
 114   delete [] fZ;
 115 }
 116
 117 void AliHBTTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
 118 {
 119   //returns position at point n
 120   if ((n<0) || (n>fN))
 121    {
 122      Error("PositionAt","Point %d out of range",n);
 123      return;
 124    }
 125
 126   x = fX[n];
 127   y = fY[n];
 128   z = fZ[n];
 129   if ( GetDebug() )
 130     {
 131       Info("AliHBTTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
 132     }
 133
 134 }
 135
 136
 137 Double_t AliHBTTrackPoints::AvarageDistance(const AliHBTTrackPoints& tr)
 138 {
 139   //returns the aritmethic avarage distance between two tracks
 140   if (fN != tr.fN)
 141    {
 142      Error("AvarageDistance","Number of points is not equal");
 143      return -1;
 144    }
 145   if ( (fN <= 0) || (tr.fN <=0) )
 146    {
 147      Error("AvarageDistance","One of tracks is empty");
 148      return -1;
 149    }
 150
 151   Double_t sum;
 152   for (Int_t i = 0; i<fN; i++)
 153    {
 154       if (fR[i] != tr.fR[i])
 155        {
 156          Error("AvarageDistance","Different radii");
 157          return -1;
 158        }
 159       Double_t dx = fX-tr.fX;
 160       Double_t dy = fY-tr.fY;
 161       Double_t dz = fZ-tr.fZ;
 162       sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 163    }
 164   return sum/((Double_t)fN);
 165 }
 166
 167 #include "AliRun.h"
 168 #include "AliRunLoader.h"
 169 #include "AliTPCtrack.h"
 170 #include "TTree.h"
 171 #include "TBranch.h"
 172 #include "TH2D.h"
 173 #include "TCanvas.h"
 174 #include "AliMagF.h"
 175
 176
 177
 178 void AliHBTTrackPoints::tp(Int_t entr)
 179 {
 180   delete gAlice;
 181   gAlice = 0x0;
 182   AliRunLoader* rl = AliRunLoader::Open();
 183   AliLoader* l = rl->GetLoader("TPCLoader");
 184   rl->LoadgAlice();
 185   AliKalmanTrack::SetConvConst(100/0.299792458/0.2/rl->GetAliRun()->Field()->Factor());
 186   l->LoadTracks();
 187   AliTPCtrack* t = new AliTPCtrack();
 188   TBranch* b=l->TreeT()->GetBranch("tracks");
 189   b->SetAddress(&t);
 190   l->TreeT()->GetEntry(entr);
 191   Int_t N = 160;
 192   AliHBTTrackPoints* tp = new AliHBTTrackPoints(N,t,1.);
 193
 194   TH2D* hxy = new TH2D("hxy","hxy",1000,150,250,1000,150,250);
 195   TH2D* hxyt = new TH2D("hxyt","hxyt",1000,150,250,1000,150,250);
 196   TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,150,250,1000,150,250);
 197
 198   TH2D* hxz = new TH2D("hxz","hxz",1000,150,250,1000,151,251);
 199   TH2D* hxzt = new TH2D("hxzt","hxzt",1000,150,250,1000,151,251);
 200   TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,150,250,1000,151,251);
 201
 202   hxyt->SetDirectory(0x0);
 203   hxy->SetDirectory(0x0);
 204   hxyTR->SetDirectory(0x0);
 205
 206   hxzt->SetDirectory(0x0);
 207   hxz->SetDirectory(0x0);
 208   hxzTR->SetDirectory(0x0);
 209
 210   Float_t x,y,z;
 211
 212   for (Int_t i = 0;i<N;i++)
 213    {
 214       Double_t r = 84.1+i;
 215       tp->PositionAt(i,x,y,z);
 216       hxy->Fill(x,y);
 217       hxz->Fill(x,z);
 218       printf("Rdemanded %f\n",r);
 219       printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
 220
 221       //BUT they are local!!!!
 222       t->PropagateTo(r);
 223       Double_t phi = t->Phi();
 224       Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius
 225
 226       Double_t alpha = t->GetAlpha();
 227       Double_t salpha = TMath::Sin(alpha);
 228       Double_t calpha = TMath::Cos(alpha);
 229       x = t->GetX()*calpha - t->GetY()*salpha;
 230       y = t->GetX()*salpha + t->GetY()*calpha;
 231       z = t->GetZ();
 232
 233       printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
 234
 235       printf("tpz - tz %f\n",z-t->GetZ());
 236       printf("\n");
 237       hxyt->Fill(x,y);
 238       hxzt->Fill(x,z);
 239
 240    }
 241
 242   rl->LoadTrackRefs();
 243   TTree* treeTR = rl->TreeTR();
 244   b = treeTR->GetBranch("TPC");
 245
 246   TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
 247   AliTrackReference* tref;
 248   b->SetAddress(&trackrefs);
 249
 250   Int_t tlab = TMath::Abs(t->GetLabel());
 251
 252   Int_t netr = treeTR->GetEntries();
 253   printf("Found %d entries in TR tree\n",netr);
 254
 255   for (Int_t e = 0; e < netr; e++)
 256    {
 257      treeTR->GetEntry(e);
 258      tref = (AliTrackReference*)trackrefs->At(0);
 259      if (tref == 0x0) continue;
 260      if (tref->GetTrack() != tlab) continue;
 261
 262      printf("Found %d entries in TR array\n",trackrefs->GetEntries());
 263
 264      for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
 265       {
 266         tref = (AliTrackReference*)trackrefs->At(i);
 267         if (tref->GetTrack() != tlab) continue;
 268         x = tref->X();
 269         y = tref->Y();
 270         z = tref->Z();
 271         printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
 272
 273         hxzTR->Fill(x,z);
 274         hxyTR->Fill(x,y);
 275         for (Int_t j = 1; j < 10; j++)
 276          {
 277             hxyTR->Fill(x, y+j*0.1);
 278             hxyTR->Fill(x, y-j*0.1);
 279             hxyTR->Fill(x+j*0.1,y);
 280             hxyTR->Fill(x-j*0.1,y);
 281
 282             hxzTR->Fill(x,z-j*0.1);
 283             hxzTR->Fill(x,z+j*0.1);
 284             hxzTR->Fill(x-j*0.1,z);
 285             hxzTR->Fill(x+j*0.1,z);
 286          }
 287       }
 288       break;
 289     }
 290   hxz->Draw("");
 291 //  hxzt->Draw("same");
 292   hxzTR->Draw("same");
 293
 294   delete rl;
 295 }