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