Check on curvature not position
[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 #include "AliTrackReference.h"
16
17 #include "AliTPCtrack.h"
18 #include "AliESDtrack.h"
19
20 #include <TMath.h>
21 #include <TClonesArray.h>
22
23 ClassImp(AliHBTTrackPoints)
24
25 Int_t AliHBTTrackPoints::fgDebug = 0;
26 AliHBTTrackPoints::AliHBTTrackPoints():
27  fN(0),
28  fX(0x0),
29  fY(0x0),
30  fZ(0x0)
31 {
32   //constructor
33 }
34 /***************************************************************/
35
36 AliHBTTrackPoints::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    }
56
57   Double_t x;
58   Double_t par[5];
59   track->GetInnerExternalParameters(x,par);     //get properties of the track
60   if (par[4] == 0)
61    {
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;
70    }
71    
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);
77   
78 }
79 /***************************************************************/
80
81 AliHBTTrackPoints::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     
109   Double_t x = 0;
110   Double_t par[5];
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
118 void 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
127   
128   Double_t y = par[0];
129   Double_t z0 = par[1];
130       
131   Double_t phi0local = TMath::ATan2(y,x);
132   Double_t phi0global = phi0local + alpha;
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
143   if (GetDebug() > 9) 
144     Info("AliHBTTrackPoints","Radius0 %f, Real Radius %f",r0,r); 
145   
146   if (GetDebug() > 5) 
147     Info("AliHBTTrackPoints","Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local);
148     
149   Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
150   
151   //this calculattions are assuming new (current) model 
152   Double_t tmp = par[2];
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;
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
182      Double_t phi = phi0global + factorPhi - factorPhi0;
183      
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;
203      fZ[i] = z0 + factorZ - factorZ0;
204      fX[i] = rc*TMath::Cos(phi);
205      fY[i] = rc*TMath::Sin(phi);
206      
207      if ( GetDebug() > 2 )
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 }
214 /***************************************************************/
215
216 AliHBTTrackPoints::~AliHBTTrackPoints()
217 {
218   //destructor
219   delete [] fX;
220   delete [] fY;
221   delete [] fZ;
222 }
223 /***************************************************************/
224
225 void 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];
237   if ( GetDebug() > 1 )
238     {
239       Info("AliHBTTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
240     }
241 }
242 /***************************************************************/
243
244 Double_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    
258   Double_t sum = 0;
259   for (Int_t i = 0; i<fN; i++)
260    {
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       } 
282    }
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;
290 }
291 /***************************************************************/
292 /***************************************************************/
293 /***************************************************************/
294 /***************************************************************/
295 /***************************************************************/
296 /***************************************************************/
297
298 #include "AliRun.h"
299 #include "AliESD.h"
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
310
311 void AliHBTTrackPoints::testesd(Int_t entr,const char* fname )
312 {
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 /***************************************************************/
436
437 void AliHBTTrackPoints::testtpc(Int_t entr)
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   
453   Float_t xmin = -250;
454   Float_t xmax =  250;
455   
456   Float_t ymin = -250;
457   Float_t ymax =  250;
458
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);
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);
491 //      Double_t phi = t->Phi();
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
520   Int_t netr = (Int_t)treeTR->GetEntries();
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 }