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