PWGPP-71 - Adding ncluster statistic to optional streamers. Print additional inform...
[u/mrichter/AliRoot.git] / TPC / diganal.C
1 TH1F * h2max = 0;
2 void tpcstart(const  char * param= "Param2",const char * out="out.root",
3               Int_t nevent=0, Float_t th=3)
4 {
5
6
7   gtpc.SetIO("galice.root", "out.root");
8
9   AliTPCD * dig=new AliTPCD;
10   dig.SetName(param);
11   gtpc.GetIn()->cd();
12   dig.Read(param);
13   gtpc.SetDParam(dig);
14   gtpc.SetTree(nevent);
15
16   gtpc.SetbDelHisto(kTRUE);
17   gtpc.SetThreshold(3);
18 }
19   
20
21
22 void tpccfind(Int_t threshold, Float_t thr = 2,
23               Int_t i1 = 2, Int_t i2 = 2, 
24               Int_t tharea =20, Int_t thmax=100)
25 {  
26   ///////////////////////GRAPHICS//DECLARATION//////////////////////
27   TCanvas  * c1 = new TCanvas("padcluster","Cluster finder 1",700,900);
28   c1->cd();
29   TCanvas  * c2 = new TCanvas("padcluster2","Cluster finder 2",700,900);
30   c2->cd();
31   c1->cd();
32   TPad * pad11 = new TPad("pad11","",0.01,0.76,0.48,0.95,21);
33   pad11->Draw();
34   TPad * pad12 = new TPad("pad12","",0.51,0.76,0.95,0.95,21);
35   pad12->Draw();
36   TPad * pad21 = new TPad("pad21","",0.01,0.56,0.49,0.74,21);
37   pad21->Draw();
38   TPad * pad22 = new TPad("pad22","",0.51,0.56,0.95,0.74,21);
39   pad22->Draw();
40   TPad * pad31 = new TPad("pad31","",0.01,0.36,0.49,0.54,21);
41   pad31->Draw();
42   TPad * pad32 = new TPad("pad32","",0.51,0.36,0.95,0.54,21);
43   pad32->Draw();
44   TPad * pad41 = new TPad("pad41","",0.01,0.16,0.49,0.34,21);
45   pad41->Draw();
46   TPad * pad42 = new TPad("pad42","",0.51,0.16,0.95,0.34,21);
47   pad42->Draw();
48
49   c2->cd();
50   TPad * pad11_2 = new TPad("pad11_2","",0.01,0.76,0.48,0.95,21);
51   pad11_2->Draw();
52   TPad * pad12_2 = new TPad("pad12_2","",0.51,0.76,0.95,0.95,21);
53   pad12_2->Draw();
54   /////////////////////HISTOGRAMS///DECLARATION///////////////////////
55   pad11->cd();
56   TH1F * hsx = new TH1F("hsx","Sigma of distribution in time",40,0,2);
57   pad12->cd();
58   TH1F * hsy = new TH1F("hsy","Sigma of distribution in pads",40,0,2);
59
60   pad21->cd();
61   TProfile * hsx2 = new TProfile("hsx2","Sigma of distribution in time",
62                          20,100,500);
63   pad22->cd();
64   TProfile * hsy2 = new TProfile("hsy2","Sigma of distribution in pads",
65                          20,100,500);
66
67   pad31->cd();
68   TH1F * harea = new TH1F("harea","Area of the peak",26,0,52);
69   pad32->cd();
70   TH1F * hmax = new TH1F("hmax","Maximal amplitude in peak",30,0,150);
71  
72
73   pad41->cd();  
74   TProfile * harea2= new TProfile("harea2","Area dependence z coordinata",
75                          20,100,500);
76   pad42->cd();
77   TProfile * hmax2 = new TProfile("hmax2","Maximal amplitude dependence",
78                          20,100,500);
79   pad41->cd();  
80
81   pad11_2->cd();
82   TProfile * harea2p= new TProfile("harea2p","Area dependence on pad coordinata",
83                          20,0,100);
84   pad12_2->cd();
85   TProfile * hmax2p = new TProfile("hmax2p","Maximal amplitude dependence on pad",
86                          20,0,50);
87   //////////////////CALCULATION//////////////////////////////////////////
88
89   for (Int_t k = i1;k <=i2; k++)
90     {
91       tpcanal(1,k,10,0,kFALSE);
92       TClusterFinder * cf=new TClusterFinder(0,0,threshold,1);
93       cf->GetHisto(&gtpc.GetHis1());
94       TClonesArray * arr = cf->FindClusters();
95       cf->Delete();
96       Int_t size = arr->GetEntries();
97       
98       if ( size>0 )   
99         for (Int_t i=0 ; i<size;i++)
100           {
101             Int_t index;
102             TCluster *c=(TCluster*)arr->UncheckedAt(i);
103             hsx->Fill(TMath::Sqrt(c.fSigmaX2));
104             hsy->Fill(TMath::Sqrt(c.fSigmaY2));    
105             if  (TMath::Sqrt(c.fSigmaX2)<thr)
106               hsx2->Fill(c.fX,TMath::Sqrt(c.fSigmaX2),1);
107             if  (TMath::Sqrt(c.fSigmaY2)<thr)
108               hsy2->Fill(c.fX,TMath::Sqrt(c.fSigmaY2),1);       
109             hmax->Fill(c.fMax);
110             harea->Fill(c.fArea);
111             if (c.fArea<tharea)  harea2->Fill(c.fX,c.fArea,1);
112             if (c.fMax<thmax) hmax2->Fill(c.fX,c.fMax,1);
113             if (c.fArea<tharea)  harea2p->Fill(c.fY,c.fArea,1);
114             if (c.fMax<thmax) hmax2p->Fill(c.fY,c.fMax,1);
115             
116           }
117     } 
118   gStyle->SetOptStat(1);
119   pad11->cd();
120   hsx->Draw();
121   pad12->cd();
122   hsy->Draw();
123
124   pad21->cd();
125   hsx2->Draw();
126   pad22->cd();
127   hsy2->Draw();
128
129   pad31->cd();
130   harea->Draw();
131   pad32->cd();
132   hmax->Draw();
133
134   pad41->cd();
135   harea2->Draw();
136   pad42->cd();
137   hmax2->Draw();
138
139
140   c1->cd();
141   TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.14,"NDC");
142   comment->SetTextAlign(12);
143   comment->SetFillColor(42);   
144   comment->ReadFile("clusters.txt");
145   comment->Draw();
146   c2->cd();
147   pad11_2->cd();
148   harea2p->Draw();
149   pad12_2->cd();
150   hmax2p->Draw();
151
152
153 }
154
155
156 void tpcanal(Int_t sec, Int_t row, Int_t pad, Float_t * res=0, 
157              Bool_t bdraw = kTRUE,Float_t xmin=400,Float_t xmax=500)
158 {
159   //calculate occupancy
160   Double_t par[3];
161   gtpc.SetSecRowTime(sec,row);
162   gtpc.SetHisto(pad);
163   //fit occupancy dependence
164
165   gStyle->SetOptStat(0); 
166   g1 = new TF1("pol0_r","pol0",xmin,xmax);  
167   gtpc.GetHis3()->Fit("pol0_r","R0Q");  
168   char s[100];
169   
170   if (bdraw == kTRUE) 
171     {       
172       gtpc.Draw("box");
173       gtpc.GetPad3().cd();
174       gtpc.GetPad3().SetGridx();
175       gtpc.GetPad3().SetGridy();       
176       gtpc.GetPad3().SetLogy();          
177       gtpc.GetPad3().Draw();
178       gtpc.GetPad2().cd();
179       gtpc.GetPad2().SetGridx();
180       gtpc.GetPad2().SetGridy();   
181       fitText = new TPaveText(0.1,0.7,0.4,0.9,"NDC");
182       gtpc.GetHis3()->Draw();
183       sprintf(s,"p0 fit on interval %3.0f+- %3.0f",xmin,xmax);
184       fitText->AddText(s);
185     }
186   g1->GetParameters(&par[0]);
187   Float_t error = g1->GetParError(0);
188
189   sprintf(s,"%0.3f+- %0.3f",par[0],error);
190   if (bdraw == kTRUE)
191     {
192       gtpc.GetHis3()->Fit("pol0_r","R0Q");  
193       fitText->AddText(s);
194       fitText->Draw();
195       gtpc.GetPad2().Update();  
196
197       //plot histograms with specified options
198       //move pads to another position be possible add text
199       gtpc.GetPad1().SetPad(0.05,0.72,0.95,0.95);
200       gtpc.GetPad2().SetPad(0.05,0.47,0.95,0.70); 
201       gtpc.GetPad3().SetPad(0.05,0.22,0.95,0.45);  
202       //add comments to the histograms 
203       gtpc.GetCanvas().cd();
204       TPaveText * comment = new TPaveText(0.05,0.03,0.95,0.2,"NDC");
205       comment->SetTextAlign(12);
206       comment->SetFillColor(42);   
207       comment->ReadFile("comment.txt");
208       comment->Draw();
209     }  
210   if (res != 0)
211     {
212       res[0] = par[0];
213       res[1] = error;
214       cout<<s<<"   "<<res[0]<<"   "<<res[1]<<"\n";
215     }
216 }
217
218
219 void tpcanalall(Int_t isec =1, Int_t lstep = 1,Float_t tmin=400)
220 {
221    AliTPCParam * tpcparam = gtpc->GetParam();
222  //make window for displaying results
223   TCanvas  * c_occu = new TCanvas("coccu","Occupancy dependence",700,900);
224   c_occu->Update();
225   TPad * pad1 = new TPad("occupancy","occupancy",0.05,0.25,0.95,0.95,21);
226   pad1->Draw();
227   //add comments to the histograms 
228   TPaveText * comment = new TPaveText(0.05,0.03,0.95,0.2,"NDC");
229   comment->SetTextAlign(12);
230   comment->SetFillColor(42);
231   comment->ReadFile("comment.txt");
232   comment->Draw();
233   //prepare histogram 
234   Int_t irow = tpcparam->GetNRow(isec);
235   Float_t xmin = tpcparam->GetPadRowRadii(isec,1);
236   Float_t xmax = tpcparam->GetPadRowRadii(isec,irow);  
237   pad1->cd();
238   char s[220];
239   char sh[220];
240   sprintf(s,"occu_sector%d",isec);
241   sprintf(sh,"Occupancy in sector %d as function of pad raw",isec);  
242   TH1F * occudep = new TH1F(s,sh,300,xmin,xmax);
243   Float_t   res[20];
244   Float_t   x;
245   for (Int_t i=2;i<irow;i+=lstep)
246     { 
247       tpcanal(isec,i,10,&res[0],kFALSE,tmin);
248       x = tpcparam->GetPadRowRadii(isec,i) ;
249       Int_t index = (300*(x-xmin))/(xmax-xmin);
250       cout<<i<<"  "<<index<<"   "<<x<<"   "<<res[0]<<"   "<<res[1]<<"\n";  
251       occudep->SetBinContent(index,res[0]);
252       occudep->SetBinError(index,res[1]);      
253     }
254   //plot occupancy histogram
255   
256   pad1->SetGridx();
257   pad1->SetGridy();
258   gStyle->SetOptFit(0);       
259   occudep->Draw("error");
260   occudep->SetXTitle("pad row center position [cm]");
261   occudep->SetYTitle("occupancy"); 
262  
263   //fit occupancy dependence
264   //linear fit
265   TF1 * g1 = new TF1("pol1_r","pol1");  
266   occudep->Fit("pol1_r","+Q"); 
267   Double_t par[3];
268   Float_t error[3]; 
269   Float_t chi;
270   g1->GetParameters(&par[0]);
271   error[0]=g1->GetParError(0);
272   error[1]=g1->GetParError(1);   
273   Float_t  chi = g1->GetChisquare();
274   sprintf(s,"Linear fit     ocupancy = (%2.3f - %2.3f) +(%2.1f+- %2.1f).r   chi2 = %2.2f",
275           par[0],error[0],1000*par[1],1000*error[1],chi);
276   comment->AddText(s);
277  //(1-exp([0]1/(r*2+[1]**2)  fit
278    TF1 * g1 = new TF1("polm1",occur,1,00,1);  
279     occudep->Fit("polm1","+Q"); 
280     Double_t par[3];
281     Float_t error[3]; 
282     g1->GetParameters(&par[0]);
283     error[0]=g1->GetParError(0);
284     //    error[1]=g1->GetParError(1);
285     chi = g1->GetChisquare();
286     sprintf(s,"(1-exp(P1/(x^2) fit   P1=(%2.3f+- %2.3f)    chi2=%2.2f ",
287           par[0],error[0],chi);
288     comment->AddText(s);
289   c_occu->Update();
290     
291 }
292
293 void tpcanalspectra(Int_t isec =1, Int_t r1= 2)
294 {
295    AliTPCParam * tpcparam = gtpc->GetParam();
296  //make window for displaying results
297   TCanvas  * c_occu = new TCanvas("occuhis","Occupancy dependence",700,900);
298   c_occu->Update();
299   TPad * pad1 = new TPad("ocpad1","occupancy1",0.05,0.61,0.95,0.95,21);
300   pad1->Draw();
301   TPad * pad2 = new TPad("ocpad2","occupancy",0.05,0.61,0.95,0.95,21);
302   pad2->Draw();
303
304   TPad * pad3 = new TPad("ocpad3","occupancy1",0.05,0.25,0.95,0.60,21);
305   pad3->Draw();
306   TPad * pad4 = new TPad("ocpad4","occupancy",0.05,0.25,0.95,0.60,21);
307   pad4->Draw();
308
309   //add comments to the histograms 
310   TPaveText * comment = new TPaveText(0.05,0.03,0.95,0.2,"NDC");
311   comment->SetTextAlign(12);
312   comment->SetFillColor(42);
313   comment->ReadFile("comment.txt");
314   comment->Draw();
315   TH2F  his
316   //prepare histogram 
317   for (Int_t i=1;i<irow;i+=lstep)
318     { 
319       tpcanal(isec,i,10,&res[0],kFALSE,tmin);
320     }
321   //plot occupancy histogram
322   
323   pad1->SetGridx();
324   pad1->SetGridy();
325 }
326
327
328
329
330 void tpcdraw(Int_t sec, Int_t row, Int_t pad)
331 {
332    gStyle->SetOptStat(0); 
333   //calculate occupancy for selected sector and pad row 
334   //for selected pad is obtained signal shape 
335   Double_t par[3];
336   gtpc.SetSecRowTime(sec,row);
337   gtpc.SetHisto(pad);
338   gtpc.Draw("box");  
339   //plot histograms with specified options
340   //move pads to another position be possible add text
341   gtpc.GetPad1().SetPad(0.05,0.72,0.95,0.95);
342   gtpc.GetPad2().SetPad(0.05,0.47,0.95,0.70); 
343   gtpc.GetPad3().SetPad(0.05,0.22,0.95,0.45);  
344   //fit histogram of occupancy on specified range <150,500> 
345   gtpc.GetPad2().cd();
346   g1 = new TF1("pol0_r","pol0",150,500); 
347   gtpc.GetHis3()->Fit("pol0_r","R0Q");    
348   g1->GetParameters(&par[0]);
349   Float_t error = g1->GetParError(0);
350   fitText = new TPaveText(0.15,0.7,0.3,0.9,"NDC");
351   fitText->AddText("p0 fit on interval <150-500>");
352   char s[100];
353   sprintf(s,"%0.3f+- %0.3f",par[0],error);
354   fitText->AddText(s);
355   fitText->Draw();
356   gtpc.GetPad2().Update();     
357   //set logarithmic 
358   gtpc.GetPad3().cd();
359   gtpc.GetPad3().SetLogy();  
360    gtpc.GetPad3().Draw();    
361   //add comments to the histograms 
362   gtpc.GetCanvas().cd();
363   TPaveText * comment = new TPaveText(0.05,0.03,0.95,0.2,"NDC");
364   comment->SetTextAlign(12);
365   comment->SetFillColor(42);
366   comment->ReadFile("comment.txt");
367   comment->Draw();
368   gtpc.GetCanvas().Update();
369   
370
371 }
372
373
374 void oDependence()
375 {
376   //set plot options
377   gStyle->SetOptFit(1); 
378   gStyle->SetOptStat(1);  
379   TCanvas  * c1 = new TCanvas("canPRF","Pad response function",700,900);
380   TPad * pad1 = new TPad("pad1THR","",0.05,0.55,0.45,0.95,21);
381   pad1->Draw();
382   TPad * pad2 = new TPad("pad2PRF","",0.55,0.55,0.95,0.95,21);
383   pad2->Draw(); 
384   TPad * pad3 = new TPad("pad3PRF","",0.55,0.05,0.95,0.45,21);
385   pad3->Draw(); 
386  
387   pad1->cd();
388   pad1->SetGridx();
389   pad1->SetGridy();
390   pad2->SetGridx();
391   pad2->SetGridy();
392   pad3->SetGridx();
393   pad3->SetGridy();
394
395   //make histogram of threshold dependence
396   TH1F * hotd =new TH1F("Occupancy dependence on threshold",
397                         "Ocupancy at first pad row as function of threshold",
398                         25,0.,25.);
399
400   hotd->SetBinContent(5,0.625);
401   hotd->SetBinError(5,0.02);
402   hotd->SetBinContent(10,0.559);
403   hotd->SetBinError(10,0.02); 
404   hotd->SetBinContent(20,0.478);
405   hotd->SetBinError(20,0.02);
406   hotd->SetXTitle("Threshold   [channels]");
407   hotd->SetYTitle("occupancy");
408   hotd->Fit("pol1","+");  
409   hotd->Draw("error");
410   //make histogram of PRF  dependence
411   TH1F * hoprfd =new TH1F("Occupancy dependence on PRF width",
412                         "Occupancy at first pad row as function of generic PRF sigma for  2.05x0.35 cm pad size ",
413                         65, 0.,6.5);
414   hoprfd->SetBinContent(10,0.492);
415   hoprfd->SetBinError(10,0.02);
416
417   hoprfd->SetBinContent(20,0.524);
418   hoprfd->SetBinError(20,0.02); 
419
420   hoprfd->SetBinContent(30,0.559);
421   hoprfd->SetBinError(30,0.02);
422   hoprfd->SetXTitle("Sigma of PRF   [mm]");
423   hoprfd->SetYTitle("occupancy");
424   pad2->cd();
425   hoprfd->Fit("pol1","+");  
426   hoprfd->Draw("error");
427   pad2->Draw();
428   //pad 3 histogram  
429   pad3->cd();
430    TH1F * hoprfd88 =new TH1F("Occupancy dependence on PRF width 08x08",
431                         "Occupancy at first pad row as function of generic PRF sigma for  0.8x0.8 cm pad size ",
432                         65, 0.,6.5);
433
434   hoprfd88->SetBinContent(20,0.322);
435   hoprfd88->SetBinError(20,0.02);
436
437   hoprfd88->SetBinContent(30,0.344);
438   hoprfd88->SetBinError(30,0.02); 
439
440   hoprfd88->SetBinContent(40,0.369);
441   hoprfd88->SetBinError(40,0.02);
442  
443   hoprfd88->SetBinContent(60,0.416);
444   hoprfd88->SetBinError(60,0.02);
445   hoprfd88->SetXTitle("Sigma of PRF   [mm]");
446   hoprfd88->SetYTitle("occupancy");
447   hoprfd88->Fit("pol1","+");  
448   hoprfd88->Draw("error");
449   c1->cd();
450   TPaveText * comment = new TPaveText(0.05,0.15,0.45,0.35,"NDC");
451   comment->SetTextAlign(12);
452   comment->SetFillColor(42);  
453   comment->ReadFile("commentdep.txt");  
454   comment->Draw();
455
456 }
457
458 void produceHisto()
459
460   TH1F * hmostimp =new TH1F("number most important",
461                         "Mean value number of over threshold digits produced by \
462 most important  particle",
463                         20,0.,23.);
464   gStyle->SetOptStat(0); 
465   hmostimp->Fill(1.,33.51);
466   hmostimp->Fill(5.,34.32);
467   hmostimp->Fill(10.,35.78);
468   hmostimp->Fill(15.,36.85);
469   hmostimp->Fill(20.,37.61);
470   hmostimp->Write();
471   hmostimp->SetMinimum(32.);
472   delete hmostimp;
473   
474   TH1F * hall =new TH1F("number all3",
475                         "Mean value Number of over threshold digits produced by \
476 alll three stored particles",
477                         20,0.,23.);
478   gStyle->SetOptStat(0); 
479   hall->Fill(1.,77.05);
480   hall->Fill(5.,76.905);
481   hall->Fill(10.,71.9);
482   hall->Fill(15.,69.7);
483   hall->Fill(20.,65.15);
484   hall->Write();
485   hall->SetMinimum(32.);    
486 }
487
488 void create6()
489 {
490   TCanvas *cnumber =new TCanvas("Number","Number",600,900);    
491   TPad *pad1 =new TPad("pad1","pad1",0.05,0.7,0.45,0.95);  
492   pad1->cd();    
493   pad1->Draw();
494   gtpc.fout->Get("His_1_1").Draw(); 
495   pad1->Update();
496   TPad *pad2 =new TPad("pad2","pad2",0.55,0.7,0.95,0.95); 
497   TPad *pad3 =new TPad("pad3","pad3",0.05,0.35,0.45,0.66);      
498   TPad *pad4 =new TPad("pad4","pad4",0.55,0.35,0.95,0.66);
499   TPad *pad5 =new TPad("pad5","pad5",0.05,0.05,0.45,0.32);      
500   TPad *pad6 =new TPad("pad5","pad6",0.55,0.05,0.95,0.32);   
501 }
502
503 Double_t xtom1(Double_t *x, Double_t *par)
504
505   Double_t xm=x[0]/100.;
506   return (par[0]/(xm*xm));
507 }
508
509 Double_t xtomo(Double_t *x, Double_t *par)
510
511   Double_t xm=x[0]/100.;
512   return (par[0]/(xm**par[1]));
513 }
514
515 Double_t occur(Double_t *x, Double_t *par)
516 {
517   Double_t xm=x[0]/100.;
518   return  (1-exp(-par[0]/(xm**2)));
519
520
521 void probability(Float_t param = 1, Float_t x1 = 0.9, Float_t x2 = 1.3, 
522                  Float_t over = 2,const char * com ="", Int_t N=20)
523 {  
524
525   //create canvas for drawing
526   TCanvas  * c1 = new TCanvas("canprob","Pad response function",700,900);
527   TPad * pad1 = new TPad("Theoretical probability","",0.05,0.22,0.95,0.95,21);
528   pad1->Draw();
529    
530   //create histogram with estimated occupancy 
531   //normalised to obtained occupancy at first pad
532   Float_t y1=0;
533   Float_t y2;
534   char s[120];
535   sprintf(s,"1-exp(-[1]/x**%f)",over);
536   cout<<s<<"\n";
537   TF1 *funr1 = new TF1("funr1",s,x1,x2);
538   funr1->SetParameters(1,param);
539   sprintf(s,"Probability  according 1-exp(-%1.3f/x**%1.1f distribution)",
540           param,over);
541   pad1->cd();
542   TH1F * hOccuT = new TH1F("hOccuT",s,5*N,x1,x2);
543   Float_t x=x1;
544   Float_t y;
545
546   for (Float_t i = 0;i<N+1;i++)
547     {     
548       y = funr1->Eval(x);
549       hOccuT->Fill(x,y);
550       x+=(x2-x1)/Float_t(N);
551     };
552   //fitting calculated dependence with linear fit and with
553   //generic function
554   pad1->cd();
555   sprintf(s,"[1]/(x**%1.1f)",over);
556   TF1 *lin1 = new TF1("lin1","pol1",x1,x2);
557   lin1->SetLineColor(2);
558   lin1->SetLineWidth(5);
559   lin1->SetLineStyle(1);
560   hOccuT->Draw();
561   hOccuT->Fit("lin1","S+");
562
563   sprintf(s,"[1]/(x**%1.1f)",over);
564   TF1 *funorig = new TF1("funorig",s,x1,x2);
565   funorig->SetLineColor(3);
566   funorig->SetLineWidth(5);
567   funorig->SetLineStyle(2);
568   hOccuT->Fit("funorig","S+");
569   
570   //find minimum and maximum and scale histo  
571   if (y1 == 0)  
572     {
573       Float_t ymin,ymax;
574       y1=lin1->Eval(x2);
575       y2=lin1->Eval(x1);
576       ymin= funorig->Eval(x2);
577       ymax= funorig->Eval(x1);
578       if (ymin<y1) y1=ymin;
579       if (ymax>y2) y2=ymax;
580     }
581   gStyle->SetOptFit(0);
582   gStyle->SetOptStat(0); 
583   hOccuT->SetMaximum(y2);
584   hOccuT->SetMinimum(y1); 
585   hOccuT->SetXTitle("r position [m]");
586   hOccuT->SetYTitle("probability");  
587   //add comments to the histograms 
588   c1->cd();
589   TPaveText * comment = new TPaveText(0.05,0.03,0.95,0.20,"NDC");
590   comment->SetTextAlign(12);
591   comment->SetFillColor(42);
592   TText *title = comment->AddText("Estimation of occupancy dependence on R position");
593   title->SetTextSize(0.04);
594   comment->AddText("Observed efect of probability saturation");
595   sprintf(s,"Supposed generic flux dependence : %1.3f/(x**%1.1f)",param,over);
596   comment->AddText(s);
597   comment->AddText("Probility : 1-exp(-flux*mean particle \"pad x time\" area)");
598   comment->AddText("Full line  linear fit ");
599   comment->AddText("Dashed line : fit by generic flux function ");
600   
601   comment->AddText(com);
602   comment->Draw();
603 }
604
605
606
607 void digamp(Float_t t1, Float_t t2, Float_t p1, Float_t p2, Float_t th)
608 {
609   gStyle->SetOptStat(1); 
610   TH1F * h2max = new TH1F("all amplitudes","all amplitude", 20, 1,600);
611   for (Int_t itime = t1;itime<t2;itime++)
612   for (Int_t ipad = p1;ipad<p2;ipad++)  
613     {
614       Int_t index = gtpc.GetHis1().GetBin(itime,ipad);
615       Float_t weight = gtpc.GetHis1().GetBinContent(index);
616       //      cout<<itime<<"\t"<<ipad<<"\t"<<weight<<"\n";
617       if (weight > th) h2max->Fill(weight);
618     };  
619    h2max->Draw();
620 }            
621
622 void digampM(Float_t t1, Float_t t2, Float_t p1, Float_t p2, Float_t th)
623 {
624   gStyle->SetOptStat(1); 
625
626   //create canvas for drawing
627   //  TCanvas  * c1 = new TCanvas("dh","Amplitude",700,900);
628   //TPad * pad1 = new TPad("amplitude","",0.05,0.22,0.95,0.95,21);
629   //pad1->Draw();
630
631   if (h2max == 0) h2max = new TH1F("max amplitudes","max amplitude", 25, 1,250);
632   for (Int_t itime = t1;itime<t2;itime++)
633   for (Int_t ipad = p1;ipad<p2;ipad++)  
634     {
635       Bool_t bmax = kTRUE;
636       Int_t index = gtpc.GetHis1().GetBin(itime,ipad);
637       Float_t weight = gtpc.GetHis1().GetBinContent(index);
638       if (weight>th)
639         {
640           for (Int_t i = -1;i<2;i++) 
641             for (Int_t j = -1;j<2;j++)    
642               if (!((0==i)&&(0==j)))
643                 {
644                   index = gtpc.GetHis1().GetBin(itime+i,ipad+j);
645                   Float_t weightl = gtpc.GetHis1().GetBinContent(index);
646                   if (!(weightl<weight)) bmax = kFALSE;
647                 }     
648           if (kTRUE==bmax) h2max->Fill(weight);
649         }
650     };  
651    h2max->Draw("error");
652
653 }            
654
655
656
657 void digampMALL(Float_t t1, Float_t t2, Float_t p1, Float_t p2, Float_t th )
658 {
659   for (Int_t i=1;i<20;i++)
660     {
661       tpcanal(1,i,10);
662       digampM(t1,t2,p1,p2,th);
663     }
664 }