]>
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 | ||
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 | 28 | ClassImp(AliTPChitD) |
0fd1fc0b | 29 | /// \endcond |
cf98c13f | 30 | |
31 | void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd=0); | |
32 | AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits); | |
33 | ||
34 | void 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 | ||
96 | void 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 | ||
129 | void 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 | ||
173 | AliTPCTrackHits * 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 | ||
208 | void 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 | /* | |
293 | paw 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 | ||
304 | void 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 | ||
333 | void 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 | } |