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