CMake: Adding missing ALI_DATE flag
[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 /// Macro to compare TClonesArray hits with interpolated hits
17 /// ConvertHits1 read 
18 ///
19 /// \author MI
20
21 #include "alles.h"
22 #include "AliObjectArray.h"
23 #include "AliTPCTrackHits.h"
24 #include "AliArrayBranch.h"
25 #include "TestTPCTrackHits.h"
26  
27 /// \cond CLASSIMP
28 ClassImp(AliTPChitD)
29 /// \endcond
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 {      
36
37   /// convert - not parametrised hits stored in Clones array
38   /// to
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");
63   treeh3->Branch("TPC", &arr,64000,100);
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);
149   treeh3->Branch("TPC", &arrd,64000,100);
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 }