#include "TCanvas.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
#include <fstream.h> 
#include "Riostream.h"

Int_t npeaks=9;

void itch_56()
{

  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(1000011);
  gStyle->SetOptFit();

  TFile* foutput = new TFile("/home/bp/rootfiles/whatever.root","RECREATE");

  ofstream cal_file;
  cal_file.open("56eff.dat");

  cal_file <<"Sigma"<<"\t"<<"Energy"<<"\t"<<"Integral"<<endl<<endl;

  Int_t iggy,u,d,j;


  char ziggy[40],filename[20];

 //  cout<<"Now, enter run # you wanna look at, i.e 730: \n";
//   cin>>iggy;



// Chained667_674.root  Chained688_698.root  Chained699_706.root  Chained730_737.root  See below


  iggy=667;

  sprintf(ziggy,"/home/bp/rootfiles/Chained%d_674.root",iggy);

  sprintf(filename,"Run%i_cal.root",iggy);
  TFile* finput = new TFile(ziggy);

  finput->ls();  

//  Clover #  56Co  152Eu  60Co   137Cs

//    0         667   688    699  730
//    1         668   692    700  731
//    2         669   693    701  732
//    3         670   694    702  733
//    4         671   695    703  734
//    5         672   696    704  735
//    6         673   697    705  736
//    7         674   698    706  737

//    Starts at page 27, Logbook 4. Starting the 152Eu now.

  // TObjArray *GRUhi = (TObjArray*)finput->Get("GRUhisto;1");
  // GRUhi->SetName("GRUhi");
  

  TCanvas *can[8];

  TH1I *h;
  TH1I *h2; 
  TH1I *h3; 
  TH1F *myhist[8];
  TH1F *myhist_b[8];
  char hiname[20];
  char canname[20];
  TString tmpstr;
  TString tmpstrfile;
  TString tmpcanstr;
  TString horace;
  TString tmpo;
  //Loop over the clovers 0-7
  for(Int_t i=0;i<8; i++){
      

    // u = (i/4); 
    sprintf(canname,"can_%i",i);
    tmpcanstr = canname;

    //u = (i/4);
    //d = (i%4);

// i = 0 u = 0 and d = 0
// i = 1 u = 0 and d = 1
// i = 2 u = 0 and d = 2
// i = 3 u = 0 and d = 3
// i = 4 u = 1 and d = 0
// i = 5 u = 1 and d = 1
// i = 6 u = 1 and d = 2
// i = 7 u = 1 and d = 3

    can[i] = new TCanvas(tmpcanstr,"fitted histograms",10,10,700,700);
    can[i]->Draw();
    // can[u]->Divide(2,2);
    // can[i]->cd(1);
    
    sprintf(hiname,"hE_C%i",i);

      tmpstr=hiname;
      cout << tmpstr << endl;

      finput->cd();
      h =(TH1I*)gFile->Get(tmpstr);
      h2 = (TH1I*)h->Clone();
      h2->SetName(tmpstr);
      h2->GetXaxis()->SetRangeUser(838,3300);

      h2->ShowBackground();

      horace = "_background";
      tmpo=tmpstr;
      tmpstr+=(horace);

      //    cout<<"Raaaaaaaaaaaaaa i = "<<i<<endl;

      //   TH1F **myhist[i] = (TH1F*)gDirectory->Get(tmpo);


     iff(i==0) TH1F *myhist[0] = (TH1F*)gDirectory->Get(tmpo);
     iff(i==1) TH1F *myhist[1] = (TH1F*)gDirectory->Get(tmpo);
     iff(i==2) TH1F *myhist[2] = (TH1F*)gDirectory->Get(tmpo);
     iff(i==3) TH1F *myhist[3] = (TH1F*)gDirectory->Get(tmpo);
     iff(i==4) TH1F *myhist[4] = (TH1F*)gDirectory->Get(tmpo);
     iff(i==5) TH1F *myhist[5] = (TH1F*)gDirectory->Get(tmpo);
     iff(i==6) TH1F *myhist[6] = (TH1F*)gDirectory->Get(tmpo);
     iff(i==7) TH1F *myhist[7] = (TH1F*)gDirectory->Get(tmpo);

      myhist[i]->Draw();

      //  TH1F *myhist_b[i] = (TH1F*)gDirectory->Get(tmpstr); 

     iff(i==0) TH1F *myhist_b[0] = (TH1F*)gDirectory->Get(tmpstr); 
     iff(i==1) TH1F *myhist_b[1] = (TH1F*)gDirectory->Get(tmpstr);
     iff(i==2) TH1F *myhist_b[2] = (TH1F*)gDirectory->Get(tmpstr);
     iff(i==3) TH1F *myhist_b[3] = (TH1F*)gDirectory->Get(tmpstr);
     iff(i==4) TH1F *myhist_b[4] = (TH1F*)gDirectory->Get(tmpstr);
     iff(i==5) TH1F *myhist_b[5] = (TH1F*)gDirectory->Get(tmpstr);
     iff(i==6) TH1F *myhist_b[6] = (TH1F*)gDirectory->Get(tmpstr);
     iff(i==7) TH1F *myhist_b[7] = (TH1F*)gDirectory->Get(tmpstr);

      myhist_b[i]->Draw();
      myhist[i]->Add(myhist_b[i],-1.);
      myhist[i]->Draw();


//       TH1F *h2_background = (TH1F*)gDirectory->Get("h2_background");

      //  h2->Add(h2_background,-1.);

      //check this /\ (see tired)      

//Use TSpectrum to find the peak candidates
 TSpectrum *s = new TSpectrum(npeaks);
 //vary 0.005 ********
 Int_t nfound = s->Search(h2,1,"new",0.005);
printf("Found %d candidate peaks to fit\n",nfound);

 TF1 *fit[10];
 char funname[20];
 
 cal_file << hiname <<endl; 
 npeaks = 0;
 Float_t *xpeaks = s->GetPositionX();


 for (Int_t p=0;p<nfound;p++) {
   Float_t xp = xpeaks[p];
   Int_t bin = h->GetXaxis()->FindBin(xp);
   Float_t yp = h->GetBinContent(bin);

   npeaks++;
   sprintf(funname,"fit%i",p);
   cout << xp << endl;
   //vary the xp+-*********was 10. Makes no odds.

      h3 = (TH1I*)myhist[i]->Clone();
      h3->SetName(tmpo);
      h3->GetXaxis()->SetRangeUser(xp-11,xp+11);

   fit[p] = new TF1(funname,"gaus",xp-11,xp+11);
   h3->Fit(fit[p],"R+");
   h3->GetXaxis()->SetTitle("chan");
  
   Double_t noddy=(Double_t)(fit[p]->GetParameter(2));

   Double_t bigears=(fit[p]->GetParameter(1));
    
   h3->GetXaxis()->SetRangeUser(bigears-(3*noddy),bigears+(3*noddy));

   //h3->GetXaxis()->SetRangeUser(bigears-(noddy),bigears+(noddy));

  Double_t prob=h3->Integral();

   cout<<"returned integral is "<<prob<<endl;

   cout<<"sigma is "<<noddy<<endl; 
   cout<<"mean is "<<bigears<<endl;
   //  h3->Write();

 cal_file << noddy << "\t" ;
 cal_file << bigears << "\t" ;
 cal_file  << prob << endl;
 }
 
 cal_file << endl;
 foutput->cd();

 h3->Write();
 //h2->Write
    cout<<"nfound = "<<nfound<<endl;
   
    foutput->cd();
    can[i]->Write();
    // tmpstrfile = "calib/new"+tmpstr+".ps";
    //can[i]->SaveAs(tmpstrfile);

    }

  // cout<<"nfound = "<<nfound<<endl;
  cal_file.close();
  
}