// 2D gaussian
Double_t g2(Double_t *x, Double_t *par)
{
   Double_t r1 = Double_t((x[0]-par[1])/par[2]);
   Double_t r2 = Double_t((x[1]-par[3])/par[4]);
   return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2));
}

// First time hnew=true to create, save and handle histos
// Then hnew to read and handle histos
void histo_tmp(Bool_t hnew=true)
{
  char FileName[158];
  char dDir[158];
  char title[158];

  gROOT->LoadMacro("EPTFitDistribution.cc+");

  EPTFitDistribution *fitHisto = new EPTFitDistribution();
  const Int_t npart = 4;   // particles
  const Int_t totDAM = 11; // detectors
  Int_t colors[] = {kRed,kGreen,kBlue,kMagenta};
  TString tcolors[] = {"Red","Green","Blue","Magenta"};
  TString particles[] = {"electron","proton","alpha","gamma"};

  // Histos limits
  Float_t maxS13[] = {30,12,5.5,4,3,3,3,2.5,2.5,2.2,2.2};
  Float_t maxS2[]  = {30,30,6.5,5,4,3,3,2.5,2.5,2.2,2.2};

  // Create histograms
  TH2F **tH2F = new TH2F*[totDAM];
  for (Int_t d=0;d<totDAM;d++) // Detector loop 0 => No Det hit, 1 => First Det ...
  { tH2F[d] = new TH2F[npart]; }

  const Int_t npar = 15;
  TF2 *f2 = new TF2("f2",g2,0,30,0,30, npar);
  Double_t f2params[npar];

  TCanvas *c1 = new TCanvas("c1","Particles distribution",0,0,400,400);

// ========== Plot Parameters ==========
  gPad->SetFillColor(10); // 10 => white
  gPad->SetFillStyle(4000); // 4000 => transparent
  gPad->SetLeftMargin(.15);
  gPad->SetRightMargin(.10);
  gPad->SetBottomMargin(.1);
  gPad->SetTopMargin(0.15);
  gPad->SetFrameBorderSize(0);
  gPad->SetFrameBorderMode(0);
  gPad->SetFrameFillColor(10);
  gPad->SetFrameFillStyle(4000);
  gPad->SetLogx(0);
  gPad->SetLogy(1);


  //=============== if new file create histos and save them
  if (hnew)
  { cout << " generating histos " << endl;
    for (Int_t d=0;d<totDAM;d++) // Detector loop 0 => No Det hit, 1 => First Det ...
    { 
      for (Int_t p=0;p<npart;p++) // particle loop e,p,alpha,gamma
      { tH2F[d][p].SetBins(512, 0, maxS2[d], 512, 0, maxS13[d]);
        sprintf(title,"DAM%02d",d);
        tH2F[d][p].SetTitle(title);

        TRandom *randgen = new TRandom();
        f2params[0] = 1.;
        f2params[1] = maxS2[d]/(4*(p+1))  + randgen->Uniform(maxS2[d]/(8*(p+1)));
        f2params[2] = f2params[1]/4       + randgen->Uniform(f2params[1]/8);
        f2params[3] = maxS13[d]/(4*(p+1)) + randgen->Uniform(maxS13[d]/(8*(p+1)));
        f2params[4] = f2params[3]/4       + randgen->Uniform(f2params[3]/8);
        cout << f2params[0] <<" " << f2params[1] <<" " << f2params[2] <<" " << f2params[3] <<" " << f2params[4] <<" "<<endl;
        f2->SetParameters(f2params);

        tH2F[d][p].FillRandom("f2",1000+TMath::Floor(randgen->Uniform(100)));
      }
    }

    cout << "Saving histos" <<endl;
    SaveHistos(tH2F,totDAM, npart,"epthistos1.root");
  }
  else //============== read histos ================
  { cout << "Reading histos" <<endl;
    ReadHistos(tH2F,totDAM, npart,"epthistos1.root");
  }

  // ==== Find limits ==========

  TAxis *tYAxis;
  TAxis *tXAxis;
  // S2 - S13
  for(Int_t proj=0; proj<2;proj++) // X and Y projection
  { cout << "proj "<<proj<<endl;
    for(Int_t n=1;n<11;n++) // detector loop
    { TH1D *tH1D[3];
      TLine lineb[3];
      TLine linet[3];
      Double_t xmin[3];
      Double_t xmax[3];
      cout << "det "<<n<<endl;
      Int_t first = -1;
      Double_t hmaximorum = 0;

      cout << "================================================="<<endl;
      for (Int_t p=0; p<3; p++)
      { // Projection X or Y with out under and overflow
        // !!!! added different projection names !!!!
        if(proj == 0) 
        { sprintf(title,"%dXp%s",n,particles[p].Data());
          tH1D[p] = tH2F[n][p].ProjectionX(title,1,tH2F[n][p].GetNbinsX()); 
        }
        else
        { sprintf(title,"%dYp%s",n,particles[p].Data());
          tH1D[p] = tH2F[n][p].ProjectionY(title,1,tH2F[n][p].GetNbinsY()); 
        }
        cout <<"title " <<tH1D[p]->GetTitle()<<endl;
        //sprintf(title,"%s %s",tcolors[p].Data(),particles[p].Data());
        cout << title<<endl;

//==================== No problem any more    ===============================
        // Get ROOT calculated max, mean and RMS
        Double_t rhmax  = tH1D[p]->GetMaximum();
        Double_t rmean  = tH1D[p]->GetMean();
        Double_t rrms   = tH1D[p]->GetRMS();
        cout<<"                 MAX      MEAN    RMS"<<endl;
        cout<<"ROOT Calculation "<<rhmax<<" "<<rmean<<" "<<rrms<<endl;

        // Calculate max, mean and RMS
        tYAxis = tH1D[p]->GetYaxis();
        tXAxis = tH1D[p]->GetXaxis();

        Double_t hmax = 0;
        Double_t sum = 0;
        Double_t sum2 = 0;
        Double_t tot = 0;
        for (Int_t binx=1; binx <= tH1D[p]->GetNbinsX(); binx++)
        { Double_t cont = tH1D[p]->GetBinContent(binx);
          Double_t center = tXAxis->GetBinCenter(binx);
          hmax = hmax > cont ? hmax : cont;
          tot += cont;
          sum+= center*cont;
          sum2+= (center*center)*cont;
        }
        hmaximorum = hmaximorum > hmax ? hmaximorum : hmax; 

        Double_t mean = sum/tot;
        Double_t rms  = TMath::Sqrt(TMath::Abs(sum2/tot-(mean*mean)));
       
        cout<<"MY   Calculation "<<hmax<<" "<<mean<<" "<<rms<<endl;

//==================== Here finish the problem ===============================

        xmin[p] = mean-rms*2;
        xmax[p] = mean+rms*2;
        tH1D[p]->SetLineColor(colors[p]);

        tH1D[p]->SetStats(kFALSE);
        if(first<0) { tH1D[p]->Draw("AL"); first=p; }
        else        { tH1D[p]->Draw("samel"); }

      }
      cout << "================================================="<<endl;

      // Draw Limits
      for (Int_t p=0; p<3; p++)
      { // low limit
        lineb[p] = TLine(xmin[p],1,xmin[p],hmaximorum);
        lineb[p].Draw();
        // high limit
        linet[p] = TLine(xmax[p],1,xmax[p],hmaximorum);
        linet[p].Draw();

      }
      tH1D[0]->Draw("sameaxis");
      gPad->Modified();
      gPad->Update();
      c1->Update();

      pause();
    }
  }
}

void SaveHistos(TH2F** histo,Int_t ni, Int_t nj,TString fileName)
{ TFile *tFile = new TFile(fileName,"RECREATE","particle distributions");
  char histoName[20];
  
  for (Int_t i=0; i<ni; i++)
  { for (Int_t j=0; j<nj; j++)
    { sprintf(histoName,"h%02d%02d",i,j);
      histo[i][j].Write(histoName);
    }
  }
  tFile->Close();
}

void ReadHistos(TH2F** histo,Int_t ni, Int_t nj,TString fileName)
{ TFile *tFile = new TFile(fileName);
  char histoName[20];

  for (Int_t i=0; i<ni; i++)
  { for (Int_t j=0; j<nj; j++)
    { sprintf(histoName,"h%02d%02d",i,j);
      histo[i][j] = *(TH2F*) tFile->Get(histoName);
      histo[i][j].SetDirectory(0);
    }
  }
 
  tFile->Close();
}
void pause()
{ cout << "Press ENTER to continue.";
  cin.get();
}