Track Points Added. Class needed for Anti-Merging cut
[u/mrichter/AliRoot.git] / HBTAN / AliHBTTrackPoints.cxx
1 #include "AliHBTTrackPoints.h"
2 #include "AliTPCtrack.h"
3 #include <TMath.h>
4 #include "AliTrackReference.h"
5 #include <TClonesArray.h>
6
7 ClassImp(AliHBTTrackPoints)
8
9 AliHBTTrackPoints::AliHBTTrackPoints():
10  fN(0),
11  fX(0x0),
12  fY(0x0),
13  fZ(0x0)
14 {
15   //constructor
16 }
17
18 AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
19  fN(n),
20  fX(new Float_t[fN]),
21  fY(new Float_t[fN]),
22  fZ(new Float_t[fN])
23 {
24   //constructor
25   //r0 starting radius
26   //dr - calculate points every dr cm, default every 30cm
27   if (track == 0x0)
28    {
29      Error("AliHBTTrackPoints","TPC track is null");
30      fN = 0;
31      delete [] fX;
32      delete [] fY;
33      delete [] fZ;
34      fX = fY = fZ = 0x0;
35      return;
36    }
37   track->PropagateTo(r0);
38
39  //* This formation is now fixed in the following way:                         *
40  //*      external param0:   local Y-coordinate of a track (cm)                *
41  //*      external param1:   local Z-coordinate of a track (cm)                *
42  //*      external param2:   local sine of the track momentum azimuth angle    *
43  //*      external param3:   tangent of the track momentum dip angle           *
44  //*      external param4:   1/pt (1/(GeV/c))                                  *
45     
46   Double_t xk;
47   Double_t par[5];
48   track->GetExternalParameters(xk,par);     //get properties of the track
49   
50   Double_t x = track->GetX();
51   Double_t y = track->GetY();
52   Double_t z0 = track->GetZ();
53       
54   Double_t phi0local = TMath::ATan2(y,x);
55   Double_t phi0global = phi0local + track->GetAlpha();
56   
57   if (phi0local<0) phi0local+=2*TMath::Pi();
58   if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
59   
60   if (phi0global<0) phi0global+=2*TMath::Pi();
61   if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
62   
63   Double_t r = TMath::Hypot(x,y);
64   
65
66   if (GetDebug()) 
67     Info("AliHBTTrackPoints","Radius0 %f, Real Radius %f",r0,r); 
68   
69   if (GetDebug()) 
70     Info("AliHBTTrackPoints","Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local);
71   Double_t c=track->GetC();
72   Double_t eta = track->GetEta();
73   
74   //this calculattions are assuming new (current) model 
75   Double_t tmp = c*x - eta;
76   tmp = 1. - tmp*tmp;
77   tmp = c*y + TMath::Sqrt(tmp);
78   Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
79
80   //Here is old model Cold=Cnew/2.
81   Double_t dcasq = dca*dca;
82   Double_t c2 = c/2.;
83   Double_t cst1 = (1.+c2*dca)*dca;//first constant
84   Double_t cst2 = 1. + 2.*c2*dca;//second constant
85       
86   Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
87   Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
88  
89   for(Int_t i = 0; i<fN; i++)
90    {
91      Double_t rc = r0 + i*dr;
92      Double_t factorPhi = TMath::ASin( (c2*rc + cst1/rc)/cst2 );//factor phi od rc
93      Double_t phi = phi0global + factorPhi - factorPhi0;
94      
95      Double_t factorZ = TMath::ASin(c2*TMath::Sqrt((rc*rc-dcasq)/cst2))*par[3]/c2;
96      fZ[i] = z0 + factorZ - factorZ0;
97      fX[i] = rc*TMath::Cos(phi);
98      fY[i] = rc*TMath::Sin(phi);
99      
100      if ( GetDebug() )
101       {
102         Info("AliHBTTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
103              fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i]));
104       }
105    }
106 }
107
108
109 AliHBTTrackPoints::~AliHBTTrackPoints()
110 {
111   //destructor
112   delete [] fX;
113   delete [] fY;
114   delete [] fZ;
115 }
116
117 void AliHBTTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
118 {
119   //returns position at point n
120   if ((n<0) || (n>fN))
121    {
122      Error("PositionAt","Point %d out of range",n);
123      return;
124    }
125   
126   x = fX[n];
127   y = fY[n];
128   z = fZ[n];
129   if ( GetDebug() )
130     {
131       Info("AliHBTTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
132     }
133
134 }
135
136
137 Double_t AliHBTTrackPoints::AvarageDistance(const AliHBTTrackPoints& tr)
138 {
139   //returns the aritmethic avarage distance between two tracks
140   if (fN != tr.fN)
141    {
142      Error("AvarageDistance","Number of points is not equal");
143      return -1;
144    }
145   if ( (fN <= 0) || (tr.fN <=0) )
146    {
147      Error("AvarageDistance","One of tracks is empty");
148      return -1;
149    }
150    
151   Double_t sum;
152   for (Int_t i = 0; i<fN; i++)
153    {
154       if (fR[i] != tr.fR[i])
155        {
156          Error("AvarageDistance","Different radii");
157          return -1;
158        }
159       Double_t dx = fX-tr.fX;
160       Double_t dy = fY-tr.fY;
161       Double_t dz = fZ-tr.fZ;
162       sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
163    }
164   return sum/((Double_t)fN);
165 }
166
167 #include "AliRun.h"
168 #include "AliRunLoader.h"
169 #include "AliTPCtrack.h"
170 #include "TTree.h"
171 #include "TBranch.h"
172 #include "TH2D.h"
173 #include "TCanvas.h"
174 #include "AliMagF.h"
175
176
177
178 void AliHBTTrackPoints::tp(Int_t entr)
179 {
180   delete gAlice;
181   gAlice = 0x0;
182   AliRunLoader* rl = AliRunLoader::Open();
183   AliLoader* l = rl->GetLoader("TPCLoader");
184   rl->LoadgAlice();
185   AliKalmanTrack::SetConvConst(100/0.299792458/0.2/rl->GetAliRun()->Field()->Factor());
186   l->LoadTracks();
187   AliTPCtrack* t = new AliTPCtrack();
188   TBranch* b=l->TreeT()->GetBranch("tracks");
189   b->SetAddress(&t);
190   l->TreeT()->GetEntry(entr);  
191   Int_t N = 160;
192   AliHBTTrackPoints* tp = new AliHBTTrackPoints(N,t,1.);
193   
194   TH2D* hxy = new TH2D("hxy","hxy",1000,150,250,1000,150,250);
195   TH2D* hxyt = new TH2D("hxyt","hxyt",1000,150,250,1000,150,250);
196   TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,150,250,1000,150,250);
197
198   TH2D* hxz = new TH2D("hxz","hxz",1000,150,250,1000,151,251);
199   TH2D* hxzt = new TH2D("hxzt","hxzt",1000,150,250,1000,151,251);
200   TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,150,250,1000,151,251);
201
202   hxyt->SetDirectory(0x0);   
203   hxy->SetDirectory(0x0);   
204   hxyTR->SetDirectory(0x0);   
205   
206   hxzt->SetDirectory(0x0);   
207   hxz->SetDirectory(0x0);   
208   hxzTR->SetDirectory(0x0);   
209
210   Float_t x,y,z;
211    
212   for (Int_t i = 0;i<N;i++)
213    {  
214       Double_t r = 84.1+i;
215       tp->PositionAt(i,x,y,z);
216       hxy->Fill(x,y);
217       hxz->Fill(x,z);
218       printf("Rdemanded %f\n",r);
219       printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
220       
221       //BUT they are local!!!!
222       t->PropagateTo(r);
223       Double_t phi = t->Phi();
224       Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius 
225       
226       Double_t alpha = t->GetAlpha();
227       Double_t salpha = TMath::Sin(alpha);
228       Double_t calpha = TMath::Cos(alpha);
229       x = t->GetX()*calpha - t->GetY()*salpha;
230       y = t->GetX()*salpha + t->GetY()*calpha;
231       z = t->GetZ();
232       
233       printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
234       
235       printf("tpz - tz %f\n",z-t->GetZ());
236       printf("\n");
237       hxyt->Fill(x,y);
238       hxzt->Fill(x,z);
239       
240    }
241
242   rl->LoadTrackRefs();
243   TTree* treeTR = rl->TreeTR();
244   b = treeTR->GetBranch("TPC");
245   
246   TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
247   AliTrackReference* tref;
248   b->SetAddress(&trackrefs);
249   
250   Int_t tlab = TMath::Abs(t->GetLabel());
251
252   Int_t netr = treeTR->GetEntries();
253   printf("Found %d entries in TR tree\n",netr);
254   
255   for (Int_t e = 0; e < netr; e++)
256    { 
257      treeTR->GetEntry(e);
258      tref = (AliTrackReference*)trackrefs->At(0);
259      if (tref == 0x0) continue;
260      if (tref->GetTrack() != tlab) continue;
261     
262      printf("Found %d entries in TR array\n",trackrefs->GetEntries());
263      
264      for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
265       {
266         tref = (AliTrackReference*)trackrefs->At(i);
267         if (tref->GetTrack() != tlab) continue;
268         x = tref->X();
269         y = tref->Y();
270         z = tref->Z();
271         printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
272
273         hxzTR->Fill(x,z);
274         hxyTR->Fill(x,y);
275         for (Int_t j = 1; j < 10; j++)
276          {
277             hxyTR->Fill(x, y+j*0.1);
278             hxyTR->Fill(x, y-j*0.1);
279             hxyTR->Fill(x+j*0.1,y);
280             hxyTR->Fill(x-j*0.1,y);
281
282             hxzTR->Fill(x,z-j*0.1);
283             hxzTR->Fill(x,z+j*0.1);
284             hxzTR->Fill(x-j*0.1,z);
285             hxzTR->Fill(x+j*0.1,z);
286          }
287       }
288       break; 
289     }  
290   hxz->Draw("");
291 //  hxzt->Draw("same");
292   hxzTR->Draw("same");
293   
294   delete rl;
295 }