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