//
//   Given a calibration and a data file this program create an ntuple with LEVEL1 calorimeter variables - Emiliano Mocchiutti
//
//   FCaloLEVEL1.cxx    version 1.01  (2006-03-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-03-03): read unique YODA file and compile correctly with ROOT classes etc.
//
//   0.00 - 1.00 (2006-03-03): clone of CaloLEVEL1.c v 1.21.
//
#include <fstream>
#include <sstream>
#include <TTree.h>
#include <TBranch.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>
extern const char *pathtocalibration();
#include <PamelaRun.h>
#include <physics/calorimeter/CalorimeterEvent.h>
#include <physics/trigger/TriggerEvent.h>
#include <CalibCalPedEvent.h>
#include <fcalostructs.h>
//
extern char *getLEVname(TString, TString, TString, TString);
extern TString whatnamewith(TString, Int_t);
extern void stringcopy(TString&, const TString&, Int_t, Int_t);
extern int CaloPede(TString, Int_t, Int_t, Calib &);
extern void CaloFindBaseRaw(Calib &, Evento &, Int_t, Int_t, Int_t);
extern void CaloCompressData(Calib &, Evento &, Int_t, Int_t, Int_t);
extern int CaloFindCalibs(TString &, TString &, Int_t &, Calib &);
extern bool existfile(TString);
//
#include <caloclassesfun.h>

