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 }