//
//   Plot the energy [MIP], nstrip and qtot distributions - Emiliano Mocchiutti
//
//   CaloMIP.c      version 2.12  (2005-11-10)
//
//   The only input needed is the path to the directory created by YODA for the data file you want to analyze.
//
//   Changelog:
//
//   2.11 - 2.12 (2005-11-10): compiled!
//
//   2.10 - 2.11 (2005-08-16): 64 bit arch bugs fixed.
//
//   2.09 - 2.10 (2005-07-25): don't load anymore yodaUtility.c use instead clone of routines
//
//   2.08 - 2.09 (2005-07-12): small changes for compilation.
//
//   2.07 - 2.08 (2005-06-28): changed to be inserted in package CaloQLOOK.
//
//   2.06 - 2.07 (2005-06-07): check if file is older than 050515_007.
//
//   2.05 - 2.06 (2005-03-09): saveas input parameter introduced (required by Maurizio).
//
//   2.04 - 2.05 (2005-03-02): Print out the position of the peak of the fitted function.
//
//   2.03 - 2.04 (2005-02-24): Changed colour background (required by Marco Casolino).
//
//   2.02 - 2.03 (2005-02-24): Changed variable definition from C/C++ style to ROOT style (int->Int_t).
//
//   2.01 - 2.02 (2005-02-07): Added fit of the MIP. Maxevent bug fixed. Changed name of default printed figures.
//
//   2.00 - 2.01 (2005-01-25): Closing the header file at the end it will clear figures on canvas, fixed.
//
//   1.02 - 2.00 (2005-01-20): Now it will access directly LEVEL1 ntuple instead of calibrating everything everytime.
//
//   1.01 - 1.02 (2005-01-17): Cleanup of the code.
//
//   1.00 - 1.01 (2004-12-15): Include also yodaUtility.c .
//
#include <fstream>
#include <sstream>
#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>
#if !defined (__CINT__)				
#include <event/PamelaRun.h>
#include <event/physics/calorimeter/CalorimeterEvent.h>
#include <event/physics/trigger/TriggerEvent.h>
#include <event/CalibCalPedEvent.h>
extern short int CaloLEVEL1(TString, TString, Int_t);
#include <caloclasses.h>
#endif
#include <CaloFunctions.h>


void CaloMIP(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"){
    gROOT->GetListOfCanvases()->Delete("a");
    gDirectory->GetList()->Delete();
    const char* startingdir = gSystem->WorkingDirectory();
#if defined (__CINT__)	
    emicheckLib();
    stringstream llib;
    const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
    llib << pamlib << "/caloclasses_h.so"; 
    gSystem->Load(llib.str().c_str());
#endif
    gStyle->SetOptStat(1111);
    if ( outDir == "" ) outDir = filename;
    //
    // Open LEVEL1 ntuple, if it doesn't exist call CaloLEVEL1.
    //
    TFile       *headerFile = 0;
    TFile       *caloFile = 0;
    headerFile = emigetFile(filename, "Physics", "Header");
    if ( !headerFile ){
    	printf("No header file, exiting...\n");
    	return;
    };
    caloFile = emigetFile(filename, "Physics.Level1", "Calorimeter");
    if ( !caloFile ){
	printf("\n WARNING: No calorimeter LEVEL1 file! \n\n I am going to call CaloLEVEL1, it could take some time to generate level1 data. \n\n");	
	headerFile->Close();
	gSystem->Exec("sleep 3"); 
	short int ERR = CaloLEVEL1(filename,"",0);	
	if ( ERR ) {
	    printf(" I am not able to calibrate data, please call CaloLEVEL1 by hand \n\n Exiting... \n\n");
	    return;
	} else {
	    printf(" OK, back to CaloMIP! \n\n");
	    headerFile = emigetFile(filename, "Physics", "Header");
	    caloFile = emigetFile(filename, "Physics.Level1", "Calorimeter");
	};
    };
    //Takes the tree of the header file
    TTree *otr = (TTree*)headerFile->Get("Pscu");
    otr->AddFriend("CaloLevel1",caloFile);
    //    
    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);
    //
    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());
    };
    //
    //    caloFile->Close();
    //headerFile->Close();
    gSystem->ChangeDirectory(startingdir);
}