short int FCaloLEVEL1(TString filename, TString OutDir="", TString calcalibfile = "", Int_t FORCE = 0){
    const char* startingdir = gSystem->WorkingDirectory();
    if ( !strcmp(OutDir.Data(),"") ) OutDir = startingdir;
    stringstream calfile;
    const char *pam_calib = pathtocalibration();
    calfile.str("");
    calfile << pam_calib << "/FCaloADC2MIP.dat";
    //
    if ( !existfile(filename) ){
      printf(" %s :no such file or directory \n",filename.Data());
      return(1);
    };
    TFile *File = new TFile(filename.Data());
    TTree *otr = (TTree*)File->Get("Physics");
    otr->SetBranchStatus("*",0); //disable all branches
    otr->SetBranchStatus("iev*",1);
    otr->SetBranchStatus("dexy*",1);
    otr->SetBranchStatus("stwerr*",1);
    otr->SetBranchStatus("perror*",1);
    otr->SetBranchStatus("cal*",1);
    otr->SetBranchStatus("base*",1);
    otr->SetBranchStatus("Pscu*",1);
    otr->SetBranchStatus("Calorimeter*",1);
    otr->SetBranchStatus("Header*",1);
    //
    pamela::calorimeter::CalorimeterEvent *de = 0;   
    pamela::PscuHeader *ph = 0;
    pamela::EventHeader *eh = 0;
    //
    otr->SetBranchAddress("Calorimeter", &de);
    otr->SetBranchAddress("Header", &eh);
    //
    // calibration file
    //
    if ( calcalibfile == "" ) calcalibfile = filename;
    //
    struct Evento evento;
    struct Evento evento2;
    struct Calib calib;
    evento.emin = 0.7;
    Int_t tot0 = 0;
    Int_t tot1 = 0;
    Int_t tot2 = 0;
    Int_t baseDIFFfull = 0;
    //
    // Define variables
    // 
    Long64_t nevents    = otr->GetEntries();
    if ( nevents < 1 ) {
	printf("The file is empty!\n");
	return(1);
    };
    Float_t mip[2][22][96];
    Int_t b[4];
    Int_t etime,evno,stwerr[4];
    //    Int_t etime,pl,evno,stwerr[4];
    //    Float_t ener, basel,sbase[2][22][6],sdexy[2][22][96],sdexyc[2][22][96];
    Float_t ener, basel,sdexy[2][22][96],sdexyc[2][22][96];
    Float_t esup[2][22][96], einf[2][22][96], edelta[2][22][96];
    //
    // create the ntuple level1 name
    //
    //    file = "dw_000000_00000.Physics.Level1.Calorimeter.Event.root"; // just to make space...
    TString numb = "1";  // level 
    TString detc = "Calorimeter";  // detector   
    char *file = getLEVname(filename,detc,numb,"root");
    //printf("file> %s\n",file);
    //return;
    //
    //    char *file2; // the full path and name of the level1 ntuple
    stringstream file2;
    const char *file3;
    file3 = OutDir;
    file2.str("");
    file2 << file3 << "/";
    file2 << file;
    //  file2 = Form("%s/Physics/Level1/%s",file3,file);  // add the path where to save
    printf("\n Filename will be: \n %s \n\n",file2.str().c_str());
    //
    // check if Level1 directory exists, if not create it.
    //    
    stringstream savedir;
    savedir.str("");
    savedir << file3 << "/";
	//    char *savedir;
	//savedir = Form("%s/Physics/Level1");
    //    Int_t ERR = gSystem->OpenDirectory(savedir.str().c_str());
    //
    // check if Level1 directory exists, if not create it.
    //    
    Int_t ERR = gSystem->MakeDirectory(savedir.str().c_str());    
    //
    if ( !ERR ) {
	printf(" LEVEL1 directory doesn't exist! Creating LEVEL1 directory \n\n");
	gSystem->MakeDirectory(savedir.str().c_str());
    } else {
	//
	// directory exists, check if file exists (if we are not in the force mode)
	//    
	if ( !FORCE ) {
	    printf(" Not in FORCE mode, check the existance of LEVEL1 data: \n\n");
	    TFile tfile(file2.str().c_str());
	    if ( !tfile.IsZombie() ) {
		printf(" ERROR: file already exists! Use FORCE = 1 to override \n\n");
		return(3);
	    } else {
		printf("\n OK, I will create it!\n");
	    };
	};
    };
    char *type;
    type = "NEW";
    if ( FORCE ) type = "RECREATE";
    TFile *hfile = 0;
    hfile = new TFile(file2.str().c_str(),type,"Calorimeter LEVEL1 data");
    //
    // variables which will fill the ntuple
    //
    Float_t perror[4];
    Float_t nstrip;
    Float_t qtot;
    Float_t calevnum[4];   
    Float_t estrip[2][22][96];
    Float_t estripnull[2][22][96];
    Float_t diffbas[2][22][6];
    // 
    // class CalorimeterLevel1
    //
    CalorimeterLevel1 *calo = new CalorimeterLevel1();
    //
    // Create a ROOT Tree
    TTree *tree = 0;
    //    calo = new CalorimeterLevel1();
    tree = new TTree("CaloLevel1","PAMELA Level1 calorimeter data");
    tree->Branch("Event","CalorimeterLevel1",&calo);
    //TBranch *branch = tree->Branch("Event",&calo,64000,0);
    //    
    // this is the value of the mip for each strip. To be changed when we will have the real values
    //
    for (Int_t s=0; s<4;s++){
	for (Int_t d = 0; d<51; d++){
	    calib.ttime[s][d] = 0 ;
	    calib.time[s][d] = 0 ;
	};
    };
    //
    Int_t okcalo = 0;
    printf(" ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
    FILE *f;
    f = fopen(calfile.str().c_str(),"rb");
    if ( !f ){
	printf(" WARNING: no calorimeter ADC to MIP file! \n Using 26 as conversion factor for all strips. \n");
    } else {
	okcalo = 1;
    };
    //
    if ( okcalo ) {
	for (Int_t m = 0; m < 2 ; m++ ){
	    for (Int_t k = 0; k < 22; k++ ){
		for (Int_t l = 0; l < 96; l++ ){
		  fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);
		  calib.calped[m][k][l] = 0. ;
		  evento.dexy[m][k][l] = 0. ;
		  evento.dexyc[m][k][l] = 0. ;
		  estrip[m][k][l] = 0.;
		  estripnull[m][k][l] = 0.;
		  einf[m][k][l] = 32000.;
		  esup[m][k][l] = 0.;
		  edelta[m][k][l] = 0.;
		};
	    };
	};
    } else {
	for (Int_t m = 0; m < 2 ; m++ ){
	    for (Int_t k = 0; k < 22; k++ ){
		for (Int_t l = 0; l < 96; l++ ){
		  mip[m][k][l] = 26. ;
		  calib.calped[m][k][l] = 0. ;
		  evento.dexy[m][k][l] = 0. ;
		  evento.dexyc[m][k][l] = 0. ;
		  estrip[m][k][l] = 0.;
		  einf[m][k][l] = 32000.;
		  esup[m][k][l] = 0.;
		  edelta[m][k][l] = 0.;
		};
	    };
	};
    };
    if ( okcalo ) fclose(f);
    //    chfile->Close();
    //
    // first of all find the calibrations in the file 
    //
    Int_t wused = 0;
    CaloFindCalibs(filename, calcalibfile, wused, calib);        
    if ( wused == 1 ) calcalibfile = filename;
    //
    // print on the screen the results:   
    //
    const char *ffile;
    Int_t done = 0;
    Int_t rdone = 0;
    Int_t fcheck = 0;
    Int_t fdone = 0;
    //    Int_t nn = 0;
    ffile = filename;		
    printf(" ------ %s ------- \n \n",ffile);	
    Int_t calibex = 0;
    TString pfile;
    for (Int_t s=0; s<4;s++){
	printf(" ** SECTION %i **\n",s);
	for (Int_t d = 0; d<51; d++){
	    if ( calib.ttime[s][d] != 0 ) {
		calibex++;
		if ( calib.fcode[s][d] != 10 ){
		    TString file2f = "";
		    stringcopy(file2f,calcalibfile,0,0);
		    pfile = (TString)whatnamewith(file2f,calib.fcode[s][d]);		    
		} else {
		    pfile = (TString)calcalibfile;
		};
		const char *ffile;
		ffile = pfile;		
		printf(" - from time %i to time %i use calibration at\n time %i, file: %s \n",calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d],ffile);
		if ( !strcmp(ffile,"wrong") ) calibex--;
	    };
	};
	printf("\n");	
    };
    printf(" ----------------------------------------------------------------------- \n \n");
    //
    if ( calibex < 4 ) {
	printf("No full calibration data in this file, sorry!\n\n ");
	//
	// remove the empty file!
	//
	stringstream remfile;
	remfile.str("");
	remfile << "rm -f " << file2.str().c_str();
	//	char *remfile;
	//remfile = Form("rm -f %s",file2.str().c_str());
	gSystem->Exec(remfile.str().c_str());
	return(4);
    };
    //
    // calibrate before starting
    //
    for (Int_t s = 0; s < 4; s++){
	b[s]=0;
	if ( calib.fcode[s][b[s]] != 10 ){
	    TString file2f = "";
	    stringcopy(file2f,calcalibfile,0,0);
	    TString pfile = whatnamewith(file2f,calib.fcode[s][b[s]]);		    
	} else {
	    TString pfile(calcalibfile);
	};
	CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
    };	    
    //
    // run over all the events
    //
    printf("\n Processed events: \n\n");
    //
    Int_t se = 5;
    bool isCOMP = 0;
    bool isFULL = 0;
    bool isRAW = 0;
    Int_t upnn = 0;
    Double_t prima = 0.;
    Double_t prmndp = 0.;
    calo = new CalorimeterLevel1();
    for (Int_t i = 0; i < nevents; i++){    
	//for (Int_t i = 0; i < 1000; i++){    
	//calo = new pamela::calorimeter::CalorimeterLevel1();
	evno = i;
	if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);
	//
	//
	//
	qtot = 0.;
	nstrip = 0.;
	//
	// read from the header of the event the absolute time at which it was recorded
	//
	otr->GetEntry(i);
	ph = eh->GetPscuHeader(); 
	etime = ph->GetOrbitalTime();
	//
	// for each event check that the calibration we are using are still within calibration limits, if not call the next calibration
	//
	if ( !calib.obtjump ) {
	    for (Int_t s = 0; s < 4; s++){
		if ( calib.ttime[s][b[s]] ){
		    while ( etime > calib.time[s][b[s]+1] ){				
			printf(" CALORIMETER: \n" );
			printf(" - Section %i, event at time %i while old calibration time limit at %i. Use new calibration at time %i -\n",s,etime,calib.time[s][b[s]+1],calib.ttime[s][b[s]+1]);
			printf(" END CALORIMETER. \n\n" );
			b[s]++;
			if ( calib.fcode[s][b[s]] != 10 ){
			    TString file2f = "";
			    stringcopy(file2f,calcalibfile,0,0);
			    TString pfile = whatnamewith(file2f,calib.fcode[s][b[s]]);		    
			} else {
			    TString pfile(calcalibfile);
			};
			CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
		    };
		};
	    };
	};
	//
	// do whatever you want with data
	//
	evento.iev = de->iev;
	//
	// run over views and planes
	//
	for (Int_t l = 0; l < 2; l++){
	    for (Int_t m = 0; m < 22; m++){
		//
		// determine the section number
		//
	        se = 5;
		if (l == 0 && m%2 == 0) se = 3;
		if (l == 0 && m%2 != 0) se = 2;
		if (l == 1 && m%2 == 0) se = 1;
		if (l == 1 && m%2 != 0) se = 0;		
		//
		// determine what kind of event we are going to analyze
		//
		isCOMP = 0;
		isFULL = 0;
		isRAW = 0;
		if ( de->stwerr[se] & (1 << 16) ) isCOMP = 1; 
		if ( de->stwerr[se] & (1 << 17) ) isFULL = 1; 
		if ( de->stwerr[se] & (1 << 3) ) isRAW = 1; 
		//
		// save the prevoius energy deposit and calibration in sbase, sdexy, sdexyc
		//
		Int_t pre = -1; 
		if ( isRAW ){
		    for (Int_t nn = 0; nn < 96; nn++){		    
			if ( nn%16 == 0 ) pre++;		    
			evento.base[l][m][pre] = calib.calbase[l][m][pre];
			sdexy[l][m][nn] = evento.dexy[l][m][nn];
			evento.dexy[l][m][nn] = de->dexy[l][m][nn] ; 
			sdexyc[l][m][nn] = evento.dexyc[l][m][nn];
			evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ; 
		    };
		};
		//
		// run over strips
		//
		pre = -1;
		for (Int_t n =0 ; n < 96; n++){		    
		    if ( n%16 == 0 ) {
			pre++;
			done = 0;
			rdone = 0;
			fdone = 0;
		    };	
		    // 
		    // baseline check and calculation
		    //
		    if ( !isRAW ) {
			//
			// if it isn't raw and we haven't checked yet, check that the baseline is fine, if not calculate it with a relaxed algorithm.
			//
			if ( !done ){
			    evento.base[l][m][pre] = de->base[l][m][pre] ;   
			    evento.dexyc[l][m][n] = de->dexyc[l][m][n] ; 
			};
		    } else {
			//
			// if it is a raw event and we haven't checked yet, calculate the baseline. Then check that the baseline is fine,
			// if not calculate it with relaxed algorithm.
			//
			if ( !rdone ){
			    CaloFindBaseRaw(calib,evento,l,m,pre);
			    tot1++;
			    rdone = 1;			
			};
		    };					    
		    //
		    // no suitable new baseline, use old ones and compress data!
		    //
		    if ( !done && (evento.base[l][m][pre] == 31000. || evento.base[l][m][pre] == 0.) ){
			tot0++;
			evento.base[l][m][pre] = calib.sbase[l][m][pre];
			upnn = n+16;
			if ( upnn > 96 ) n = 96;
			for ( Int_t nn = n; nn<upnn; nn++ ){
			    evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
			};
			CaloCompressData(calib,evento,l,m,pre);
			done = 1;			
		    };
		    //
		    // if here we don't have a valid baseline we will skip the event.
		    //
		    if ( evento.base[l][m][pre] == 31000. ) tot2++;
		    //
		    //  check the baseline calculation in FULL mode: baseline recorded from DSP must be identical to the one calculated using RAW data.
		    //
		    if ( isFULL ) {
			if ( !fdone) {
			    fdone = 1;
			    if ( de->base[l][m][pre] > 0. && de->base[l][m][pre] < 32000. ) {
				fcheck = 0;
				for (Int_t nn = pre*16; nn < (pre+1)*16 ; nn++){
				    if ( de->dexy[l][m][nn] > 0. &&  de->dexy[l][m][nn] < 32000. ) fcheck++;   
				};	
				if ( fcheck ) {
				    memcpy(&evento2, &evento, sizeof(evento));
				    prima = evento.base[l][m][pre];
				    CaloFindBaseRaw(calib,evento2,l,m,pre);	
				    diffbas[l][m][pre] = 0.;
				    if ( evento2.base[l][m][pre] != 31000. ) {
				        prmndp = prima - evento2.base[l][m][pre];
					if ( prmndp > 100000000. || prmndp < 100000000. ) prmndp = 0.;
					diffbas[l][m][pre] = (float)prmndp;
					calo->diffbas[l][m][pre] = (float)prmndp;
					if ( prmndp != 0. ) {
					    baseDIFFfull++;
					    //printf("differenza! %f %i %i %i %i %i \n",prmndp,i,l,m,pre,nn);
					    //return;
					};
				    };				    
				};			    
			    };
			};
		    };
		    //
		    // CALIBRATION ALGORITHM
		    // 
		    basel = evento.base[l][m][pre];		    
		    ener = evento.dexyc[l][m][n];
		    estrip[l][m][n] = 0.;
		    if ( ener > esup[l][m][n] && ener < 32000 ) esup[l][m][n] = ener;
		    if ( ener < einf[l][m][n] && ener > 0 ) einf[l][m][n] = ener;		    
		    if ( basel>0 && basel < 30000. && ener > 0. ){
			estrip[l][m][n] = (ener - calib.calped[l][m][n] - basel)/mip[l][m][n] ;
			//
			// OK, now in estrip we have the energy deposit in MIP of all the strips for this event (at the end of loops of course)
			//
			//if (estrip[l][m][n]<0.1 ) printf(" %i %i %i here %f \n",l,m,n,estrip[l][m][n]);
			calo->good[l][m][n] = (int)calib.calgood[l][m][n];
			if ( estrip[l][m][n] > evento.emin && calib.calgood[l][m][n] == 0 ) {
			    qtot += estrip[l][m][n];
			    nstrip += 1.;
			};
		    };
		    calib.sbase[l][m][pre] = evento.base[l][m][pre];
		};		
	    };
	};
	//
	// per ogni evento...
	//
	calo->qtot = qtot;	
	calo->nstrip = nstrip;	
	calo->evno = evno;
	calo->nobase = tot2;
	memcpy(calo->stwerr, de->stwerr, sizeof(stwerr));
	memcpy(calo->perror, de->perror, sizeof(perror));
	memcpy(calo->calevnum, de->calevnum, sizeof(calevnum));
	memcpy(calo->estrip, estrip, sizeof(estrip));
	tree->Fill();
	memcpy(estrip, estripnull, sizeof(estrip));
    };
    hfile->Write();
    hfile->Close();
    printf("\n");
    gSystem->ChangeDirectory(startingdir);
    return(0);
}