/*******************************************
     July 2007    V.L. Morgunov

   Muon calibration using pure Gaussian fit

We should move to the adc channels to get 
MIP pick position, because of no gain, no intercalibration.

Please run it as:  .x gfind_mip.C+ to invoke real compiler
*******************************************/

#include <iostream>
#include <fstream>

#include "TFile.h"
#include "TTree.h"
#include "TBrowser.h"
#include "TRandom.h"
#include "TCanvas.h"
#include <TStyle.h>
#include <TF1.h>
#include <TH1F.h>

using namespace std;

//     Function to fit histograms
TF1 *mg  = new TF1("mg","gaus",-1000.0,5000.0);
//================================================================
Float_t fastfit(TH1F* h1, Float_t &m, Float_t &s){
//================================================================
//      Fit by Gaussian top of MIP distribution 
//    It destroyes of histogram content in memory
// Histogram should consists of track spectrum (ADC-PED) one to one
//     This fit will not work for small gain/light.
// MIP peak shoulde be at least as far as 90 ADC channels from pedestal
//  because of pedestal width is about 100 channels at base
//
// OUTPUT: Chi-sq, Mean, Sigma
//
// Efficiency is better than 99.9 %, see ratios for different voltages
//================================================================
//  Prepare initial histogram
  h1->Rebin(10);
  Float_t entr = h1->Integral();
  h1->Sumw2();
  h1->Scale(1.0/entr);
//------------ Main trik is here  ------------
//  First remove pedestal peak from histogram, but partialy
  for(int k=1;k<=h1->GetNbinsX();k++)
    if(h1->GetBinCenter(k)<50.) // pedestal cut
      h1->SetBinContent(k,0.0);
//  Then find region to fit and remove all other 
  Float_t ymax = 0.3*h1->GetMaximum(); // magic number
  for(int k=1;k<=h1->GetNbinsX();k++)
    if(h1->GetBinContent(k)<ymax)
      h1->SetBinContent(k,0.0);
  Float_t me = h1->GetMean();
  Float_t rm = h1->GetRMS();
  Float_t fit_low  = me-1.5*rm; // magic number
  Float_t fit_high = me+1.0*rm; // magic number
// and Fit the rest of histogram by gaussian function
  mg->SetParameter(0,0.01);
  mg->SetParameter(1,me);
  mg->SetParameter(2,rm);
  h1->Fit(mg,"QR","",fit_low,fit_high);
  m = mg->GetParameter(1); 
  s = mg->GetParameter(2);
  Float_t gchisq = mg->GetChisquare();
  Float_t gndf   = mg->GetNDF();
  if(gndf<1) gndf=1;
  return gchisq/gndf;
}

//================================================================
void read_file(const char *fname, Float_t *result){
//================================================================
  char histo_name[30];
  Float_t me,si,chi;
  TH1F *his = 0;
  TFile *f = TFile::Open(fname);
  if (!f || f->IsZombie()) {
    printf("Cannot open file: %s\n",fname);
    return;
  }
  cout<<" Open  File:  "<<fname<<endl;  
  for(int j=0;j<38;j++){
    cout<<" Fitting of Layer: "<<j<<endl;
    for(int i=0;i<216;i++){
      me  = 0.0;
      si  = 0.0;
      chi = 0.0;
      sprintf(histo_name,"Layer[%02d_%03d]",j,i);
      his = (TH1F*) f->Get(histo_name);
      if (his!=0){
       	chi = fastfit(his, me, si); // Make fit
	if(me<90.0 || me>3000.0){   // strict condition
	  me  = 0.0;
	  si  = 0.0;
	  chi = 0.0;
	}
	result[j*216+i] = me;
      }
    }
  }
  f->Close("R");
  delete f;
  cout<<"           Close File:   "<<fname<<endl;  
}
//================================================================
void gfind_mip() { //  Main program
//================================================================
// Reads files with track histograms in ADC channels and find MIP MPV
//================================================================
  Float_t cf1[38*216];
  Float_t cf2[38*216];
  Float_t cf3[38*216];
  Float_t cf4[38*216];
  Float_t cf5[38*216]; 
  for(int i=0;i<38*216;i++){
    cf1[i] = 0.0;
    cf2[i] = 0.0;
    cf3[i] = 0.0;
    cf4[i] = 0.0;
    cf5[i] = 0.0;
  }
//  At the moment we remember only mean values (easy to change)
  read_file("./minus_03_trk_spect_hist_in_adc.root",cf1);
  read_file("./nominal_trk_spect_hist_in_adc.root", cf2);
  read_file("./plus_03_trk_spect_hist_in_adc.root", cf3);
  read_file("./plus_05_trk_spect_hist_in_adc.root", cf4);
  read_file("./modified_trk_spect_hist_in_adc.root",cf5);

  ofstream fd; // File to write results
  fd.open("all_mu_coeffs.vec");
  cout<<"  Open output File:   all_mu_coeffs.vec"<<endl;
  char s[256];
  for(int j=0;j<38;j++){
    for(int i=0;i<216;i++){
      if(j>29 && (i>71 && i<147))
	continue;
      int ad = j*216+i;
      sprintf(s," %02d  %03d  %8.1f  %8.1f  %8.1f  %8.1f %8.1f\n",
	      j,i,cf1[ad],cf2[ad],cf3[ad],cf4[ad],cf5[ad]);
      fd<<s;
      //      cout<<s;
    }
  }
 fd.close(); 
}