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