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