//
//   Plot the energy [MIP], nstrip and qtot distributions - Emiliano Mocchiutti
//
//   FCaloMIP.cxx    version 1.01  (2006-06-03)
//
//   The only input needed is the path to the directory created by YODA for the data file you want to analyze.
//
//   Changelog:
//
//   1.00 - 1.01 (2006-06-03): Flight version, read unique YODA file.
//
//   0.00 - 1.00 (2006-06-03): Clone of CaloMIP.c v2r12
//
#include <iostream>
#include <sstream>
#include <fstream>
#include <TTree.h>
#include <TClassEdit.h>
#include <TObject.h>
#include <TList.h>
#include <TSystem.h>
#include <TSystemDirectory.h>
#include <TString.h>
#include <TFile.h>
#include <TClass.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TLatex.h>
#include <TPad.h>
#include <TPaveLabel.h>
#include <TChain.h>
#include <TStyle.h>
#include <TF1.h>
//
extern bool existfile(TString);
extern void stringcopy(TString&, const TString&, Int_t, Int_t);
TF1 *langaufit(TH1F *, Double_t *, Double_t *, Double_t *, Double_t *, Double_t *, Double_t *, Double_t *, Int_t *, TString);
Double_t langaupro(Double_t *);
//
#include <caloclassesfun.h>
//
using namespace std;

