1 #include "AliHBTTrackPoints.h"
2 //_________________________________
3 ////////////////////////////////////////////////////////////
5 // class AliHBTTrackPoints //
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) //
14 ////////////////////////////////////////////////////////////
15 #include "AliTrackReference.h"
17 #include "AliTPCtrack.h"
18 #include "AliESDtrack.h"
21 #include <TClonesArray.h>
23 ClassImp(AliHBTTrackPoints)
25 Int_t AliHBTTrackPoints::fgDebug = 0;
26 AliHBTTrackPoints::AliHBTTrackPoints():
34 /***************************************************************/
36 AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
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
48 Error("AliHBTTrackPoints","ESD track is null");
59 track->GetInnerExternalParameters(x,par); //get properties of the track
62 Error("AliHBTTrackPoints","This ESD track does not contain TPC information");
68 Error("AliHBTTrackPoints","Zero Magnetic field passed as parameter.");
72 Double_t alpha = track->GetInnerAlpha();
73 Double_t cc = 1000./0.299792458/mf;//conversion constant
76 MakePoints(dr,r0,x,par,c,alpha);
79 /***************************************************************/
81 AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
89 //dr - calculate points every dr cm, default every 30cm
92 Error("AliHBTTrackPoints","TPC track is null");
100 track->PropagateTo(r0);
102 //* This formation is now fixed in the following way: *
103 //* external param0: local Y-coordinate of a track (cm) *
104 //* external param1: local Z-coordinate of a track (cm) *
105 //* external param2: local sine of the track momentum azimuth angle *
106 //* external param3: tangent of the track momentum dip angle *
107 //* external param4: 1/pt (1/(GeV/c)) *
111 track->GetExternalParameters(x,par); //get properties of the track
113 Double_t alpha = track->GetAlpha();
114 Double_t c=track->GetC();
115 MakePoints(dr,r0,x,par,c,alpha);
118 void AliHBTTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
120 //Calculates points starting at radius r0
121 //spacing every dr (in radial direction)
122 // according to track parameters
123 // x - position in sector local reference frame. x i parallel to R and sector is symmetric with respect to x
124 // par - track external parameters; array with 5 elements; look at AliTPCtrack.h or AliESDtrack.h for their meaning
125 // c - track curvature
126 // alpha - sector's rotation angle (phi) == angle needed for local to global transformation
129 Double_t z0 = par[1];
131 Double_t phi0local = TMath::ATan2(y,x);
132 Double_t phi0global = phi0local + alpha;
134 if (phi0local<0) phi0local+=2*TMath::Pi();
135 if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
137 if (phi0global<0) phi0global+=2*TMath::Pi();
138 if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
140 Double_t r = TMath::Hypot(x,y);
144 Info("AliHBTTrackPoints","Radius0 %f, Real Radius %f",r0,r);
147 Info("AliHBTTrackPoints","Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local);
149 Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
151 //this calculattions are assuming new (current) model
152 Double_t tmp = par[2];
154 tmp = c*y + TMath::Sqrt(tmp);
155 Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
157 //Here is old model Cold=Cnew/2.
158 Double_t dcasq = dca*dca;
160 Double_t cst1 = (1.+c2*dca)*dca;//first constant
161 Double_t cst2 = 1. + 2.*c2*dca;//second constant
163 Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
164 Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
166 for(Int_t i = 0; i<fN; i++)
168 Double_t rc = r0 + i*dr;
169 Double_t factorPhi = TMath::ASin( (c2*rc + cst1/rc)/cst2 );//factor phi od rc
170 Double_t phi = phi0global + factorPhi - factorPhi0;
172 Double_t factorZ = TMath::ASin(c2*TMath::Sqrt((rc*rc-dcasq)/cst2))*par[3]/c2;
173 fZ[i] = z0 + factorZ - factorZ0;
174 fX[i] = rc*TMath::Cos(phi);
175 fY[i] = rc*TMath::Sin(phi);
177 if ( GetDebug() > 2 )
179 Info("AliHBTTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
180 fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i]));
184 /***************************************************************/
186 AliHBTTrackPoints::~AliHBTTrackPoints()
193 /***************************************************************/
195 void AliHBTTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
197 //returns position at point n
200 Error("PositionAt","Point %d out of range",n);
207 if ( GetDebug() > 1 )
209 Info("AliHBTTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
212 /***************************************************************/
214 Double_t AliHBTTrackPoints::AvarageDistance(const AliHBTTrackPoints& tr)
216 //returns the aritmethic avarage distance between two tracks
219 Error("AvarageDistance","Number of points is not equal");
222 if ( (fN <= 0) || (tr.fN <=0) )
224 Error("AvarageDistance","One of tracks is empty");
229 for (Int_t i = 0; i<fN; i++)
233 // Float_t r1sq = fX[i]*fX[i]+fY[i]*fY[i];
234 // Float_t r2sq = tr.fX[i]*tr.fX[i]+tr.fY[i]*tr.fY[i];
235 Float_t r1sq = TMath::Hypot(fX[i],fY[i]);
236 Float_t r2sq = TMath::Hypot(tr.fX[i],tr.fY[i]);
237 Info("AvarageDistance","radii: %f %f",r1sq,r2sq);
241 Double_t dx = fX[i]-tr.fX[i];
242 Double_t dy = fY[i]-tr.fY[i];
243 Double_t dz = fZ[i]-tr.fZ[i];
244 sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
248 Info("AvarageDistance","Diff: x ,y z: %f , %f, %f",dx,dy,dz);
249 Info("AvarageDistance","xxyyzz %f %f %f %f %f %f",
250 fX[i],tr.fX[i],fY[i],tr.fY[i],fZ[i],tr.fZ[i]);
254 Double_t retval = sum/((Double_t)fN);
257 Info("AvarageDistance","Avarage distance is %f.",retval);
261 /***************************************************************/
262 /***************************************************************/
263 /***************************************************************/
264 /***************************************************************/
265 /***************************************************************/
266 /***************************************************************/
270 #include "AliRunLoader.h"
271 #include "AliTPCtrack.h"
281 void AliHBTTrackPoints::testesd(Int_t entr,const char* fname )
285 AliRunLoader* rl = AliRunLoader::Open();
288 Float_t mf = rl->GetAliRun()->Field()->SolenoidField();
291 TFile* fFile = TFile::Open(fname);
295 printf("testesd: There is no suche a ESD file\n");
298 AliESD* esd = dynamic_cast<AliESD*>(fFile->Get("0"));
299 AliESDtrack *t = esd->GetTrack(entr);
302 ::Error("testesd","Can not get track %d",entr);
308 AliHBTTrackPoints* tp = new AliHBTTrackPoints(N,t,mf,1.);
319 TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
320 TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
321 TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
323 TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
324 TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
325 TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
327 hxyt->SetDirectory(0x0);
328 hxy->SetDirectory(0x0);
329 hxyTR->SetDirectory(0x0);
331 hxzt->SetDirectory(0x0);
332 hxz->SetDirectory(0x0);
333 hxzTR->SetDirectory(0x0);
337 for (Int_t i = 0;i<N;i++)
340 tp->PositionAt(i,x,y,z);
343 printf("Rdemanded %f\n",r);
344 printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
349 TTree* treeTR = rl->TreeTR();
350 TBranch* b = treeTR->GetBranch("TPC");
352 TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
353 AliTrackReference* tref;
354 b->SetAddress(&trackrefs);
356 Int_t tlab = TMath::Abs(t->GetLabel());
358 Int_t netr = (Int_t)treeTR->GetEntries();
359 printf("Found %d entries in TR tree\n",netr);
361 for (Int_t e = 0; e < netr; e++)
364 tref = (AliTrackReference*)trackrefs->At(0);
365 if (tref == 0x0) continue;
366 if (tref->GetTrack() != tlab) continue;
368 printf("Found %d entries in TR array\n",trackrefs->GetEntries());
370 for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
372 tref = (AliTrackReference*)trackrefs->At(i);
373 if (tref->GetTrack() != tlab) continue;
377 printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
381 for (Int_t j = 1; j < 10; j++)
383 hxyTR->Fill(x, y+j*0.1);
384 hxyTR->Fill(x, y-j*0.1);
385 hxyTR->Fill(x+j*0.1,y);
386 hxyTR->Fill(x-j*0.1,y);
388 hxzTR->Fill(x,z-j*0.1);
389 hxzTR->Fill(x,z+j*0.1);
390 hxzTR->Fill(x-j*0.1,z);
391 hxzTR->Fill(x+j*0.1,z);
397 // hxzt->Draw("same");
403 /***************************************************************/
404 /***************************************************************/
405 /***************************************************************/
407 void AliHBTTrackPoints::testtpc(Int_t entr)
411 AliRunLoader* rl = AliRunLoader::Open();
412 AliLoader* l = rl->GetLoader("TPCLoader");
414 AliKalmanTrack::SetConvConst(100/0.299792458/0.2/rl->GetAliRun()->Field()->Factor());
416 AliTPCtrack* t = new AliTPCtrack();
417 TBranch* b=l->TreeT()->GetBranch("tracks");
419 l->TreeT()->GetEntry(entr);
421 AliHBTTrackPoints* tp = new AliHBTTrackPoints(N,t,1.);
432 TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
433 TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
434 TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
436 TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
437 TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
438 TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
440 hxyt->SetDirectory(0x0);
441 hxy->SetDirectory(0x0);
442 hxyTR->SetDirectory(0x0);
444 hxzt->SetDirectory(0x0);
445 hxz->SetDirectory(0x0);
446 hxzTR->SetDirectory(0x0);
450 for (Int_t i = 0;i<N;i++)
453 tp->PositionAt(i,x,y,z);
456 printf("Rdemanded %f\n",r);
457 printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
459 //BUT they are local!!!!
461 // Double_t phi = t->Phi();
462 Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius
464 Double_t alpha = t->GetAlpha();
465 Double_t salpha = TMath::Sin(alpha);
466 Double_t calpha = TMath::Cos(alpha);
467 x = t->GetX()*calpha - t->GetY()*salpha;
468 y = t->GetX()*salpha + t->GetY()*calpha;
471 printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
473 printf("tpz - tz %f\n",z-t->GetZ());
481 TTree* treeTR = rl->TreeTR();
482 b = treeTR->GetBranch("TPC");
484 TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
485 AliTrackReference* tref;
486 b->SetAddress(&trackrefs);
488 Int_t tlab = TMath::Abs(t->GetLabel());
490 Int_t netr = (Int_t)treeTR->GetEntries();
491 printf("Found %d entries in TR tree\n",netr);
493 for (Int_t e = 0; e < netr; e++)
496 tref = (AliTrackReference*)trackrefs->At(0);
497 if (tref == 0x0) continue;
498 if (tref->GetTrack() != tlab) continue;
500 printf("Found %d entries in TR array\n",trackrefs->GetEntries());
502 for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
504 tref = (AliTrackReference*)trackrefs->At(i);
505 if (tref->GetTrack() != tlab) continue;
509 printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
513 for (Int_t j = 1; j < 10; j++)
515 hxyTR->Fill(x, y+j*0.1);
516 hxyTR->Fill(x, y-j*0.1);
517 hxyTR->Fill(x+j*0.1,y);
518 hxyTR->Fill(x-j*0.1,y);
520 hxzTR->Fill(x,z-j*0.1);
521 hxzTR->Fill(x,z+j*0.1);
522 hxzTR->Fill(x-j*0.1,z);
523 hxzTR->Fill(x+j*0.1,z);
529 // hxzt->Draw("same");