void gammarate()
{  

   gROOT->Reset();
   c1 = new TCanvas("c1","Gamma rate",700,300,500,500);
   c1->SetGridx();
   c1->SetGridy();
   c1->SetLogy(1);
   c1->SetLogx(1);
 
   TFile f("area.root");
   MHMcCollectionArea* mhmc = (MHMcCollectionArea*)f.Get("MHMcCollectionArea");
   TH1D* h1 = (TH1D*)mhmc->GetHist();
   
    
   h1->GetXaxis();
   Float_t xmin = h1->GetXaxis()->GetXmin();
   Float_t xmax = h1->GetXaxis()->GetXmax();
   Int_t nbin = h1->GetNbinsX();

   // gamma spectrum:
   fun1 = new TF1("fun1","10000*2.706e-11*pow(x,-2.47-0.11*log(x))",xmin/1000.,xmax/1000.);

   TH1F *h3f = (TH1F*)h1->Clone("h3f");
   h3f->SetTitle("Integrale per bin");
   h3f->Reset();
   h3f->SetMinimum(1.e-13);
   h3f->SetMaximum(1.e-7);

   TH1F *h4f = (TH1F*)h1->Clone("h4f");

   h4f->SetTitle("Rate");
   h4f->Reset();
   h4f->SetMinimum(1.e-4);
   h4f->SetMaximum(1.e-1);


   Float_t Error = 0.;

   cout << "Bin\tArea\tCenter\tFlux\t\tMolt\tInt\tRate"<<endl; 
   for (Int_t i = 1; i <= nbin; i++)
   {     
     //area->GetXaxis()->FindBin(i)
     Float_t valArea = h1->GetBinContent(i);      
     Float_t bincenter = h1->GetBinCenter(i);
     Float_t valFlusso = fun1->Eval(bincenter/1000., 0., 0.);
     Float_t Integral = fun1->Integral(h1->GetBinLowEdge(i)/1000.,h1->GetBinLowEdge(i+1)/1000.);
     h3f->SetBinContent(i,Integral);
     Float_t Rate = Integral * valArea;
     Float_t ErArea = h1->GetBinError(i);
     Float_t ErRate = ErArea * Integral;
     Error += ErRate*ErRate;
     h4f->SetBinContent(i,Rate);
     h4f->SetBinError(i,ErRate);
     cout << i << "\t" << bincenter << "\t" << Rate << "\t"<< endl;
    }
   //h1->DrawCopy(); 
   //h3f->DrawCopy();
   h4f->SetStats(kFALSE);
   h4f->SetTitle("Gamma rate vs energy");
   h4f->GetYaxis()->SetTitle("Rate");
   h4f->DrawCopy("pe");

   cout << "Rate (Hz):" << h4f->Integral(0,100)<< " +- " <<pow(Error, .5)<<endl;

   h4f->Delete();
   h3f->Delete();
   h1->Delete();
}