1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.3 2000/11/02 10:22:50 kowal2
25 macro to compare TClonesArray hits with interpolated hits
30 #include "AliObjectArray.h"
31 #include "AliTPCTrackHits.h"
32 #include "AliArrayBranch.h"
33 #include "TestTPCTrackHits.h"
37 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd=0);
38 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
40 void ConvertHits(const char * benchmark, Bool_t debug)
43 //convert - not parametrised hits stored in Clones array
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);
54 TTree * treeP = (TTree*)f.Get("TreeK0");
55 TBranch *branchP = treeP->GetBranch("Particles");
56 TClonesArray *arrp = new TClonesArray("TParticle",100);
57 branchP->SetAddress(&arrp);
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));
67 TFile f3("treehcl.root","recreate");
68 f3.SetCompressionLevel(2);
69 TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
70 treeh3->Branch("TPC", &arr,64000,100);
72 gBenchmark->Start(benchmark);
73 Int_t trsize = (Int_t)treeCl->GetEntries();
74 for (Int_t i=0;i<300;i++){
76 Int_t size = branch->GetEvent(i);
78 if ((size>0)&& arr->GetEntriesFast()>0) {
81 myhits =MakeTrack(arr,arrp,myhits);
82 CompareHits(arr,myhits,debug);
87 if ((i%10)==0) cout<<i<<"\n";
97 gBenchmark->Stop(benchmark);
98 gBenchmark->Show(benchmark);
103 void CompareHits(const char * benchmark, Bool_t debug)
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);
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);
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);
126 if ((arr->GetEntriesFast()>0) && size>0) {
127 CompareHits(arr,myhits,debug);
128 if ((i%10)==0) cout<<i<<"\n";
131 gBenchmark->Stop(benchmark);
132 gBenchmark->Show(benchmark);
136 void CompareHitsG(const char * benchmark, Bool_t debug)
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);
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);
152 TFile f3("treehdelta.root","recreate");
153 f3.SetCompressionLevel(2);
154 TTree * treeh3 = new TTree("DelataH","DeltaH");
155 TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
156 treeh3->Branch("TPC", &arrd,64000,100);
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);
165 if ((arr->GetEntriesFast()>0) && size>0) {
166 CompareHits(arr,myhits,debug,arrd);
168 if ((i%10)==0) cout<<i<<"\n";
173 gBenchmark->Stop(benchmark);
174 gBenchmark->Show(benchmark);
180 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)
183 //make track wit hits according
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);
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);
203 myhits->SetHitPrecision(0.002);
204 myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
209 myhits->FlushHitStack();
215 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
218 // if debug option kTRUE
219 // compare hits and write result to the stdoutput
220 AliTPChit * hit, *hit2;
221 if (arrd) arrd->Clear();
223 for (Int_t i=0;i<arr->GetEntriesFast();i++){
224 hit = (AliTPChit*)arr->At(i);
226 if (i==0) myhits->First();
228 hit2 = myhits->GetHit();
231 if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
237 AliTrackHitsParam *param= myhits->GetParam();
238 AliHitInfo * info = myhits->GetHitInfo();
241 TClonesArray &larrd = *arrd;
242 AliTPChitD * h = new(larrd[i]) AliTPChitD;
246 h->SetTrack(hit->GetTrack());
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;
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);
282 //Filter for paw visualisation of results
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
290 //myfilter C1 dd.txt >out1.txt
291 //myfilter C2 dd.txt >out2.txt
292 //myfilter C3 dd.txt >out3.txt
294 // sed -n 1,50000p dd.txt | myfilter C1 | sed -n 1,50000p >out1.txt
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
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
304 Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r z fi alfa theta dr1 ddr ddz ddrfi dd
307 Nt/Create 3 'Test of Match2' 2 ! ! event stacksize
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)
317 //recalc parameters not fixing origin point
319 n* (fSumX2*fSumX4-fSumX3*fSumX3) -
320 fSumX* (fSumX*fSumX4-fSumX3*fSumX2)+
321 fSumX2* (fSumX*fSumX3-fSumX2*fSumX2);
323 if (TMath::Abs(det)>0){
325 (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
326 fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
327 fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
329 (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
330 fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
331 fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
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";
340 void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
342 Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
343 fSumY = fSumYX = fSumYX2 = fSumX = fSumX2 = fSumX3 = fSumX4 =0;
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;
356 Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);