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