void getlandaupeak(TH2* hist, int idsta=1, TString opt=""){
  /*
    ROOT Macro for obtaining Landau peak values from 2d histogram (id vs Qav)
    input: 2d-hist
           option="nstop" -> wait is skipped
    output: landaupeak.txt  -> peak values
  */

  Int_t nbin=hist->GetXaxis()->GetNbins();
  Int_t num=0;
  vector<double> peak;//Landau peak value
  if (nbin>1000){
    cout<<"Too much X-axis bins. Please check the input 2D-histogram."<<endl;
    exit(0);
  }

  for (Int_t i=idsta-1;i<nbin;i++){// begin of loop ------->
    num++;
    Int_t id = i+1;
    TH1 *h=MakeBanY(hist,id);

    Double_t p[3];
    FitHist(h,p);
    peak.push_back(p[1]);

    // Draw
    h->Draw();
    gPad->Update();
    cout<<"ID="<<id<<" peak="<<p[1]<<endl;

    if (opt!="nstop")
      if(Wait()) break;

  }// end of loop <--------------

  // file write
  ofstream ofs1;
  ofs1.open("landaupeak.txt");
  vector<double>::iterator it=peak.begin();
  while(it!=peak.end()){
    ofs1<<*it<<endl;
    it++;
  }
  ofs1.close();

  cout<<num<<" peaks are fitted, END"<<endl;
}
//___________________________________________________________________________
void FitHist(TH1 *h, double p[]){

  int searchedbin[2];
  SearchBins(h,searchedbin);// search local max. and min.

  TF1 *func = new TF1("myfunc","landau");
  double mean = h->GetMean();
  double rms = h->GetRMS();
  func->SetParameters( h->GetMaximum(),
		       h->GetBinCenter(searchedbin[0]),
		       h->GetRMS()      );
  // fitting range
  Double_t xmin,xmax;
  xmin = h->GetBinCenter(searchedbin[1]);
  xmax = h->GetBinCenter(searchedbin[0])+rms;

  h->Fit(func,"Q","",xmin,xmax);
  func->GetParameters(p);
}
//___________________________________________________________________________
void SearchBins(TH1* h, int searchedbin[]){
  /*
    search bin at local maximum (Landau peak) and local minimum 
    from high energy side
  */
  const int range=20;  // for taking average of contents
  int thre = h->GetMaximum()*0.5;// Threshold value 

  Int_t nbin=h->GetXaxis()->GetNbins();
  int isearch=0;
  int sum,sumold=0;
  for(int i=nbin;i>range;i--){
    if (h->GetBinContent(i)<thre) continue;
    sum=h->Integral(i-range,i);
    if (isearch==0 && sum<sumold){// Local maximum
      searchedbin[0]=i;
      isearch=1;
    } else if (isearch==1 && sum>sumold){// Local minimum
      searchedbin[1]=i;
      break;
    }
    sumold=sum;
  }
  if (isearch==0){
    searchedbin[0]=h->GetMaximumBin();
    searchedbin[1]=0;
  }
}
//___________________________________________________________________________
TH1* MakeBanY(TH2* hist, Int_t id){
  TString str=hist->GetName();
  ostringstream oss;
  oss<<id;
  str+="_bny"+oss.str();
  return hist->ProjectionY(str.Data(),id,id);
}
//___________________________________________________________________________
bool Wait(){
  //check and quit if q is typed
  printf("[enter]:continue, [q]:quit >");
  char s[1];
  gets(s);
  if(s[0]=='q') {
    cout<<"Break fit loop"<<endl;
    return true;
  }
  return false;
}
//___________________________________________________________________________