CMake: Changing DA detector name to DA0 to match DAQ settings
[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 /// \class TestTPCTrackHits
17 /// Macro to compare TClonesArray hits with interpolated hits
18 /// ConvertHits1 read
19 ///
20 /// \author MI
21
22 #include "alles.h"
23 #include "AliObjectArray.h"
24 #include "AliTPCTrackHits.h"
25 #include "AliArrayBranch.h"
26 #include "TestTPCTrackHits.h"
27  
28 /// \cond CLASSIMP
29 ClassImp(AliTPChitD)
30 /// \endcond
31  
32 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits,  Bool_t debug, TClonesArray *arrd=0);
33 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
34
35 void ConvertHits(const char * benchmark, Bool_t debug)
36 {      
37
38   /// convert - not parametrised hits stored in Clones array
39   /// to
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");
64   treeh3->Branch("TPC", &arr,64000,100);
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
97 void CompareHits(const char * benchmark, Bool_t debug)
98 {        
99   /// compare persistent hits
100
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
131 void CompareHitsG(const char * benchmark, Bool_t debug)
132 {        
133   /// compare persistent hits
134
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);
152   treeh3->Branch("TPC", &arrd,64000,100);
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    
176 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)  
177 {
178   /// make track wit hits  according
179
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
211 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
212 {     
213   /// if debug option  kTRUE
214   /// compare hits and write result  to the stdoutput
215
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 /*
296 paw 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
307 void   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 {
312   /// recalc parameters not fixing origin point
313
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
336 void 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 }