Default Branch split level set to 99
[u/mrichter/AliRoot.git] / TPC / TestTPCTrackHits.cxx
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
16 /*
17 $Log$
18 Revision 1.3  2000/11/02 10:22:50  kowal2
19 Logs added
20
21 */
22
23 /*
24   Author : MI
25   macro to compare TClonesArray hits with interpolated hits
26   ConvertHits1 read 
27 */
28
29 #include "alles.h"
30 #include "AliObjectArray.h"
31 #include "AliTPCTrackHits.h"
32 #include "AliArrayBranch.h"
33 #include "TestTPCTrackHits.h"
34  
35 ClassImp(AliTPChitD)
36  
37 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits,  Bool_t debug, TClonesArray *arrd=0);
38 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
39
40 void ConvertHits(const char * benchmark, Bool_t debug)
41 {      
42   //
43   //convert - not parametrised hits stored in Clones array
44   //to
45
46
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);
53   //particle  
54   TTree * treeP = (TTree*)f.Get("TreeK0");
55   TBranch *branchP = treeP->GetBranch("Particles"); 
56   TClonesArray *arrp = new TClonesArray("TParticle",100);
57   branchP->SetAddress(&arrp);
58   branchP->GetEvent(0);
59   //
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));
66
67   TFile f3("treehcl.root","recreate");
68   f3.SetCompressionLevel(2);
69   TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
70   treeh3->Branch("TPC", &arr,64000,100);
71
72   gBenchmark->Start(benchmark); 
73   Int_t trsize = (Int_t)treeCl->GetEntries();
74   for (Int_t i=0;i<300;i++){    
75     arr->Clear(); 
76     Int_t size = branch->GetEvent(i);     
77     //if (size>0){
78     if ((size>0)&& arr->GetEntriesFast()>0) {    
79       f.cd();
80       myhits->Clear();
81       myhits =MakeTrack(arr,arrp,myhits);            
82       CompareHits(arr,myhits,debug);
83       f2.cd();
84       treeh2->Fill();
85       f3.cd();
86       treeh3->Fill();      
87       if ((i%10)==0) cout<<i<<"\n"; 
88     }    
89   }
90   f3.cd();
91   treeh3->Write();
92   f2.cd();
93   treeh2->Write();
94   f2.Close();
95   f.Close();
96   delete myhits;
97   gBenchmark->Stop(benchmark);
98   gBenchmark->Show(benchmark); 
99 }  
100
101
102
103 void CompareHits(const char * benchmark, Bool_t debug)
104 {        
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);
111
112     
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);
118    
119   cout<<"Lets go!\n"; 
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);
125     //if (arr){
126     if ((arr->GetEntriesFast()>0) && size>0) {
127       CompareHits(arr,myhits,debug);
128       if ((i%10)==0) cout<<i<<"\n"; 
129     }    
130   }
131   gBenchmark->Stop(benchmark);
132   gBenchmark->Show(benchmark); 
133 }  
134
135
136 void CompareHitsG(const char * benchmark, Bool_t debug)
137 {        
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);
144
145     
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);
151
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);
157    
158   cout<<"Lets go!\n"; 
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);
164     //if (arr){
165     if ((arr->GetEntriesFast()>0) && size>0) {
166       CompareHits(arr,myhits,debug,arrd);      
167     }  
168     if ((i%10)==0) cout<<i<<"\n";
169     treeh3->Fill();         
170   }
171   treeh3->Write();
172   f3.Close();
173   gBenchmark->Stop(benchmark);
174   gBenchmark->Show(benchmark); 
175 }  
176
177
178  
179    
180 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)  
181 {
182   //
183   //make track wit hits  according 
184   AliTPChit * hit;
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);    
192     if (hit){
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);
201       }
202       else {
203         myhits->SetHitPrecision(0.002);
204         myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
205                              hit->Z(), hit->fQ);
206       }
207     }
208   }
209   myhits->FlushHitStack();
210   return myhits;  
211 }
212
213
214
215 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
216 {     
217   //
218   // if debug option  kTRUE
219   // compare hits and write result  to the stdoutput
220   AliTPChit * hit, *hit2;
221   if (arrd) arrd->Clear();
222
223   for (Int_t i=0;i<arr->GetEntriesFast();i++){
224     hit = (AliTPChit*)arr->At(i);  
225     if (hit){
226       if (i==0) myhits->First();
227       else myhits->Next();
228       hit2 = myhits->GetHit();   
229       if (!hit2) {
230         hit2=0;
231         if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
232       }
233
234       if (hit&&hit2){
235         // if (hit){
236         
237         AliTrackHitsParam *param= myhits->GetParam();
238         AliHitInfo * info = myhits->GetHitInfo();
239         //      
240         if (arrd) {
241           TClonesArray &larrd = *arrd;
242           AliTPChitD * h = new(larrd[i]) AliTPChitD;
243           h->SetX(hit->X());
244           h->SetY(hit->Y());
245           h->SetZ(hit->Z());
246           h->SetTrack(hit->GetTrack());
247           h->fQ = hit->fQ;
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;
256         }
257
258         if (debug){
259           Float_t dd = 
260             TMath::Sqrt(
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);
273         }
274       }        
275     }
276   }
277
278
279
280
281
282 //Filter for paw visualisation of results
283 //
284 //sed filter  
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 
288
289 //filter   
290 //myfilter C1 dd.txt >out1.txt
291 //myfilter C2 dd.txt >out2.txt
292 //myfilter C3 dd.txt >out3.txt
293 // sed -n 1,50000p
294 // sed -n 1,50000p dd.txt | myfilter C1  | sed -n 1,50000p >out1.txt
295
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 
299 /*
300 paw visualisation
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
302   Nt/Read 1 out1.txt 
303  nt 
304  Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r  z fi alfa theta dr1 ddr ddz ddrfi dd
305   Nt/Read 2 out2.txt
306  
307  Nt/Create 3 'Test of Match2' 2 ! ! event stacksize 
308   Nt/Read 3 out3.txt     
309 */ 
310
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)
315 {
316   //
317   //recalc parameters not fixing origin point
318   Double_t det = 
319     n* (fSumX2*fSumX4-fSumX3*fSumX3) -
320     fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
321     fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
322     
323   if (TMath::Abs(det)>0){    
324     a = 
325       (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
326        fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
327        fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det; 
328     b=
329       (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
330       fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
331       fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
332     c=
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";
337   }
338 }
339
340 void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
341 {
342   Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
343   fSumY = fSumYX = fSumYX2 = fSumX =  fSumX2 = fSumX3 =  fSumX4 =0;
344   Double_t a2,b2,c2;
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;
348     fSumY+=y;
349     fSumYX+=y*x;
350     fSumYX2+=y*x*x;
351     fSumX += x;
352     fSumX2 += x*x;
353     fSumX3 += x*x*x;
354     fSumX4 += x*x*x*x;
355   }
356   Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);  
357 }