void FCaloMIP(const TString filename, Int_t view = 0, Int_t plane = 0, Int_t strip = 0,  Int_t fromevent = 0, Int_t toevent = 0, TString outDir = "", TString saveas = "gif"){
    gStyle->SetOptStat(1111);
    const char* startingdir = gSystem->WorkingDirectory();
    if ( !strcmp(outDir.Data(),"") ) outDir = startingdir;
    if ( !existfile(filename) ){
      printf(" %s :no such file or directory \n",filename.Data());
      printf("\n ERROR: No calorimeter LEVEL1 file! \n\n Run FCaloLEVEL1 to generate level1 data. \n\n");	
      return;
    };
    //
    TFile *File = new TFile(filename.Data());
    TTree *otr = (TTree*)File->Get("CaloLevel1");
    if ( !otr ){
      printf(" %s are you sure? I need a LEVEL1 file! \n",filename.Data());
      return;
    };    
    //    
    CalorimeterLevel1 *calo = new CalorimeterLevel1();
    otr->SetBranchAddress("Event", &calo);
    //
    Long64_t nevents    = otr->GetEntries();
    //
    // input limits
    //
    Int_t minview;
    Int_t maxview;
    if ( view > 2 || view < 0 ) {
	printf("You can choose view = 1 (X-views), view = 2 (Y-views) or view = 0 (both)\n");
	return;
    };
    if ( view == 0 ) {
	minview = 0;
	maxview = 2;
    } else {
	minview = view - 1;
	maxview = view;
    };

    if ( plane > 22 || plane < 0 ) {
	printf("You can choose plane between 1 and 22 or plane = 0 (all)\n");
	return;
    };
    Int_t minplane;
    Int_t maxplane;
    if ( plane == 0 ) {
	minplane = 0;
	maxplane = 22;
    } else {
	minplane = plane - 1;
	maxplane = plane;
    };
    Int_t minstrip;
    Int_t maxstrip;
    if ( strip > 96 || strip < 0 ) {
	printf("You can choose strip between 0 and 96 or strip = 0 (all)\n");
	return;
    };
    if ( strip == 0 ) {
	minstrip = 0;
	maxstrip = 96;
    } else {
	minstrip = strip - 1;
	maxstrip = strip;
    };

    if ( fromevent > toevent && toevent ){
	printf("It must be fromevent < toevent \n");
	return;
    };


    if ( fromevent > nevents+1 || fromevent < 0 ) {
	printf("You can choose fromevent between 0 (all) and %i \n",(int)nevents+1);
	return;
    };

    if ( toevent > nevents+1 || toevent < 0 ) {
	printf("You can choose toevent between 0 (all) and %i \n",(int)nevents+1);
	return;
    };
    Int_t minevent;
    Int_t maxevent;
    if ( fromevent == 0 ) {
	minevent = 0;
	maxevent = nevents;
    } else {
	minevent = fromevent - 1;
	if ( toevent > 0 ){
	    maxevent = toevent;
	} else {
	    maxevent = fromevent - 1;
	};
    };

    //
    // Book the histograms:
    //
    //
    TH1F *gmip;
    gmip = new TH1F("the MIP","",55,-0.19,1.9);
    //gmip->SetBit(TH1F::kCanRebin);
    //
    TH1F *diffbase;
    diffbase = new TH1F("diffbase","",55,-0.19,0.19);
    diffbase->SetBit(TH1F::kCanRebin);
    //
    TH1F *Qtot;
    Qtot = new TH1F("qtot","",70,0.,60.);
    Qtot->SetBit(TH1F::kCanRebin);
    //
    TH1F *Nstrip;
    Nstrip = new TH1F("nstrip","",70,0.,60.);
    Nstrip->SetBit(TH1F::kCanRebin);
    //
    // run over all the events
    //
    //
    Int_t tot2 = 0;
    Int_t baseDIFFfull = 0;
    printf("\n Processed events: \n\n");
    //
    for (Int_t i = minevent; i < maxevent; i++){
	if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);

	//
	// read from the header of the event the absolute time at which it was recorded
	//
	otr->GetEntry(i);

	//
	// run over views and planes
	//
	for (Int_t l = minview; l < maxview; l++){
	    for (Int_t m = minplane; m < maxplane; m++){
		Int_t pre = -1; 
		for (Int_t n = minstrip; n < maxstrip; n++){
		    //
		    if ( n%16 == 0 ) pre++;		    
		    if ( n == minstrip && pre == -1 ){
			for (Int_t nm = minstrip; pre == -1 ; nm--){
			    if ( nm%16 == 0 ) pre++;
			};
		    };
		    if ( pre == -1 ) {
			printf("Unable to determine the preamplifier number. Quit.\n");
			return;
		    };
		    //		    
		    if ( calo->nobase ) tot2++;
		    if ( calo->diffbas[l][m][pre] != 0. ) {						
			//			printf("differenza %f %i %i %i %i %i \n",calo->diffbas[l][m][pre],i,l,m,pre,n);
			baseDIFFfull = 1;
			diffbase->Fill(calo->diffbas[l][m][pre]);					    
			//			return;
		    };
		    if ( calo->estrip[l][m][n] != 0. ) gmip->Fill(calo->estrip[l][m][n]);
		};
	    };
	};	    
	if ( calo->nstrip > 0. ) Nstrip->Fill(calo->nstrip);
	if ( calo->qtot > 0. ) Qtot->Fill(calo->qtot);
    };
    if ( tot2 ){
	printf("WARNING: events with no baselines (event lost): %i\n",tot2);
    };
    //
    // figures:
    //
    // THE MIP
    TCanvas *figura = new TCanvas("The_calorimeter_MIP", "The_calorimeter_MIP", 1100, 900);
    figura->SetFillColor(10);
    figura->Range(0,0,100,100);
    figura->cd();
    gPad->SetLogy();
    gmip->SetXTitle("MIP");
    gmip->SetYTitle("Number of events");
    gPad->SetTicks();
    printf("\n Fitting...\n\n");
    // Setting fit range and start values
    Double_t fr[2];
    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
    fr[0] = 0.7;
    fr[1] = 2.0;
    pllo[0]=0.01; pllo[1]=0.5; pllo[2]=1.0; pllo[3]=0.05;
    plhi[0]=5.0; plhi[1]=5.0; plhi[2]=1000000.0; plhi[3]=1.0;
    sv[0]=1.; sv[1]=1.0; sv[2]=500.0; sv[3]=0.2;    
    Double_t chisqr;
    Int_t    ndf;
    TF1 *fitsnr = langaufit(gmip,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf,"QR0");
    //
    Double_t SNRPeak = langaupro(fp);
    printf("\n Conversion factor: %f \n",SNRPeak);
    //    
    //    gStyle->SetOptStat(1111);
    gStyle->SetOptFit(111);
    printf("\n ...done!\n");
    gmip->Draw();
    fitsnr->Draw("lsame");
    // Differences between DSP and calculated baselines in full mode
    TCanvas *base = 0;
    if ( baseDIFFfull ){
	base = new TCanvas("WARNING", "WARNING_full mode,_differences_between_calculated_and_DSP_baselines!!", 1100, 900);
	base->SetFillColor(10);
	base->Range(0,0,100,100);
	base->cd();
	gStyle->SetOptStat(0);
	gPad->SetLogy();
	diffbase->Draw();	
	gStyle->SetOptStat(1111);
    };
    // nstrip and qtot  WARNING: only GOOD strips considered!
    TCanvas *figura2 = new TCanvas("nstrip_and_qtot", "nstrip_and_qtot", 1100, 900);
    figura2->SetFillColor(10);
    figura2->Range(0,0,100,100);
    TPad *pd1;
    TPad *pd2;
    pd1 = new TPad("pd1","This is pad1",0.02,0.02,0.49,0.98,10);
    pd2 = new TPad("pd2","This is pad2",0.51,0.02,0.98,0.98,10);
    figura2->cd();
    gStyle->SetOptStat(1111);
    pd1->Range(0,0,100,100);
    pd2->Range(0,0,100,100);
    pd1->SetGrid();
    pd2->SetGrid();
    pd1->SetTicks();
    pd2->SetTicks();
    pd1->Draw();
    pd2->Draw();
    pd1->cd();
    gPad->SetLogy();
    Nstrip->SetXTitle("Number of hit (E>0.7 MIP)");
    Nstrip->SetYTitle("Number of events");
    Nstrip->Draw();
    pd2->cd();
    gPad->SetLogy();
    Qtot->SetXTitle("MIP");
    Qtot->SetYTitle("Number of events");
    Qtot->Draw();
    printf("\n");
    //    
    const string fil = (const char*)filename;
    Int_t posiz = fil.find("dw_");
    if ( posiz == -1 ) posiz = fil.find("DW_");
    Int_t posiz2 = posiz+13;
    TString file2;
    stringcopy(file2,filename,posiz,posiz2);
    //
    const char *figrec = file2;
    const char *outdir = outDir;
    stringstream figsave;    
    const char *format = saveas;
    figsave.str("");
    figsave << outdir << "/" ;
    figsave << figrec << "_mip1.";
    figsave << format;
    //
    figura->SaveAs(figsave.str().c_str());
    figsave.str("");
    figsave << outdir << "/" ;
    figsave << figrec << "_mip2.";
    figsave << format;
    //
    figura2->SaveAs(figsave.str().c_str());
    if ( baseDIFFfull ){
	figsave.str("");
	figsave << outdir << "/" ;
	figsave << figrec << "_mip3.";
	figsave << format;
	//
	base->SaveAs(figsave.str().c_str());
    };
    //
    gSystem->ChangeDirectory(startingdir);
}