doxy: all TPC/Attic converted
[u/mrichter/AliRoot.git] / TPC / TestTPCTrackHits.cxx
CommitLineData
b6895dd8 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
0fd1fc0b 16/// Macro to compare TClonesArray hits with interpolated hits
17/// ConvertHits1 read
18///
19/// \author MI
cf98c13f 20
21#include "alles.h"
22#include "AliObjectArray.h"
23#include "AliTPCTrackHits.h"
24#include "AliArrayBranch.h"
25#include "TestTPCTrackHits.h"
26
0fd1fc0b 27/// \cond CLASSIMP
cf98c13f 28ClassImp(AliTPChitD)
0fd1fc0b 29/// \endcond
cf98c13f 30
31void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd=0);
32AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
33
34void ConvertHits(const char * benchmark, Bool_t debug)
35{
cf98c13f 36
0fd1fc0b 37 /// convert - not parametrised hits stored in Clones array
38 /// to
cf98c13f 39
40 TFile f("galice.root","update");
41 TClonesArray *arr = new TClonesArray("AliTPChit",100);
42 TTree * treeCl =(TTree*)f.Get("TreeH0");
43 TBranch *branch = treeCl->GetBranch("TPC");
44 AliTPCTrackHits * myhits = new AliTPCTrackHits;
45 branch->SetAddress(&arr);
46 //particle
47 TTree * treeP = (TTree*)f.Get("TreeK0");
48 TBranch *branchP = treeP->GetBranch("Particles");
49 TClonesArray *arrp = new TClonesArray("TParticle",100);
50 branchP->SetAddress(&arrp);
51 branchP->GetEvent(0);
52 //
53 TFile f2("treeh.root","recreate");
54 f2.SetCompressionLevel(10);
55 TTree * treeh2 = new TTree("TreeTPCH","TreeTPCH");
56 //treeh2->AliBranch("TPC","AliTPCTrackHits", &myhits,4096);
57 treeh2->GetListOfBranches()->Add(new AliObjectBranch("TPC","AliTPCTrackHits",
58 &myhits,treeh2,4096,1));
59
60 TFile f3("treehcl.root","recreate");
61 f3.SetCompressionLevel(2);
62 TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
d0f40f23 63 treeh3->Branch("TPC", &arr,64000,100);
cf98c13f 64
65 gBenchmark->Start(benchmark);
66 Int_t trsize = (Int_t)treeCl->GetEntries();
67 for (Int_t i=0;i<300;i++){
68 arr->Clear();
69 Int_t size = branch->GetEvent(i);
70 //if (size>0){
71 if ((size>0)&& arr->GetEntriesFast()>0) {
72 f.cd();
73 myhits->Clear();
74 myhits =MakeTrack(arr,arrp,myhits);
75 CompareHits(arr,myhits,debug);
76 f2.cd();
77 treeh2->Fill();
78 f3.cd();
79 treeh3->Fill();
80 if ((i%10)==0) cout<<i<<"\n";
81 }
82 }
83 f3.cd();
84 treeh3->Write();
85 f2.cd();
86 treeh2->Write();
87 f2.Close();
88 f.Close();
89 delete myhits;
90 gBenchmark->Stop(benchmark);
91 gBenchmark->Show(benchmark);
92}
93
94
95
96void CompareHits(const char * benchmark, Bool_t debug)
97{
98 //compare persistent hits
99 TFile f2("treeh.root");
100 TTree * treeh2 = (TTree*)f2.Get("TreeTPCH");
101 AliTPCTrackHits * myhits = new AliTPCTrackHits ;
102 AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC");
103 mybranch->SetAddress(&myhits);
104
105
106 TClonesArray *arr = new TClonesArray("AliTPChit",100);
107 TFile f3("treehcl.root");
108 TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
109 TBranch *branch = treeh3->GetBranch("TPC");
110 branch->SetAddress(&arr);
111
112 cout<<"Lets go!\n";
113 gBenchmark->Start(benchmark);
114 Int_t trsize = (Int_t)treeh3->GetEntries();
115 for (Int_t i=0;i<300;i++){
116 Int_t size = branch->GetEvent(i);
117 mybranch->GetEvent(i);
118 //if (arr){
119 if ((arr->GetEntriesFast()>0) && size>0) {
120 CompareHits(arr,myhits,debug);
121 if ((i%10)==0) cout<<i<<"\n";
122 }
123 }
124 gBenchmark->Stop(benchmark);
125 gBenchmark->Show(benchmark);
126}
127
128
129void CompareHitsG(const char * benchmark, Bool_t debug)
130{
131 //compare persistent hits
132 TFile f2("galice.root");
133 TTree * treeh2 = (TTree*)f2.Get("TreeH0");
134 AliTPCTrackHits * myhits = new AliTPCTrackHits ;
135 AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC2");
136 mybranch->SetAddress(&myhits);
137
138
139 TClonesArray *arr = new TClonesArray("AliTPChit",100);
140 //TFile f3("treehcl.root");
141 //TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
142 TBranch *branch = treeh2->GetBranch("TPC");
143 branch->SetAddress(&arr);
144
145 TFile f3("treehdelta.root","recreate");
146 f3.SetCompressionLevel(2);
147 TTree * treeh3 = new TTree("DelataH","DeltaH");
148 TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
d0f40f23 149 treeh3->Branch("TPC", &arrd,64000,100);
cf98c13f 150
151 cout<<"Lets go!\n";
152 gBenchmark->Start(benchmark);
153 Int_t trsize = treeh2->GetEntries();
154 for (Int_t i=0;i<trsize;i++){
155 Int_t size = branch->GetEvent(i);
156 mybranch->GetEvent(i);
157 //if (arr){
158 if ((arr->GetEntriesFast()>0) && size>0) {
159 CompareHits(arr,myhits,debug,arrd);
160 }
161 if ((i%10)==0) cout<<i<<"\n";
162 treeh3->Fill();
163 }
164 treeh3->Write();
165 f3.Close();
166 gBenchmark->Stop(benchmark);
167 gBenchmark->Show(benchmark);
168}
169
170
171
172
173AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)
174{
175 //
176 //make track wit hits according
177 AliTPChit * hit;
178 // AliTPCTrackHits * myhits= new AliTPCTrackHits;
179 myhits->SetHitPrecision(0.002);
180 myhits->SetStepPrecision(0.003);
181 myhits->SetMaxDistance(100);
182 myhits->FlushHitStack(kTRUE);
183 for (Int_t i=0;i<arr->GetEntriesFast();i++){
184 hit = (AliTPChit*)arr->At(i);
185 if (hit){
186 TParticle *p = (TParticle*)arrp->At(hit->GetTrack());
187 Float_t momentum = TMath::Sqrt(p->Px()*p->Px()+p->Py()*p->Py());
188 Float_t ran= 100.*gRandom->Rndm();
189 if (ran<1.) myhits->FlushHitStack(kTRUE);
190 if (momentum<0.01) {//if small momentum - not precise
191 myhits->SetHitPrecision(0.05);
192 myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
193 hit->Z(), hit->fQ+1000);
194 }
195 else {
196 myhits->SetHitPrecision(0.002);
197 myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
198 hit->Z(), hit->fQ);
199 }
200 }
201 }
202 myhits->FlushHitStack();
203 return myhits;
204}
205
206
207
208void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
209{
210 //
211 // if debug option kTRUE
212 // compare hits and write result to the stdoutput
213 AliTPChit * hit, *hit2;
214 if (arrd) arrd->Clear();
215
216 for (Int_t i=0;i<arr->GetEntriesFast();i++){
217 hit = (AliTPChit*)arr->At(i);
218 if (hit){
219 if (i==0) myhits->First();
220 else myhits->Next();
221 hit2 = myhits->GetHit();
222 if (!hit2) {
223 hit2=0;
224 if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
225 }
226
227 if (hit&&hit2){
228 // if (hit){
229
230 AliTrackHitsParam *param= myhits->GetParam();
231 AliHitInfo * info = myhits->GetHitInfo();
232 //
233 if (arrd) {
234 TClonesArray &larrd = *arrd;
235 AliTPChitD * h = new(larrd[i]) AliTPChitD;
236 h->SetX(hit->X());
237 h->SetY(hit->Y());
238 h->SetZ(hit->Z());
239 h->SetTrack(hit->GetTrack());
240 h->fQ = hit->fQ;
241 h->fSector = hit->fSector;
242 AliTPChit * hitd = h->GetDelta();
243 hitd->SetX(hit->X()-hit2->X());
244 hitd->SetY(hit->Y()-hit2->Y());
245 hitd->SetZ(hit->Z()-hit2->Z());
246 hitd->SetTrack(hit->GetTrack()-hit2->GetTrack());
247 hitd->fQ = hit->fQ-hit2->fQ;
248 hitd->fSector = hit->fSector-hit2->fSector;
249 }
250
251 if (debug){
252 Float_t dd =
253 TMath::Sqrt(
254 (hit->X()-hit2->X())*(hit->X()-hit2->X())+
255 (hit->Y()-hit2->Y())*(hit->Y()-hit2->Y())+
256 (hit->Z()-hit2->Z())*(hit->Z()-hit2->Z()));
257 printf("C1\t%d\t%d\t%d\t%d\t",
258 hit->fSector,hit2->fSector,hit->GetTrack(),hit2->GetTrack());
259 printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
260 hit->X(),hit2->X(), hit->Y(),hit2->Y(),hit->Z(),hit2->Z());
261 printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
262 dd, param->fR,param->fZ,param->fFi,param->fAn,
263 param->fAd,param->fTheta);
264 printf("%d\t%d\t%d\n",
265 (Int_t)info->fHitDistance, (Int_t)hit->fQ, (Int_t)hit2->fQ);
266 }
267 }
268 }
269 }
270}
271
272
273
274
275//Filter for paw visualisation of results
276//
277//sed filter
278 //sed '/*/d' dd.txt | sed '/C/d'| cat -n > out.txt
279 // sed -n '/C2/p' dd.txt | sed 's/C1//' | cat -n >out2.txt
280 // sed -n '/C3/p' dd.txt | sed 's/C2//' | cat -n >out3.txt
281
282//filter
283//myfilter C1 dd.txt >out1.txt
284//myfilter C2 dd.txt >out2.txt
285//myfilter C3 dd.txt >out3.txt
286// sed -n 1,50000p
287// sed -n 1,50000p dd.txt | myfilter C1 | sed -n 1,50000p >out1.txt
288
289//myfilter C1 dd.txt | sed -n 1,50000p >out1.txt
290//myfilter C2 dd.txt | sed -n 1,50000p >out2.txt
291//myfilter C3 dd.txt | sed -n 1,50000p >out3.txt
292/*
293paw visualisation
294 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
295 Nt/Read 1 out1.txt
296 nt
297 Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r z fi alfa theta dr1 ddr ddz ddrfi dd
298 Nt/Read 2 out2.txt
299
300 Nt/Create 3 'Test of Match2' 2 ! ! event stacksize
301 Nt/Read 3 out3.txt
302*/
303
304void Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
305 Double_t fSumX, Double_t fSumX2, Double_t fSumX3,
306 Double_t fSumX4, Int_t n,
307 Double_t &a, Double_t &b, Double_t &c)
308{
309 //
310 //recalc parameters not fixing origin point
311 Double_t det =
312 n* (fSumX2*fSumX4-fSumX3*fSumX3) -
313 fSumX* (fSumX*fSumX4-fSumX3*fSumX2)+
314 fSumX2* (fSumX*fSumX3-fSumX2*fSumX2);
315
316 if (TMath::Abs(det)>0){
317 a =
318 (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
319 fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
320 fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
321 b=
322 (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
323 fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
324 fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
325 c=
326 (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
327 fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
328 fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
329 cout<<a<<"\t"<<b<<"\t"<<c<<"\n";
330 }
331}
332
333void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
334{
335 Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
336 fSumY = fSumYX = fSumYX2 = fSumX = fSumX2 = fSumX3 = fSumX4 =0;
337 Double_t a2,b2,c2;
338 for (Int_t i=0;i<n; i++){
339 Float_t x = gRandom->Rndm();
340 Float_t y = a+b*x+c*x*x;
341 fSumY+=y;
342 fSumYX+=y*x;
343 fSumYX2+=y*x*x;
344 fSumX += x;
345 fSumX2 += x*x;
346 fSumX3 += x*x*x;
347 fSumX4 += x*x*x*x;
348 }
349 Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);
350}