Fixing Coding violations
[u/mrichter/AliRoot.git] / HBTAN / AliHBTTrackPoints.cxx
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 ////////////////////////////////////////////////////////////
14
15 #include <TClonesArray.h>
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"
23
24 ClassImp(AliHBTTrackPoints)
25
26 Int_t AliHBTTrackPoints::fgDebug = 0;
27 AliHBTTrackPoints::AliHBTTrackPoints():
28  fN(0),
29  fX(0x0),
30  fY(0x0),
31  fZ(0x0)
32 {
33   //constructor
34 }
35 /***************************************************************/
36
37 AliHBTTrackPoints::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    }
57
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   
74   Double_t x;
75   Double_t par[5];
76   track->GetInnerExternalParameters(x,par);     //get properties of the track
77   if (par[4] == 0)
78    {
79      Error("AliHBTTrackPoints","This ESD track seem not to contain TPC information (curv is 0)");
80      return;
81    }
82
83   if (mf == 0.0)
84    {
85      Error("AliHBTTrackPoints","Zero Magnetic field passed as parameter.");
86      return;
87    }
88    
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);
94   
95 }
96 /***************************************************************/
97
98 AliHBTTrackPoints::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     
126   Double_t x = 0;
127   Double_t par[5];
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
135 void 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
144   
145   Double_t y = par[0];
146   Double_t z0 = par[1];
147       
148   Double_t phi0local = TMath::ATan2(y,x);
149   Double_t phi0global = phi0local + alpha;
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
160   if (GetDebug() > 9) 
161     Info("AliHBTTrackPoints","Radius0 %f, Real Radius %f",r0,r); 
162   
163   if (GetDebug() > 5) 
164     Info("AliHBTTrackPoints","Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local);
165     
166   Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
167   
168   //this calculattions are assuming new (current) model 
169   Double_t tmp = par[2];
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;
186      Double_t ftmp = (c2*rc + cst1/rc)/cst2;
187      if (ftmp > 1.0) 
188       {
189         if (GetDebug() > 1) 
190           Warning("AliHBTTrackPoints","ASin argument > 1 %f:",ftmp);
191         ftmp=1.0;
192       }
193      else if (ftmp < -1.0) 
194       {
195         if (GetDebug() > 1) 
196           Warning("AliHBTTrackPoints","ASin argument < -1 %f:",ftmp);
197         ftmp=-1.0;
198       }
199       
200      Double_t factorPhi = TMath::ASin( ftmp );//factor phi od rc
201      Double_t phi = phi0global + factorPhi - factorPhi0;
202      
203      ftmp = (rc*rc-dcasq)/cst2;
204      if (ftmp < 0.0) 
205       {
206         if (GetDebug() > 1) 
207           Warning("AliHBTTrackPoints","Sqrt argument < 0: %f",ftmp);
208         ftmp=0.0;
209       }
210      
211      ftmp = c2*TMath::Sqrt(ftmp);
212      if (ftmp > 1.0) 
213       {
214         if (GetDebug() > 1) 
215           Warning("AliHBTTrackPoints","ASin argument > 1: %f",ftmp);
216         ftmp=1.0;
217       }
218      else if (ftmp < -1.0) 
219       {
220         if (GetDebug() > 1) 
221           Warning("AliHBTTrackPoints","ASin argument < -1: %f",ftmp);
222         ftmp=-1.0;
223       }
224      Double_t factorZ = TMath::ASin(ftmp)*par[3]/c2;
225      fZ[i] = z0 + factorZ - factorZ0;
226      fX[i] = rc*TMath::Cos(phi);
227      fY[i] = rc*TMath::Sin(phi);
228      
229      if ( GetDebug() > 2 )
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 }
236 /***************************************************************/
237
238 AliHBTTrackPoints::~AliHBTTrackPoints()
239 {
240   //destructor
241   delete [] fX;
242   delete [] fY;
243   delete [] fZ;
244 }
245 /***************************************************************/
246
247 void 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];
259   if ( GetDebug() > 1 )
260     {
261       Info("AliHBTTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
262     }
263 }
264 /***************************************************************/
265
266 Double_t AliHBTTrackPoints::AvarageDistance(const AliHBTTrackPoints& tr)
267 {
268   //returns the aritmethic avarage distance between two tracks
269 //  Info("AvarageDistance","Entered");
270   if ( (fN <= 0) || (tr.fN <=0) )
271    {
272      if (GetDebug()) Warning("AvarageDistance","One of tracks is empty");
273      return -1;
274    }
275
276   if (fN != tr.fN)
277    {
278      Warning("AvarageDistance","Number of points is not equal");
279      return -1;
280    }
281    
282   Double_t sum = 0;
283   for (Int_t i = 0; i<fN; i++)
284    {
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      
300      if (GetDebug()>1)
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       } 
306    }
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;
314 }
315 /***************************************************************/
316 /***************************************************************/
317 /***************************************************************/
318 /***************************************************************/
319 /***************************************************************/
320 /***************************************************************/
321
322 #include "AliRun.h"
323 #include "AliESD.h"
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
334
335 void AliHBTTrackPoints::testesd(Int_t entr,const char* fname )
336 {
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 /***************************************************************/
460
461 void AliHBTTrackPoints::testtpc(Int_t entr)
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   
477   Float_t xmin = -250;
478   Float_t xmax =  250;
479   
480   Float_t ymin = -250;
481   Float_t ymax =  250;
482
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);
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);
515 //      Double_t phi = t->Phi();
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
544   Int_t netr = (Int_t)treeTR->GetEntries();
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 }