int DEBUG = 1;
double res[2];
void findpeak(TH1D* h){

  TSpectrum *speaks = new TSpectrum(5);
  Int_t nfound = speaks->Search(h, 2, "", 0.10);

  printf("Found %d candidate peaks to fit\n",nfound);
	    
  //printf("Found %d candidate peaks to fit\n",nfound);
  double *xpeaks = speaks->GetPositionX();
  double max = 0.;
  double MaxPeak = 0.;
  int MaxPN;
  for (int pp=0;pp<nfound;pp++) {
    double xp = xpeaks[pp];
    cout<<"Peak "<<pp<<": x = "<<xp<<endl;
    
    Int_t bin = h->GetXaxis()->FindBin(xp);
    double yp = h->GetBinContent(bin);
    if (xp>max){
      MaxPeak = yp;
      MaxPN = pp;
      max = xp;
    }
  }
  
  double maxc = max;
  double bw = h->GetBinWidth(1);
  double loli = maxc - 20.*bw;
  double hili = maxc + 60.*bw;
  if (loli<0){
    loli = 0.;
  }
  
  h->Fit("landau","QR","",loli,hili);
  TF1 *f1 = h->GetFunction("landau");
  double pos = f1->GetParameter(1);
  double sig = f1->GetParameter(2);
  hili = pos + 4.*sig;
  loli = pos - 2.*sig;
  h->Fit("landau","QR","",loli,hili);
  f1 = h->GetFunction("landau");
  pos = f1->GetParameter(1);
  sig = f1->GetParameter(2);
  hili = pos + 4.*sig;
  loli = pos - 2.*sig;
  h->Fit("landau","QR","",loli,hili);
  f1 = h->GetFunction("landau");
  pos = f1->GetParameter(1);      
  sig = f1->GetParameter(2);


  if (DEBUG){
    h->Draw();
    gPad->Update();
    getchar();
  }
  
 
  res[0] = pos;
  res[1] = sig;

  return res;

}


void findmpvs(int RunNum, int det){

  char fnam[128];
  sprintf(fnam,"localdir/run%d/tagger_data_run%d.root",RunNum,RunNum);
  TFile *RF = new TFile(fnam, "READ");

  
  int mic = 0;
  int hod = 0;

  if (det == 1){
    hod = 1;
  } else if (det == 2){
    mic = 1;
  }

  if (mic) {
    
    TH2D *micr = (TH2D*)RF->Get("tagmPeak");
    
    for (int k=1;k<105;k++){
      
      
      TH1D *h = micr->ProjectionX("h",k,k);
      if (h->GetEntries()<100){
	continue;
      }
      cout<<"Counter # "<<k-1<<"    #hits: "<<h->GetEntries()<<endl;
     
      // try to find location of mip peak
      findpeak(h);
      
    }
    
  }

  double DataH[500][2];
  if (hod) {

   TH2D *hodo = (TH2D*)RF->Get("taghPeak");
    
    for (int k=1;k<350;k++){
      DataH[k][0] = 0;
      DataH[k][1] = 0;
      
      TH1D *h = hodo->ProjectionX("h",k,k);
      if (h->GetEntries()<100){
	continue;
      }

      h->GetXaxis()->SetRangeUser(200.,3500.);

      // try to find location of mip peak
      findpeak(h);
      DataH[k][0] = res[0];
      DataH[k][1] = res[1];

    }    
  }

  char outf[128];
  sprintf(outf,"localdir/run%d/hodoscope_peak.dat",RunNum); 
  ofstream OUTF(outf);
  for (int k=1;k<350;k++){
    OUTF<<k<<"   "<<DataH[k][0]<<"   "<<DataH[k][1]<<endl;
  }
  OUTF.close();

}