AliRoot Core  v5-06-30 (35d6c57)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TestTPCTrackHits.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
21 
22 #include "alles.h"
23 #include "AliObjectArray.h"
24 #include "AliTPCTrackHits.h"
25 #include "AliArrayBranch.h"
26 #include "TestTPCTrackHits.h"
27 
31 
32 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd=0);
33 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
34 
35 void ConvertHits(const char * benchmark, Bool_t debug)
36 {
37 
40 
41  TFile f("galice.root","update");
42  TClonesArray *arr = new TClonesArray("AliTPChit",100);
43  TTree * treeCl =(TTree*)f.Get("TreeH0");
44  TBranch *branch = treeCl->GetBranch("TPC");
45  AliTPCTrackHits * myhits = new AliTPCTrackHits;
46  branch->SetAddress(&arr);
47  //particle
48  TTree * treeP = (TTree*)f.Get("TreeK0");
49  TBranch *branchP = treeP->GetBranch("Particles");
50  TClonesArray *arrp = new TClonesArray("TParticle",100);
51  branchP->SetAddress(&arrp);
52  branchP->GetEvent(0);
53  //
54  TFile f2("treeh.root","recreate");
55  f2.SetCompressionLevel(10);
56  TTree * treeh2 = new TTree("TreeTPCH","TreeTPCH");
57  //treeh2->AliBranch("TPC","AliTPCTrackHits", &myhits,4096);
58  treeh2->GetListOfBranches()->Add(new AliObjectBranch("TPC","AliTPCTrackHits",
59  &myhits,treeh2,4096,1));
60 
61  TFile f3("treehcl.root","recreate");
62  f3.SetCompressionLevel(2);
63  TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
64  treeh3->Branch("TPC", &arr,64000,100);
65 
66  gBenchmark->Start(benchmark);
67  Int_t trsize = (Int_t)treeCl->GetEntries();
68  for (Int_t i=0;i<300;i++){
69  arr->Clear();
70  Int_t size = branch->GetEvent(i);
71  //if (size>0){
72  if ((size>0)&& arr->GetEntriesFast()>0) {
73  f.cd();
74  myhits->Clear();
75  myhits =MakeTrack(arr,arrp,myhits);
76  CompareHits(arr,myhits,debug);
77  f2.cd();
78  treeh2->Fill();
79  f3.cd();
80  treeh3->Fill();
81  if ((i%10)==0) cout<<i<<"\n";
82  }
83  }
84  f3.cd();
85  treeh3->Write();
86  f2.cd();
87  treeh2->Write();
88  f2.Close();
89  f.Close();
90  delete myhits;
91  gBenchmark->Stop(benchmark);
92  gBenchmark->Show(benchmark);
93 }
94 
95 
96 
97 void CompareHits(const char * benchmark, Bool_t debug)
98 {
100 
101  TFile f2("treeh.root");
102  TTree * treeh2 = (TTree*)f2.Get("TreeTPCH");
103  AliTPCTrackHits * myhits = new AliTPCTrackHits ;
104  AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC");
105  mybranch->SetAddress(&myhits);
106 
107 
108  TClonesArray *arr = new TClonesArray("AliTPChit",100);
109  TFile f3("treehcl.root");
110  TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
111  TBranch *branch = treeh3->GetBranch("TPC");
112  branch->SetAddress(&arr);
113 
114  cout<<"Lets go!\n";
115  gBenchmark->Start(benchmark);
116  Int_t trsize = (Int_t)treeh3->GetEntries();
117  for (Int_t i=0;i<300;i++){
118  Int_t size = branch->GetEvent(i);
119  mybranch->GetEvent(i);
120  //if (arr){
121  if ((arr->GetEntriesFast()>0) && size>0) {
122  CompareHits(arr,myhits,debug);
123  if ((i%10)==0) cout<<i<<"\n";
124  }
125  }
126  gBenchmark->Stop(benchmark);
127  gBenchmark->Show(benchmark);
128 }
129 
130 
131 void CompareHitsG(const char * benchmark, Bool_t debug)
132 {
134 
135  TFile f2("galice.root");
136  TTree * treeh2 = (TTree*)f2.Get("TreeH0");
137  AliTPCTrackHits * myhits = new AliTPCTrackHits ;
138  AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC2");
139  mybranch->SetAddress(&myhits);
140 
141 
142  TClonesArray *arr = new TClonesArray("AliTPChit",100);
143  //TFile f3("treehcl.root");
144  //TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
145  TBranch *branch = treeh2->GetBranch("TPC");
146  branch->SetAddress(&arr);
147 
148  TFile f3("treehdelta.root","recreate");
149  f3.SetCompressionLevel(2);
150  TTree * treeh3 = new TTree("DelataH","DeltaH");
151  TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
152  treeh3->Branch("TPC", &arrd,64000,100);
153 
154  cout<<"Lets go!\n";
155  gBenchmark->Start(benchmark);
156  Int_t trsize = treeh2->GetEntries();
157  for (Int_t i=0;i<trsize;i++){
158  Int_t size = branch->GetEvent(i);
159  mybranch->GetEvent(i);
160  //if (arr){
161  if ((arr->GetEntriesFast()>0) && size>0) {
162  CompareHits(arr,myhits,debug,arrd);
163  }
164  if ((i%10)==0) cout<<i<<"\n";
165  treeh3->Fill();
166  }
167  treeh3->Write();
168  f3.Close();
169  gBenchmark->Stop(benchmark);
170  gBenchmark->Show(benchmark);
171 }
172 
173 
174 
175 
176 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)
177 {
179 
180  AliTPChit * hit;
181  // AliTPCTrackHits * myhits= new AliTPCTrackHits;
182  myhits->SetHitPrecision(0.002);
183  myhits->SetStepPrecision(0.003);
184  myhits->SetMaxDistance(100);
185  myhits->FlushHitStack(kTRUE);
186  for (Int_t i=0;i<arr->GetEntriesFast();i++){
187  hit = (AliTPChit*)arr->At(i);
188  if (hit){
189  TParticle *p = (TParticle*)arrp->At(hit->GetTrack());
190  Float_t momentum = TMath::Sqrt(p->Px()*p->Px()+p->Py()*p->Py());
191  Float_t ran= 100.*gRandom->Rndm();
192  if (ran<1.) myhits->FlushHitStack(kTRUE);
193  if (momentum<0.01) {//if small momentum - not precise
194  myhits->SetHitPrecision(0.05);
195  myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
196  hit->Z(), hit->fQ+1000);
197  }
198  else {
199  myhits->SetHitPrecision(0.002);
200  myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
201  hit->Z(), hit->fQ);
202  }
203  }
204  }
205  myhits->FlushHitStack();
206  return myhits;
207 }
208 
209 
210 
211 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
212 {
215 
216  AliTPChit * hit, *hit2;
217  if (arrd) arrd->Clear();
218 
219  for (Int_t i=0;i<arr->GetEntriesFast();i++){
220  hit = (AliTPChit*)arr->At(i);
221  if (hit){
222  if (i==0) myhits->First();
223  else myhits->Next();
224  hit2 = myhits->GetHit();
225  if (!hit2) {
226  hit2=0;
227  if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
228  }
229 
230  if (hit&&hit2){
231  // if (hit){
232 
233  AliTrackHitsParam *param= myhits->GetParam();
234  AliHitInfo * info = myhits->GetHitInfo();
235  //
236  if (arrd) {
237  TClonesArray &larrd = *arrd;
238  AliTPChitD * h = new(larrd[i]) AliTPChitD;
239  h->SetX(hit->X());
240  h->SetY(hit->Y());
241  h->SetZ(hit->Z());
242  h->SetTrack(hit->GetTrack());
243  h->fQ = hit->fQ;
244  h->fSector = hit->fSector;
245  AliTPChit * hitd = h->GetDelta();
246  hitd->SetX(hit->X()-hit2->X());
247  hitd->SetY(hit->Y()-hit2->Y());
248  hitd->SetZ(hit->Z()-hit2->Z());
249  hitd->SetTrack(hit->GetTrack()-hit2->GetTrack());
250  hitd->fQ = hit->fQ-hit2->fQ;
251  hitd->fSector = hit->fSector-hit2->fSector;
252  }
253 
254  if (debug){
255  Float_t dd =
256  TMath::Sqrt(
257  (hit->X()-hit2->X())*(hit->X()-hit2->X())+
258  (hit->Y()-hit2->Y())*(hit->Y()-hit2->Y())+
259  (hit->Z()-hit2->Z())*(hit->Z()-hit2->Z()));
260  printf("C1\t%d\t%d\t%d\t%d\t",
261  hit->fSector,hit2->fSector,hit->GetTrack(),hit2->GetTrack());
262  printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
263  hit->X(),hit2->X(), hit->Y(),hit2->Y(),hit->Z(),hit2->Z());
264  printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
265  dd, param->fR,param->fZ,param->fFi,param->fAn,
266  param->fAd,param->fTheta);
267  printf("%d\t%d\t%d\n",
268  (Int_t)info->fHitDistance, (Int_t)hit->fQ, (Int_t)hit2->fQ);
269  }
270  }
271  }
272  }
273 }
274 
275 
276 
277 
278 //Filter for paw visualisation of results
279 //
280 //sed filter
281  //sed '/*/d' dd.txt | sed '/C/d'| cat -n > out.txt
282  // sed -n '/C2/p' dd.txt | sed 's/C1//' | cat -n >out2.txt
283  // sed -n '/C3/p' dd.txt | sed 's/C2//' | cat -n >out3.txt
284 
285 //filter
286 //myfilter C1 dd.txt >out1.txt
287 //myfilter C2 dd.txt >out2.txt
288 //myfilter C3 dd.txt >out3.txt
289 // sed -n 1,50000p
290 // sed -n 1,50000p dd.txt | myfilter C1 | sed -n 1,50000p >out1.txt
291 
292 //myfilter C1 dd.txt | sed -n 1,50000p >out1.txt
293 //myfilter C2 dd.txt | sed -n 1,50000p >out2.txt
294 //myfilter C3 dd.txt | sed -n 1,50000p >out3.txt
295 /*
296 paw visualisation
297  Nt/Create 1 'Test of Match' 21 ! ! event v1 v2 t1 t2 x1 x2 y1 y2 z1 z2 dd r z fi an ad theta dist q1 q2
298  Nt/Read 1 out1.txt
299  nt
300  Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r z fi alfa theta dr1 ddr ddz ddrfi dd
301  Nt/Read 2 out2.txt
302 
303  Nt/Create 3 'Test of Match2' 2 ! ! event stacksize
304  Nt/Read 3 out3.txt
305 */
306 
307 void Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
308  Double_t fSumX, Double_t fSumX2, Double_t fSumX3,
309  Double_t fSumX4, Int_t n,
310  Double_t &a, Double_t &b, Double_t &c)
311 {
313 
314  Double_t det =
315  n* (fSumX2*fSumX4-fSumX3*fSumX3) -
316  fSumX* (fSumX*fSumX4-fSumX3*fSumX2)+
317  fSumX2* (fSumX*fSumX3-fSumX2*fSumX2);
318 
319  if (TMath::Abs(det)>0){
320  a =
321  (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
322  fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
323  fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
324  b=
325  (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
326  fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
327  fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
328  c=
329  (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
330  fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
331  fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
332  cout<<a<<"\t"<<b<<"\t"<<c<<"\n";
333  }
334 }
335 
336 void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
337 {
338  Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
339  fSumY = fSumYX = fSumYX2 = fSumX = fSumX2 = fSumX3 = fSumX4 =0;
340  Double_t a2,b2,c2;
341  for (Int_t i=0;i<n; i++){
342  Float_t x = gRandom->Rndm();
343  Float_t y = a+b*x+c*x*x;
344  fSumY+=y;
345  fSumYX+=y*x;
346  fSumYX2+=y*x*x;
347  fSumX += x;
348  fSumX2 += x*x;
349  fSumX3 += x*x*x;
350  fSumX4 += x*x*x*x;
351  }
352  Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);
353 }
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void CompareHits(TClonesArray *arr, AliTPCTrackHits *myhits, Bool_t debug, TClonesArray *arrd=0)
TBenchmark * gBenchmark
void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
TFile f("CalibObjects.root")
TTree * treeP
ClassImp(TPCGenInfo)
Definition: AliTPCCmpNG.C:254
Macro to compare TClonesArray hits with interpolated hits.
void CompareHitsG(const char *benchmark, Bool_t debug)
void Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2, Double_t fSumX, Double_t fSumX2, Double_t fSumX3, Double_t fSumX4, Int_t n, Double_t &a, Double_t &b, Double_t &c)
AliTPChit * GetDelta()
void ConvertHits(const char *benchmark, Bool_t debug)
AliTPCTrackHits * MakeTrack(TClonesArray *arr, TClonesArray *arrp, AliTPCTrackHits *myhits)
Int_t debug