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