// do a lognormal fit to an histogram
void lognormal_fit() {

gROOT -> ForceStyle();
gStyle -> SetStatFont(32);
gStyle -> SetStatW(0.20);
gStyle -> SetStatH(0.20);
gStyle -> SetOptStat(1111111);

  // fit to a log normal

  //TFile * f = TFile::Open("mass_T1.root");   // for topology 1
  //TFile * f = TFile::Open("mass_T2.root"); // for topology 2
  TFile * f = TFile::Open("mass_T3.root"); // for topology 3
  //TFile * f = TFile::Open("mass_T4.root"); // for topology 4

  TCanvas * c;
  //f -> GetObject("Topology1 1", c); // for topology 1
  //f -> GetObject("Topology1 2", c); // for topology 2
  f -> GetObject("Topology1 3", c); // for topology 3
  //f -> GetObject("Topology1 4", c); // for topology 4

  delete f; // no longer needed
  c -> Draw(); // must be drawn first

  for (Int_t i = 1; i < 51; i++) { // 50 pads
    cout << "--- Histogram " << i << "---" << endl; 
    c -> cd(i);
    TH1F * h1 = ((TH1F*)(gPad -> FindObject("h1"))); // histogram in "current" pad
    if (h1) {
      h1 -> SetName(h1 -> GetName() + TString("_") + i); // set unique name
      h1 -> SetDirectory(gROOT); // connect to "ROOT Memory"
      //h1 -> Print();

      // fit to a log normal
      TF1 * f1 = new TF1("f1", "[0]*ROOT::Math::lognormal_pdf(x,[1],[2]) ");

      // set initial parameters
      double p[3];
      p[0] = h1->GetEntries() * h1->GetXaxis()->GetBinWidth(1); // area of the histogram
      double area=p[0];

      // find median of histogram
      double prob[] = {0.5};
      double q[1];
      h1 -> GetQuantiles(1, q, prob);
      double median = q[0];

      // find mode of histogram
      double mode = h1 -> GetBinCenter(h1 -> GetMaximumBin());

      std::cout << "histogram mode=" << mode << " median=" << median << " area=" << area << std::endl;

      // m is log(median)
      //p[1] = std::log(median);
      p[1] = log(135.);
      // s2 is  log(median) - log(mode)
      //p[2] = std::sqrt(std::log(median / mode));
      p[2] = 0.15;

      f1 -> SetParameters(p);
      f1 -> SetParName(0, "A");
      f1 -> SetParName(1, "m");
      f1 -> SetParName(2, "s");

      f1->SetParLimits(1, log(110.), log(170.)); // mu
      f1->SetParLimits(2, 0.05, 0.3); // sigma

      //f1->FixParameter(1,log(135));
      //f1->FixParameter(2,0.05);

      //f1 -> Print();
      h1 -> Draw();

      //h1 -> Fit(f1,"V");
      h1->Fit("f1", "QRI", "", 0., 250.);

      double a = f1 -> GetParameter(0);
      double m = f1 -> GetParameter(1);
      double s = f1 -> GetParameter(2);
      cout << "Fitted a=" << a  << " m=" << m <<  " s=" << s << endl;
      cout << "=> mean " << exp(m-s*s/2)  << " median=" << exp(m) <<  " mode=" << exp(m-s*s) << endl;


    } //if h1

  } // i loop

  c -> cd(0);

}//lognormal_fit