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