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