.so cleanup: more less-obvious cleanup
[u/mrichter/AliRoot.git] / TPC / AliTPCTestMerge.C
1 #ifndef __CINT__
2 #include "alles.h"
3 #include "AliRunDigitizer.h"
4 #include "AliTPCDigitizer.h"
5 #include "AliH2F.h"
6 #include "TTree.h"
7
8
9 #endif
10
11 // Marian Ivanov
12 // test of the tpc merging using AliRunDigitizer and 
13 // TPC  Hits2Digits, Hits2SDigits and SDigits2Digits macros
14
15 // preparation  
16 // 0. make 2 directorys - ev1 and ev2
17
18 // 1.make hits, digits,sdigits and sdigits-digits in both directories 
19 //  1.a  galice -b -q grun.C and produce hits
20 //  1.b. cp galice.root galice.root.hits
21 //  1.c  run AliTPCHits2Digits.C
22 //  1.d  cp galice.root galice.root.digits
23 //  1.e  copy back cp galice.root.hits galice.root
24 //  1.f  run AliTPCSDigits2Digits.C
25 //  1.g  cp galice.root  galice.root.sdigits
26 //  1.h  run AliTPCSDigits2Digit.C 
27 //  1.i  cp galice.root galice.root.dig2
28
29 // 2. cp ev1/galice.root/galice.root.sdigit galice.root
30 // 3. load this macro and run testmerge()
31
32 // 4. run test function bellow to compare merged digits with original one
33 // 5. to be  noticed output ftom testx() function - should be bigger than
34 //    noise because in proces of digitisation we use different seed
35 //    of random numbers for diffusion and gas gain
36 //    -anyway in place where we have the signal should occur the signal in both casses
37 //    - only the amplitude should be different - by factor of sqrt
38
39       
40 void testmerge()
41 {
42   // merge two example events
43   //
44   //it merge two events -one from current directory -second from directory ev2
45   
46   if(gAlice) delete gAlice;
47   AliRunDigitizer * manager = new AliRunDigitizer(2,1);
48   manager->SetTreeDTPCBaseName("TreeD_75x40_100x60_150x60_");
49   manager->SetInputTreeTPCSBaseName("TreeS_75x40_100x60_150x60_");
50   manager->SetInputStream(0,"galice.root");
51   manager->SetInputStream(1,"ev2/galice.root.sdigits");
52   AliTPCDigitizer *dTPC = new AliTPCDigitizer(manager);
53   manager->SetNrOfEventsToWrite(1); 
54   TStopwatch timer;
55   timer.Start();
56
57   manager->Exec(""); 
58   timer.Stop();
59   timer.Print();
60
61 }
62
63 void drawmerged(Int_t sec, Int_t row, Int_t x1=-1, Int_t x2=-1, Int_t y1=-1, Int_t y2=-1)
64 {
65   //if you think that there is memory leak -
66   //you are tru but othervise graphic doesn't work
67   // sec=0; row =0;
68   TFile * f = new TFile("galice.root");
69   TFile * f1= new TFile("ev1/galice.root.digits");
70   TFile * f2= new TFile("ev2/galice.root.digits");
71   TTree * tree = (TTree*)f->Get("TreeD_75x40_100x60_150x60_0");
72   TTree * tree1 = (TTree*)f1->Get("TreeD_75x40_100x60_150x60_0");
73   TTree * tree2 = (TTree*)f2->Get("TreeD_75x40_100x60_150x60_0");
74   //
75   AliSimDigits *dig=0;
76   AliSimDigits *dig1=0;
77   AliSimDigits *dig2=0;
78   //
79   tree->GetBranch("Segment")->SetAddress(&dig);
80   tree1->GetBranch("Segment")->SetAddress(&dig1);
81   tree2->GetBranch("Segment")->SetAddress(&dig2);
82   AliTPCParam * param =(AliTPCParam*) f->Get("75x40_100x60");
83   if(param){
84     delete param;
85     param=new AliTPCParamSR();
86   }
87   else param=(AliTPCParam*) f->Get("75x40_100x60_150x60");
88   Int_t index = param->GetIndex(sec,row);
89   tree->GetEvent(index);
90   tree1->GetEvent(index);
91   tree2->GetEvent(index);
92
93   TCanvas * c = new TCanvas(" Test merged digits", "test",600,900);
94   c->Divide(1,3);
95   //
96   c->cd(1);
97   AliH2F * his = dig->DrawDigits("cont1",x1,x2,y1,y2);
98   his->SetTitle("MergedDigits");
99   his->SetName("Merged Digits");
100   his->GetXaxis()->SetTitle("time");
101   his->GetYaxis()->SetTitle("pad");
102   his->DrawClone("cont1");
103   //
104   c->cd(2);
105   AliH2F * his1 =dig1->DrawDigits("cont1",x1,x2,y1,y2); 
106   his1->SetTitle("background");
107   his1->SetName("background");
108   his1->GetXaxis()->SetTitle("time");
109   his1->GetYaxis()->SetTitle("pad");
110   his1->DrawClone("cont1"); 
111   //
112   c->cd(3);
113   AliH2F * his2 =dig2->DrawDigits("cont1",x1,x2,y1,y2); 
114   his2->SetTitle("signal");
115   his2->SetName("signal");
116   his2->GetXaxis()->SetTitle("time");
117   his2->GetYaxis()->SetTitle("pad");
118   his2->DrawClone("cont1");  
119 }
120
121
122 void drawd(TFile * f, Int_t amp1, Int_t amp2)
123 {
124   TTree * tree = (TTree*)f->Get("TreeD_75x40_100x60_150x60_0");
125   AliSimDigits *dig=0;
126   tree->GetBranch("Segment")->SetAddress(&dig); 
127   TH1F * his = new TH1F("his","his",amp2-amp1,amp1,amp2);
128   for (Int_t i=0;i<60;i++){  
129     tree->GetEvent(i);   
130     dig->ExpandBuffer();
131     Int_t nrows = dig->GetNRows();
132     Int_t ncols = dig->GetNCols();
133     for (Int_t rows=0;rows<nrows; rows++)
134       for (Int_t col=0;col<ncols; col++){    
135         Int_t d  = dig->GetDigitFast(rows,col);
136         his->Fill(d);
137       }
138   }
139   his->Draw();
140 }
141
142 void test1(){
143   //test of the merged digits
144   //compare merged digits with standard digits
145   TFile f("galice.root");
146   TFile f1("ev1/galice.root.digits");
147   TFile f2("ev2/galice.root.digits");
148   TTree * tree = (TTree*)f.Get("TreeD_75x40_100x60_150x60_0");
149   TTree * tree1 = (TTree*)f1.Get("TreeD_75x40_100x60_150x60_0");
150   TTree * tree2 = (TTree*)f2.Get("TreeD_75x40_100x60_150x60_0");
151   //
152   AliSimDigits *dig=0;
153   AliSimDigits *dig1=0;
154   AliSimDigits *dig2=0;
155   //
156   tree->GetBranch("Segment")->SetAddress(&dig);
157   tree1->GetBranch("Segment")->SetAddress(&dig1);
158   tree2->GetBranch("Segment")->SetAddress(&dig2);
159   //
160   for (Int_t i=0;i<6000;i++){
161     tree->GetEvent(i);
162     tree1->GetEvent(i);
163     tree2->GetEvent(i);
164     dig->ExpandBuffer();
165     dig1->ExpandBuffer();
166     dig2->ExpandBuffer();
167     //
168     Int_t nrows = dig->GetNRows();
169     Int_t ncols = dig->GetNCols();
170     for (Int_t rows=0;rows<nrows; rows++)
171       for (Int_t col=0;col<ncols; col++){
172         Int_t d  = dig->GetDigitFast(rows,col);
173         Int_t d1 = dig1->GetDigitFast(rows,col);
174         Int_t d2 = dig2->GetDigitFast(rows,col);
175         
176         if (abs(d-(d1+d2))>4)
177           printf("%d\t%d\t%d\t%d\t%d\t%d\n",i,rows,col,d,d1,d2);
178       }
179   }
180 }
181
182 void test5(){
183   //
184   //compare merged digits with digits obtained hits2sdig->sdigtodig
185   TFile f("galice.root");
186   TFile f1("ev1/galice.root.dig2");
187   TFile f2("ev2/galice.root.dig2");
188   TTree * tree = (TTree*)f.Get("TreeD_75x40_100x60_150x60_0");
189   TTree * tree1 = (TTree*)f1.Get("TreeD_75x40_100x60_150x60_0");
190   TTree * tree2 = (TTree*)f2.Get("TreeD_75x40_100x60_150x60_0");
191
192   AliSimDigits *dig=0;
193   AliSimDigits *dig1=0;
194   AliSimDigits *dig2=0;
195   
196   tree->GetBranch("Segment")->SetAddress(&dig);
197   tree1->GetBranch("Segment")->SetAddress(&dig1);
198   tree2->GetBranch("Segment")->SetAddress(&dig2);
199
200
201
202
203   for (Int_t i=0;i<6000;i++){
204     tree->GetEvent(i);
205     tree1->GetEvent(i);
206     tree2->GetEvent(i);
207     dig->ExpandBuffer();
208     dig1->ExpandBuffer();
209     dig2->ExpandBuffer();
210     
211     Int_t nrows = dig->GetNRows();
212     Int_t ncols = dig->GetNCols();
213
214     for (Int_t rows=0;rows<nrows; rows++)
215       for (Int_t col=0;col<ncols; col++){
216         Int_t d  = dig->GetDigitFast(rows,col);
217         Int_t d1 = dig1->GetDigitFast(rows,col);
218         Int_t d2 = dig2->GetDigitFast(rows,col);
219         
220         if (abs(d-d1-d2)>4)
221         //if (d2>5)
222           printf("%d\t%d\t%d\t%d\t%d\t%d\n",i,rows,col,d,d1,d2);
223       }
224   }
225 }
226
227 void test3(){
228   //test of the merged digits
229   TFile f("galice.root");
230   TFile f1("ev1/galice.root.sdigits");
231   TFile f2("ev2/galice.root.sdigits");
232   TTree * tree = (TTree*)f.Get("TreeD_75x40_100x60_150x60_0");
233   TTree * tree1 = (TTree*)f1.Get("TreeS_75x40_100x60_150x60_0");
234   TTree * tree2 = (TTree*)f2.Get("TreeS_75x40_100x60_150x60_0");
235   //
236   AliSimDigits *dig=0;
237   AliSimDigits *dig1=0;
238   AliSimDigits *dig2=0;
239   //  
240   tree->GetBranch("Segment")->SetAddress(&dig);
241   tree1->GetBranch("Segment")->SetAddress(&dig1);
242   tree2->GetBranch("Segment")->SetAddress(&dig2);
243   //
244   for (Int_t i=0;i<6000;i++){
245     tree->GetEvent(i);
246     tree1->GetEvent(i);
247     tree2->GetEvent(i);
248     if ( (dig1->GetID()!=i) || (dig2->GetID()!=i) ) {
249       printf("missed segments\n");
250     }
251     //
252     dig->ExpandBuffer();
253     dig1->ExpandBuffer();
254     dig2->ExpandBuffer();
255     //
256     Int_t nrows = dig->GetNRows();
257     Int_t ncols = dig->GetNCols();
258     //
259     for (Int_t rows=0;rows<nrows; rows++)
260       for (Int_t col=0;col<ncols; col++){
261         Int_t d  = dig->GetDigitFast(rows,col);
262         Int_t d1 = dig1->GetDigitFast(rows,col)/16;
263         Int_t d2 = dig2->GetDigitFast(rows,col)/16;     
264         if (abs(d-d1-d2)>4)
265         //if (d2>5)
266           printf("%d\t%d\t%d\t%d\t%d\t%d\n",i,rows,col,d,d1,d2);
267       }
268   }
269 }
270
271
272 void TestSDigitsDig2(){
273   //test of the digits produced by the Hits2Digits 
274   //and Hits2SDigits - SDigits2Digits chain
275   TFile f1("galice.root.digits");
276   TFile f2("galice.root.dig2");
277   //
278   TTree * tree1 = (TTree*)f1.Get("TreeD_75x40_100x60_150x60_0");
279   TTree * tree2 = (TTree*)f2.Get("TreeD_75x40_100x60_150x60_0");
280   //
281   AliSimDigits *dig1=0;
282   AliSimDigits *dig2=0;
283   //  
284   tree1->GetBranch("Segment")->SetAddress(&dig1);
285   tree2->GetBranch("Segment")->SetAddress(&dig2);
286   //
287   for (Int_t i=0;i<6000;i++){
288     //tree->GetEvent(i);
289     tree1->GetEvent(i);
290     tree2->GetEvent(i);
291     //dig->ExpandBuffer();
292     if ( (dig1->GetID()!=i) || (dig2->GetID()!=i) ) {
293       printf("miised semgnets\n");
294     }
295     //        
296     dig1->ExpandBuffer();
297     dig2->ExpandBuffer();
298     dig1->ExpandTrackBuffer();
299     dig2->ExpandTrackBuffer();
300     //
301     Int_t nrows = dig1->GetNRows();
302     Int_t ncols = dig1->GetNCols();
303     //
304     for (Int_t rows=0;rows<nrows; rows++)
305       for (Int_t col=0;col<ncols; col++){
306         Int_t d1 = dig1->GetDigitFast(rows,col);
307         Int_t d2 = dig2->GetDigitFast(rows,col);
308         Int_t t1_1 =dig1->GetTrackIDFast(rows,col,0);
309         Int_t t1_2 =dig1->GetTrackIDFast(rows,col,1);
310         Int_t t1_3 =dig1->GetTrackIDFast(rows,col,2);
311         //
312         Int_t t2_1 =dig2->GetTrackIDFast(rows,col,0);
313         Int_t t2_2 =dig2->GetTrackIDFast(rows,col,1);
314         Int_t t2_3 =dig2->GetTrackIDFast(rows,col,2);
315         //
316         if ( (d2>0) && (d1>0))
317           if ( ( ( d2>2) || (d1>2)) && (t2_1!=t1_1)) {
318             printf("%d\t%d\t%d\t\t%d\t%d\n",i,rows,col,t1_1,t2_1);
319             printf("\t\t\t\t%d\t%d\n",d1,d2);
320           }
321       }
322   }
323 }
324
325 void TestSDigitsDig1(){
326   //test of the digits produced by the Hits2Digits 
327   //and Hits2SDigits - SDigits2Digits chain
328   TFile f1("galice.root.digits");
329   TFile f2("galice.root.dig2");
330   //
331   TTree * tree1 = (TTree*)f1.Get("TreeD_75x40_100x60_150x60_0");
332   TTree * tree2 = (TTree*)f2.Get("TreeD_75x40_100x60_150x60_0");
333   //
334   AliSimDigits *dig1=0;
335   AliSimDigits *dig2=0;
336   //  
337   tree1->GetBranch("Segment")->SetAddress(&dig1);
338   tree2->GetBranch("Segment")->SetAddress(&dig2);
339   //
340   for (Int_t i=0;i<6000;i++){
341     //tree->GetEvent(i);
342     tree1->GetEvent(i);
343     tree2->GetEvent(i);
344     //dig->ExpandBuffer();
345     if ( (dig1->GetID()!=i) || (dig2->GetID()!=i) ) {
346       printf("miised semgnets\n");
347     }
348     //        
349     dig1->ExpandBuffer();
350     dig2->ExpandBuffer();
351     dig1->ExpandTrackBuffer();
352     dig2->ExpandTrackBuffer();
353     //
354     Int_t nrows = dig1->GetNRows();
355     Int_t ncols = dig1->GetNCols();
356     //
357     for (Int_t rows=0;rows<nrows; rows++)
358       for (Int_t col=0;col<ncols; col++){
359         //      Int_t d  = dig->GetDigitFast(rows,col);
360         Int_t d1 = dig1->GetDigitFast(rows,col);
361         Int_t d2 = dig2->GetDigitFast(rows,col);
362  
363         //      
364         if ( ((d2>4) || (d1>4))  && abs(d1-d2)>4)
365           printf("%d\t%d\t%d\t\t%d\t%d\n",i,rows,col,d1,d2);
366       }
367   }
368 }
369
370
371
372 void test4(){
373   //TPC internal test
374   TFile f1("galice.root.sdigits");
375   TFile f2("galice.root.digits");
376   TTree * tree1 = (TTree*)f1.Get("TreeS_75x40_100x60_150x60_0");
377   TTree * tree2 = (TTree*)f2.Get("TreeD_75x40_100x60_150x60_0");
378   //
379   AliSimDigits *dig1=0;
380   AliSimDigits *dig2=0;
381   //
382   tree1->GetBranch("Segment")->SetAddress(&dig1);
383   tree2->GetBranch("Segment")->SetAddress(&dig2);
384   //
385   for (Int_t i=0;i<6000;i++){
386     //tree->GetEvent(i);
387     tree1->GetEvent(i);
388     tree2->GetEvent(i);
389     //dig->ExpandBuffer();
390     if ( (dig1->GetID()!=i) || (dig2->GetID()!=i) ) {
391       printf("miised semgnets\n");
392     }
393     //        
394     dig1->ExpandBuffer();
395     dig2->ExpandBuffer();
396     //
397     Int_t nrows = dig1->GetNRows();
398     Int_t ncols = dig1->GetNCols();
399     //
400     for (Int_t rows=0;rows<nrows; rows++)
401       for (Int_t col=0;col<ncols; col++){
402         Int_t d1 = dig1->GetDigitFast(rows,col)/16.;
403         Int_t d2 = dig2->GetDigitFast(rows,col);                
404         if ((d2>5) &&abs(d1-d2)>2)        printf("%d\t%d\t%d\t\t%d\t%d\n",i,rows,col,d1,d2);
405       }
406   }
407 }
408