//
//   Calorimeter useful functions and subroutines - Emiliano Mocchiutti
//
//   CaloFunctions.h      version 4.02  (2006-01-24)
//
//   Programs in this file are called by other programs and cannot be run by hand.
//
//   Changelog:
//
//   4.01 - 4.02 (2006-01-24): small bug in getfileLEVname fixed.
//
//   4.00 - 4.01 (2005-12-13): bug in CaloFindCalibs (do not search on other calibration files when fails to find any calibration): fixed.
//
//   3.10 - 4.00 (2005-11-29): preparing for the final release, changed checkemilib (PAM_YODALIB -> PAM_YODA).
//
//   3.09 - 3.10 (2005-11-15): changed langaupro since it did not work in the compiled version; changed WhatToDo to work in compiled version.
//
//   3.08 - 3.09 (2005-11-11): small changes to PrintFigure subroutine
//
//   3.07 - 3.08 (2005-10-13): classes are now in another file (caloclasses.h)
//
//   3.06 - 3.07 (2005-10-10): added subroutine ColorTOFMIP for TOF colors.
//
//   3.05 - 3.06 (2005-10-04): added subroutine getfileLEVname.
//
//   3.04 - 3.05 (2005-08-04): Fixed bugs when running on a 64 bit arch.
//
//   3.03 - 3.04 (2005-07-14): Bug in getLEVname, fixed.
//
//   3.02 - 3.03 (2005-07-13): Finally CaloFunctions.h can be compiled under ROOT! Added clones of routines in the yodautility file no needed anymore.
//
//   3.01 - 3.02 (2005-07-11): Some more tuning for the macro compilation.
//
//   3.00 - 3.01 (2005-07-06): changed getLEVname since it gave some problems when called more than once in a ROOT sesssion. Tracker structurs now placed in a tracker
//                             include file.
//
//   2.27 - 3.00 (2005-07-05): Started the tuning for the compilation of the macros. Changed order of functions, added include files for ROOT routines, fixed some
//                             implicit casting. Fixed redefinition of class constructors. Small change in the calib structure.
//
//   2.26 - 2.27 (2005-06-23): Changed class CalorimeterLevel2 and structures to add good, perr,swerr and crc infos.
//
//   2.25 - 2.26 (2005-06-16): Changed class CalorimeterLevel2 and structures to convert OBT, headc and evfile to integers.
//
//   2.24 - 2.25 (2005-06-07): Added subroutine checkifpretyodaversion to check if file is older than 050515_007 (since that file there is a new format in TMTC data).
//
//   2.23 - 2.24 (2005-06-06): Small changes in the struct "calib". Change type of CaloPede subroutine, now it is an integer and not void: return 3 if it doesn't find the 
//                             calibration file.
//
//   2.22 - 2.23 (2005-05-30): Added checkmerging and fillmerging subroutines to avoid processing twice the same event.
//
//   2.21 - 2.22 (2005-05-19): Speed-up and better behaviour of CaloFindcalibs. No need to change calls to this subroutine, all internal changes.
//
//   2.20 - 2.21 (2005-05-17): Changed CaloFindCalibs in order to check also the calibration file provided if any but first check the current file.
//
//   2.19 - 2.20 (2005-05-06): Added mip to structure calib.
//
//   2.18 - 2.19 (2005-05-05): Added calorimeter class CalorimeterLevel2 and structures CaLevel1 CaLevel2. Changed getLEVname subroutines (still compatible with older version).
//
//   2.17 - 2.18 (2005-04-19): Added the old subroutine Calo1stcalib with the name Calo1stcalibOLD (to be used in ShowEvent.c).
//
//   2.16 - 2.17 (2005-04-14): Added classes for calorimeter ADC to MIP conversion. Bug fixed in CalorimeterCalibration class. CaloFindBaseRawNC subroutine added to 
//                             calculate baseline from raw data without any compression. Added getEmiFile subroutine to find YODA files (equivalent of getFile in
//                             yodautility.c); sintax is similar to getFile but it doesn't crash after 1111 or 2222 files.
//
//   2.15 - 2.16 (2005-04-11): Bugs fixes in the CaloFindCalibs, fetchpreviousfile and whatnamewith subroutines. Consider the case of an OBT jump in data.
//
//   2.14 - 2.15 (2005-04-07): Changes in the CaloFindCalibs subroutine. Added CaloCompressData subroutine, some small bugs fixed.
//
//   2.13 - 2.14 (2005-04-06): Changes in the CaloFindCalibs subroutine. Added fetchpreviousfile and whatnamewith subroutines. Now
//                             the program looks for calorimeter calibrations also in previous files and associates them to data in
//                             a correct way.
//
//   2.12 - 2.13 (2005-03-17): Changes in the WhatToDo subroutine.
//
//   2.11 - 2.12 (2005-03-09): PrintFigure subroutine changed to accomplish the new "saveas" script style.
//
//   2.10 - 2.11 (2005-03-04): Changed structure CalorimeterCalibration to save also fit parameters.
//
//   2.09 - 2.10 (2005-03-02): Show calibration CRC errors (in CaloFindCalibs). Added langaupro subroutines.
//
//   2.08 - 2.09 (2005-02-28): Calorimeter calibration classes changed (required by the new CaloADC2MIP).
//
//   2.07 - 2.08 (2005-02-25): check the existence of calibration header file in CaloFindCalibs!
//
//   2.06 - 2.07 (2005-02-24): bug in CaloFindBase! DO NOT put to zero the ADC values if you can't determine the baseline!
//
//   2.05 - 2.06 (2005-02-24): CalorimeterLevel1 class must be cleared during initialization to avoid the so-called "gundam-bug".
//
//   2.04 - 2.05 (2005-02-24): Changed variable definition from C/C++ style to ROOT style (int->Int_t).
//
//   2.03 - 2.04 (2005-02-23): Mispelling in the WhatToDo subroutine, fixed.
//
//   2.02 - 2.03 (2005-02-18): Calo1stCalib call CaloPede also if we don't have all four the calorimeter calibrations.
//
//   2.01 - 2.02 (2005-02-10): "progressive number" instead of "event number" in WhatToDo subroutine.
//
//   2.00 - 2.01 (2005-02-07): Added langau fit subroutines and the PrintFigure subroutine (used in CaloPLANES mainly). Some bug fixed,
//                             added user interaction routine and getfilename and getlevname subroutines to retrieve file name and to give
//                             the level name correctly.
//
//   1.00 - 2.00 (2005-01-17): Cleanup of the code
//
#include <iostream>
#include <TSystem.h>
#include <TH1.h>
#include <TF1.h>
#include <TROOT.h>
#include <TStyle.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 <TKey.h>
#include <TClassTable.h>

typedef struct Chmerge {
    Int_t fromflag;
    Int_t scanflag;
    Int_t pscuOBT;
    Int_t pscucount;
    Int_t lastOBT[1000];
    Int_t lastcount[1000];
} chmerge;

typedef struct Evento {
      Int_t     iev;
      Int_t     stwerr[4];
      Float_t   emin;
      Float_t   perror[4];
      Float_t   dexy[2][22][96];
      Float_t   dexyc[2][22][96];
      Double_t  base[2][22][6];
      Float_t   calselftrig[4][7];
      Float_t   calIItrig[4];
      Float_t   calstriphit[4];
      Float_t   calDSPtaberr[4];
      Float_t   calevnum[4];
      Float_t   shift; 
} evento;    

typedef struct Calib {
    TString tablename;
    TString basepath;
    TString yodalev;
    char *db;
    Int_t DW;
    Int_t mysql;
    Int_t obtjump;    
    Int_t time[4][51];
    Int_t ttime[4][51];
    Int_t fcode[4][51];
    Int_t iev;
    Int_t cstwerr[4];
    Float_t ispaw;
    Float_t cperror[4];
    Float_t mip[2][22][96];
    Float_t calped[2][22][96];
    Float_t calgood[2][22][96];
    Float_t calthr[2][22][6];
    Float_t calrms[2][22][96];
    Float_t calbase[2][22][6];
    Float_t calvar[2][22][6];
    Float_t calpuls[2][22][96];
    Double_t sbase[2][22][6];
    Double_t al_p[5][2];
    Int_t trkchi2;
    Bool_t good2;
} calib;

void emicheckLib(){
    if (!TClassTable::GetDict("pamela::TmtcEvent")){
	char *pamyodalib = gSystem->ExpandPathName("$PAM_YODA");
	stringstream libYoda;
	libYoda.str("");
	libYoda << pamyodalib << "/lib/libyoda.so";
	gSystem->Load(libYoda.str().c_str());
  }
}

void emidigForFiles(TList& out, TSystemDirectory *tsd, string defin){
    const char *initialDir = gSystem->pwd();
    TList       *lnk        = tsd->GetListOfFiles();
    if (lnk==0) return;	
    
    TSystemFile *file       = (TSystemFile*)lnk->First();
    if (file->IsZombie()) {      
	lnk->Delete();
	gSystem->cd(initialDir);
	return;
    }
    
    string      *fileDes    = new string(file->GetName());
    unsigned int loc = 0;
    while(file){
	fileDes = new string(file->GetName());
	if(!((fileDes->compare(".") == 0) || (fileDes->compare("..") == 0))){
	    if (file->IsDirectory()){
		emidigForFiles(out, (TSystemDirectory*)file, defin);
	    } else {
		loc = fileDes->find(defin.c_str(), 0);
		if (loc != string::npos){
		    out.AddLast((TObject*)file);
		}
	    }
	    
	}
	file = (TSystemFile*)lnk->After(file);
    }
    gSystem->cd(initialDir);
}

void emimakeAllFriend(TTree* out, TList* input){
    TList           *keys;
    TKey            *key;
    TTree           *tr;
    string          *base;
    string          *className;
    unsigned int             loc;
    bool            cont = false;
    TSystemFile *tsf = (TSystemFile*)input->First();
    if (tsf->IsZombie()) return;
    while(tsf){
        base    = new string(tsf->GetTitle());
        base->append("/");
        base->append(tsf->GetName());
        TFile *tf = new TFile(base->c_str());
	if (tf->IsZombie()) {
	  tsf = (TSystemFile*)input->After(tsf);
	  continue;
	}
        keys = tf->GetListOfKeys();
        key = (TKey*)keys->First();
        while(key){
            className = new string(key->GetClassName());
            loc = className->find("Tree", 0);
            if(loc != string::npos){
            tr = (TTree*)key->ReadObj();
            if (tr->GetEntries() > 0){
                if(cont){
                    tr = (TTree*)key->ReadObj();
                    out->AddFriend(tr->GetName(), tf);
                } else {
                    out = ((TTree*)key->ReadObj());
                    cont = true;
                }
            }
            }
            key = (TKey*)keys->After(key);
        }

        tsf = (TSystemFile*)input->After(tsf);
    }
}

void emimakeAllChained(TChain& out, TList* input){
    TList           *keys;
    TKey            *key;
    TTree           *tr;
    string          *base;
    string          *className;
    unsigned int             loc;

    TSystemFile *tsf = (TSystemFile*)input->First();
    if (tsf->IsZombie()) return;

    while(tsf){
        base    = new string(tsf->GetTitle());
        base->append("/");
        base->append(tsf->GetName());
        TFile *tf = new TFile(base->c_str());
	if (tf->IsZombie()) {
	  tsf = (TSystemFile*)input->After(tsf);
	  continue;
	}	

        keys = tf->GetListOfKeys();
        key = (TKey*)keys->First();
        while(key){
            className = new string(key->GetClassName());
            loc = className->find("Tree", 0);
            if(loc != string::npos){
		tr = (TTree*)key->ReadObj();
		if (tr->GetEntries() > 0){
                    out.Add(base->c_str());
		}
            }
            key = (TKey*)keys->After(key);
        }
        tsf = (TSystemFile*)input->After(tsf);
    }
}

TFile* emigetFile(TString base, TString packetType, TString subType = "Event"){  
  TSystemDirectory *targetDir  = new TSystemDirectory("", base); 
  TList       *filesList       = new TList;
  TSystemFile *tsf = 0;
  string      *tmpString = 0;
  tmpString = new string(".");
  tmpString->append(packetType);
  tmpString->append(".");
  tmpString->append(subType);
  tmpString->append(".");  
  emidigForFiles(*filesList, targetDir, tmpString->c_str());
  if ( filesList->IsEmpty() ) {
      filesList->Delete();
      return 0;
  };
  targetDir  = new TSystemDirectory("", base);
  tsf = (TSystemFile*)filesList->First();
  tmpString    = new string(tsf->GetTitle());
  tmpString->append("/");
  tmpString->append(tsf->GetName());
  targetDir->Delete();
  filesList->Delete();
  return new TFile(tmpString->c_str());      
}

Int_t emigetLastNotZeroBin(TH1 *histo){
  Int_t  lastBin = 0;
  Int_t  tempBin = 0;
  Stat_t minVal  = 0;
  Stat_t tempVal = 0;
  Int_t  range = histo->GetNbinsX();

  tempBin = histo->GetMaximumBin();  
  tempVal = histo->GetBinContent(tempBin);
  minVal  = tempVal;
  lastBin = tempBin;  

  for (Int_t i = histo->GetMaximumBin()  + 1; i < range; i++){
    tempVal = histo->GetBinContent(i);
    if ((tempVal != 0) && (tempVal < minVal)) {
      minVal  = tempVal;
      lastBin = i;
    }
  }
  return (Int_t)(lastBin*1.10);
}

Int_t emigetFirstNotZeroBin(TH1 *histo){
  Int_t  lastBin = 0;
  Int_t  tempBin = 0;
  Stat_t minVal  = 0;
  Stat_t tempVal = 0;

  tempBin = histo->GetMaximumBin();  
  tempVal = histo->GetBinContent(tempBin);
  minVal = tempVal;
  lastBin = tempBin;  

  for (Int_t i = histo->GetMaximumBin() - 1; i > 0; i--){
    tempVal = histo->GetBinContent(i);
    if ((tempVal != 0) && (tempVal < minVal)) {
      minVal  = tempVal;
      lastBin = i;
    }
  }
  return (Int_t)(lastBin*0.90);
}

void stringcopy(TString& s1, const TString& s2, Int_t from=0, Int_t to=0){    
    if ( to == 0 ){
	Int_t t2length = s2.Length();
	s1 = "";
	to = t2length;
    };
    for (Int_t i = from; i<to; i++){
	s1.Append(s2[i],1);
    };
}

void stringappend(TString& s1, const TString& s2){    
    Int_t t2length = s2.Length();
    for (Int_t i = 0; i<t2length; i++){
	s1.Append(s2[i],1);
    };
}

TFile* getEmiFile(TString base, TString packetType, TString subType="Event"){
    //
    TString newfile;
    stringcopy(newfile,base);
    //
    const char *base2 = base;
    stringstream fil;
    fil.str("");
    fil << base2;
    Int_t posiz = fil.str().find("dw_");
    Int_t upper = 0;
    if ( posiz == -1 ) {
	posiz = fil.str().find("DW_");
	upper = 1;
    };        
    if ( posiz == -1 ) return 0;
    TString fullname;
    stringcopy(fullname,base,posiz,posiz+15);
    stringappend(fullname,".");
    stringappend(fullname,packetType);
    stringappend(fullname,".");
    stringappend(fullname,subType);
    stringappend(fullname,".root");
    //
    stringappend(newfile,"/");
    stringappend(newfile,packetType);
    stringappend(newfile,"/");
    stringappend(newfile,fullname);
    //
    if ( !gSystem->IsFileInIncludePath(newfile) ){    
	return 0;
    } else {
	return new TFile(newfile);
    };
}
        
void CaloCompressData(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
    Float_t minstrip = 100000.;
    Float_t rms = 0.;
    for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	if ( calib.calgood[l][m][e] == 0. &&  evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip &&  evento.dexyc[l][m][e] > 0.) {
	    minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
	    rms = calib.calthr[l][m][pre];
	};
    };
    //
    // compression
    //
    if ( minstrip < evento.base[l][m][pre] && minstrip != 100000.){
	for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	    if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.;
	};
    };
}

void CaloFindBase(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
    Float_t ominstrip = 100000.;
    Float_t orms = 0.;
    Float_t minstrip = 100000.;
    Float_t rms = 0.;
    evento.base[l][m][pre] = 0.;
    for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	if ( calib.calgood[l][m][e] == 0. &&  evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip &&  evento.dexyc[l][m][e] > 0.) {
	    ominstrip = minstrip;
	    minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
	    orms = rms;
	    rms = calib.calthr[l][m][pre];
	};
	if ( calib.calgood[l][m][e] == 0. &&  evento.dexyc[l][m][e]-calib.calped[l][m][e] < ominstrip &&  evento.dexyc[l][m][e]-calib.calped[l][m][e] >minstrip ) {
	    ominstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
	    orms = calib.calthr[l][m][pre];
	};
    };
    if ( minstrip != 100000. ) {
	if ( ominstrip != 10000.) {
	    minstrip = ominstrip;
	    rms = orms;
	};
	Float_t strip6s = 0.;
	for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	    if ( (evento.dexyc[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexyc[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
		strip6s += 1.;
		evento.base[l][m][pre] += (evento.dexyc[l][m][e] - calib.calped[l][m][e]);
	    };
	    //
	    // compression
	    //
	    if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.;
	};
	if ( strip6s >= 9. ){	    
	    Double_t arro = evento.base[l][m][pre]/strip6s;
	    Float_t deci = 1000.*((float)arro - float(int(arro)));		      
	    if ( deci < 500. ) {
		arro = double(int(arro));
	    } else {
		arro = 1. + double(int(arro));
	    };
	    evento.base[l][m][pre] = arro;
	} else {
	    evento.base[l][m][pre] = 31000.;
	};
    } else {
	evento.base[l][m][pre] = 31000.;
    };
}

void CaloFindBaseRaw(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
    Float_t minstrip = 100000.;
    Float_t rms = 0.;
    evento.base[l][m][pre] = 0.;
    for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	if ( calib.calgood[l][m][e] == 0. &&  evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip &&  evento.dexy[l][m][e] > 0.) {
	    minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e];
	    rms = calib.calthr[l][m][pre];
	};
    };
    if ( minstrip != 100000. ) {
	Float_t strip6s = 0.;
	for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	    if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
		strip6s += 1.;
		evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]);
	    };
	    //
	    //  compression
	    //
	    if ( abs((int)(evento.dexy[l][m][e]-calib.calped[l][m][e])) <= (minstrip+rms) ) {
		evento.dexyc[l][m][e] = 0.;
	    } else {
		evento.dexyc[l][m][e] = evento.dexy[l][m][e];
	    };
	};
	if ( strip6s >= 9. ){	    
	    Double_t arro = evento.base[l][m][pre]/strip6s;
	    Float_t deci = 1000.*((float)arro - float(int(arro)));		      
	    if ( deci < 500. ) {
		arro = double(int(arro));
	    } else {
		arro = 1. + double(int(arro));
	    };
	    evento.base[l][m][pre] = arro;
	} else {
	    evento.base[l][m][pre] = 31000.;
	    for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
		evento.dexyc[l][m][e] = evento.dexy[l][m][e];
	    };
	};
    } else {
	evento.base[l][m][pre] = 31000.;
    };
}

void CaloFindBaseRawNC(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
    Float_t minstrip = 100000.;
    Float_t rms = 0.;
    evento.base[l][m][pre] = 0.;
    for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	if ( calib.calgood[l][m][e] == 0. &&  evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip &&  evento.dexy[l][m][e] > 0.) {
	    minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e];
	    rms = calib.calthr[l][m][pre];
	};
    };
    if ( minstrip != 100000. ) {
	Float_t strip6s = 0.;
	for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
	    if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
		strip6s += 1.;
		evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]);
	    };
	    //
	    //  compression
	    //
	    evento.dexyc[l][m][e] = evento.dexy[l][m][e];	    
	};
	if ( strip6s >= 9. ){	    
	    Double_t arro = evento.base[l][m][pre]/strip6s;
	    Float_t deci = 1000.*((float)arro - float(int(arro)));		      
	    if ( deci < 500. ) {
		arro = double(int(arro));
	    } else {
		arro = 1. + double(int(arro));
	    };
	    evento.base[l][m][pre] = arro;
	} else {
	    evento.base[l][m][pre] = 31000.;
	    for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
		evento.dexyc[l][m][e] = evento.dexy[l][m][e];
	    };
	};
    } else {
	evento.base[l][m][pre] = 31000.;
    };
}

int CaloPede(TString filename, Int_t s, Int_t atime, Calib & calib){
    TFile *calheadFile = 0;
    calheadFile = getEmiFile(filename, "CalibCalPed", "Header");
    TFile *calcalibFile = 0;
    calcalibFile = getEmiFile(filename, "CalibCalPed","Event");    
    //
    const char *filen = 0;
    filen = filename;
    if ( !calheadFile || !calcalibFile ) {
	if ( calheadFile ) calheadFile->Close();
	if ( calcalibFile ) calcalibFile->Close();
	printf("\n\nERROR: not able to open file: %s \n",filen);
	return(3);
    };
    TTree *tr = (TTree*)calheadFile->Get("Pscu");
    tr->AddFriend("CalibCalPed", calcalibFile);
    pamela::CalibCalPedEvent *ce = 0;
    pamela::PscuHeader *cph = 0;
    pamela::EventHeader *ceh = 0;
    tr->SetBranchAddress("Header", &ceh);
    tr->SetBranchAddress("Event", &ce);
    Long64_t ncalibs    = tr->GetEntries();
    for (Int_t ci = 0; ci < ncalibs ; ci++){
	tr->GetEntry(ci);
	cph = ceh->GetPscuHeader(); 		
	if ( atime == cph->GetOrbitalTime()){
	    calib.iev = ce->iev;
	    if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
		for ( Int_t d=0 ; d<11 ;d++  ){
		    Int_t pre = -1;
		    for ( Int_t j=0; j<96 ;j++){
			if ( j%16 == 0 ) pre++;
			if ( s == 2 ){
			    calib.calped[0][2*d+1][j] = ce->calped[3][d][j];
			    calib.cstwerr[3] = ce->cstwerr[3];
			    calib.cperror[3] = ce->cperror[3];
			    calib.calgood[0][2*d+1][j] = ce->calgood[3][d][j];
			    calib.calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
			    calib.calrms[0][2*d+1][j] = ce->calrms[3][d][j];
			    calib.calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
			    calib.calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
			};
			if ( s == 3 ){
			    calib.calped[0][2*d][j] = ce->calped[1][d][j];
			    calib.cstwerr[1] = ce->cstwerr[1];
			    calib.cperror[1] = ce->cperror[1];
			    calib.calgood[0][2*d][j] = ce->calgood[1][d][j];
			    calib.calthr[0][2*d][pre] = ce->calthr[1][d][pre];
			    calib.calrms[0][2*d][j] = ce->calrms[1][d][j];
			    calib.calbase[0][2*d][pre] = ce->calbase[1][d][pre];
			    calib.calvar[0][2*d][pre] = ce->calvar[1][d][pre];
			};
			if ( s == 0 ){
			    calib.calped[1][2*d][j] = ce->calped[0][d][j];
			    calib.cstwerr[0] = ce->cstwerr[0];
			    calib.cperror[0] = ce->cperror[0];
			    calib.calgood[1][2*d][j] = ce->calgood[0][d][j];
			    calib.calthr[1][2*d][pre] = ce->calthr[0][d][pre];
			    calib.calrms[1][2*d][j] = ce->calrms[0][d][j];
			    calib.calbase[1][2*d][pre] = ce->calbase[0][d][pre];
			    calib.calvar[1][2*d][pre] = ce->calvar[0][d][pre];
			};
			if ( s == 1 ){
			    calib.calped[1][2*d+1][j] = ce->calped[2][d][j];
			    calib.cstwerr[2] = ce->cstwerr[2];
			    calib.cperror[2] = ce->cperror[2];
			    calib.calgood[1][2*d+1][j] = ce->calgood[2][d][j];
			    calib.calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
			    calib.calrms[1][2*d+1][j] = ce->calrms[2][d][j];
			    calib.calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
			    calib.calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
			};
		    };
		};
	    };
	};
    };    	   
    tr->Delete();
    if ( calheadFile ) calheadFile->Close();
    if ( calcalibFile ) calcalibFile->Close();
    return(0);
}
 
TString getFilename(const TString filename){
    const string fil = (const char*)filename;
    Int_t posiz = fil.find("dw_");
    if ( posiz == -1 ) posiz = fil.find("DW_");
    if ( posiz == -1 ) return 0;
    Int_t posiz2 = posiz+13;
    TString file2;
    stringcopy(file2,filename,posiz,posiz2);
    TString pdat(".dat");
    stringappend(file2,pdat);
    return file2;
}

Int_t fetchpreviousfile(TString ffile, Int_t s){
    const char *file = ffile;
    stringstream fil;
    fil.str("");
    fil << file ;
    Int_t posiz = fil.str().find("dw_");
    Int_t upper = 0;
    if ( posiz == -1 ) {
	posiz = fil.str().find("DW_");
	upper = 1;
    };
    if ( posiz == -1 ) return(1);
    //
    TString nomefile = getFilename(ffile);    
    const char *ufnome = nomefile;
    stringstream nyear;
    stringstream nmonth;
    stringstream nday;
    stringstream nfno;
    TString year;
    TString month;
    TString day;
    TString fno;
    stringcopy(year, ufnome, 3, 5);
    stringcopy(month, ufnome, 5, 7);
    stringcopy(day, ufnome, 7, 9);
    stringcopy(fno, ufnome, 10, 13);
    const char *cyear = year;
    const char *cmonth = month;
    const char *cday = day;
    const char *cfno = fno;
    Int_t jy = atoi(cyear);
    Int_t jm = atoi(cmonth);
    Int_t jd = atoi(cday);
    Int_t jn = atoi(cfno);
    Int_t inter = 0;
    char *zy;
    char *zm;
    char *zd;
    char *zn = ""; 
    char *zyn;
    Int_t jyn = 0; 
    Int_t goodthis = 0; 
    Int_t onegood = 0; 
    Long64_t ncalibs;
    TTree *tr;
    pamela::CalibCalPedEvent *ce = 0;   
    pamela::PscuHeader *cph = 0;
    pamela::EventHeader *ceh = 0; 
    if ( jn > 1 ) {
	jn--;
    } else {
	jn = 999;
	if ( jd > 1 ) {
	    jd--;
	} else {
	    if ( jm > 1 ) {
		jm--;
	    } else {
		jm = 12;
		if ( jy > 3 ) {
		    jy--;
		} else {
		    return(1);
		};
	    };
	    if ( jm==4 || jm==6 || jm==9 || jm==11 ){
		jd = 30;
	    } else {
		jd = 31;
	    };
	    if ( jm==2 ) jd = 29;
	};
    };
    Int_t gday = 0;
    Int_t retvar = 1;
    stringstream nyd;
    while ( jy>3 ) {
	while ( jm>0 ){
	    while ( jd > 0 ){
		while ( jn > 0  ){
		    retvar--;
		    nyear.str("");
		    nyear << jy;
		    if ( jy<10 ){
			zy = "0";
		    } else {
			zy = "";
		    };
		    nmonth.str("");
		    nmonth << jm;
		    if ( jm<10 ){
			zm = "0";
		    } else {
			zm = "";			
		    };
		    nday.str("");
		    nday << jd;
		    if ( jd<10 ){
			zd = "0";
		    } else {
			zd = "";
		    };
		    nfno.str("");
		    nfno << jn;
		    if ( jn>99 ) zn = ""; 
		    if ( jn<100 && jn>9 ) zn = "0";
		    if ( jn<10 ) zn = "00";
		    jyn = 0; 
		    goodthis = 0; 
		    onegood = 0; 
		    while ( jyn < 100 ){
		    begwhile:  
			nyd.str("");
			nyd << jyn;
			if ( jyn<10 ){
			    zyn = "0";
			} else {
			    zyn = "";
			};
			stringstream nfile;
			nfile.str("");
			if ( upper ){
			    nfile << "DW_" << zy;
			} else {
			    nfile << "dw_" << zy;
			};
			nfile << nyear.str().c_str();
			nfile << zm;
			nfile << nmonth.str().c_str();
			nfile << zd;
			nfile << nday.str().c_str() << "_";
			nfile << zn;
			nfile << nfno.str().c_str();
			nfile << zyn;
			nfile << nyd.str().c_str();
			//
			TString fgfile = "";
			stringcopy(fgfile,ffile,0,posiz);
			stringappend(fgfile,nfile.str().c_str());
			TFile *calheadFile = 0;
			calheadFile = getEmiFile(fgfile, "CalibCalPed", "Header");
			TFile *calcalibFile = 0;
			calcalibFile = getEmiFile(fgfile, "CalibCalPed","Event");    
			//
			if ( !calheadFile || !calcalibFile ) {
			    if ( calheadFile ) calheadFile->Close();
			    if ( calcalibFile ) calcalibFile->Close();
			    if ( jyn > 0 && onegood ) {
				goodthis = 1;
				jyn--;
				goto begwhile;
			    };			    
			    if ( jyn == 0 ) {
				onegood = 0;
				jyn = 100;
				goto salta;
			    };
			    if ( jyn > 0 && !onegood ) {
				onegood = 0;
				jyn = 100;
				goto salta;
			    };
			    goto salta;
			} else {
			    onegood = 1;
			    if ( !goodthis ) {
				calheadFile->Close();
				calcalibFile->Close();
				goto salta;
			    };
			};
			//Takes the tree of the header file
			tr = (TTree*)calheadFile->Get("Pscu");
			tr->AddFriend("CalibCalPed", calcalibFile);
			tr->SetBranchAddress("Header", &ceh);
			tr->SetBranchAddress("Event", &ce);		    
			ncalibs = tr->GetEntries();
			if ( ncalibs == 0 ) {
			    jyn = 100;
			    onegood = 0;
			    goodthis = 0;
			    calheadFile->Close();
			    calcalibFile->Close();
			    goto salta;
			};
			inter = 0;
			for (Int_t c = 0; c < ncalibs; c++){
			    tr->GetEntry(c);
			    cph = ceh->GetPscuHeader(); 
			    if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
				inter++;
			    };
			};
			if ( inter ){
			    jyn = 100;
			    calheadFile->Close();
			    calcalibFile->Close();
			    return(retvar);
			};
			calheadFile->Close();
			calcalibFile->Close();
		    salta:
			jyn++;
		    };
		    jn--;
		};
		jn = 999;
		jd--;
		gday++;
		if ( gday>1 ) return(1);
	    };
	    jm--;
	    if ( jm==4 || jm==6 || jm==9 || jm==11 ){
		jd = 30;
	    } else {
		jd = 31;
	    };
	    if ( jm==2 ) jd = 29;
	};	    
	jm = 12;
	jy--;
    };
    return(retvar);
}

TString whatnamewith(TString ffile, Int_t retvar){
    Int_t debug = 0;
    const char *file = ffile; 
    stringstream fil;
    fil.str("");
    fil << file ;
    if ( debug ) printf("ffile %s retvar %i \n",file,retvar);
    Int_t posiz = fil.str().find("dw_");
    Int_t upper = 0;
    if ( posiz == -1 ) {
	posiz = fil.str().find("DW_");
	upper = 1;
    };
    if ( posiz == -1 ) return("wrong");
    //
    //    printf("posiz %i \n",posiz);
    TString nomefile = getFilename(ffile);   
    const char *ufnome = nomefile;
    stringstream nyear;
    stringstream nmonth;
    stringstream nday;
    stringstream nfno;
    TString year;
    TString month;
    TString day;
    TString fno;
    stringcopy(year, ufnome, 3, 5);
    stringcopy(month, ufnome, 5, 7);
    stringcopy(day, ufnome, 7, 9);
    stringcopy(fno, ufnome, 10, 13);
    const char *cyear = year;
    const char *cmonth = month;
    const char *cday = day;
    const char *cfno = fno;
    Int_t jy = atoi(cyear);
    Int_t jm = atoi(cmonth);
    Int_t jd = atoi(cday);
    Int_t jn = atoi(cfno);
    Int_t jyn = 0; 
    Int_t goodthis = 0; 
    Int_t onegood = 0; 
    char *zy;
    char *zm;
    char *zd;
    char *zn = ""; 
    char *zyn;
    Long64_t ncalibs;
    TString fgfile;
    TTree *tr;
    TFile *calheadFile = 0;
    TFile *calcalibFile = 0;
    if ( jn > 1 ) {
	jn--;
    } else {
	jn = 999;
	if ( jd > 1 ) {
	    jd--;
	} else {
	    if ( jm > 1 ) {
		jm--;
	    } else {
		jm = 12;
		if ( jy > 3 ) {
		    jy--;
		} else {
		    return("wrong");
		};
	    };
	    if ( jm==4 || jm==6 || jm==9 || jm==11 ){
		jd = 30;
	    } else {
		jd = 31;
	    };
	    if ( jm==2 ) jd = 29;
	};
    };
    Int_t retval = 1;
    Int_t gday = 0;
    stringstream nyd;
    while ( jy>3 ) {
	while ( jm>0 ){
	    while ( jd > 0 ){
		while ( jn > 0  ){
		    retval--;
		    if ( retval < retvar ) return("wrong");
		    nyear.str("");
		    nyear << jy;
		    if ( jy<10 ){
			zy = "0";
		    } else {
			zy = "";
		    };
		    nmonth.str("");
		    nmonth << jm;
		    if ( jm<10 ){
			zm = "0";
		    } else {
			zm = "";			
		    };
		    nday.str("");
		    nday << jd;
		    if ( jd<10 ){
			zd = "0";
		    } else {
			zd = "";
		    };
		    nfno.str("");
		    nfno << jn;
		    if ( jn>99 ) zn = ""; 
		    if ( jn<100 && jn>9 ) zn = "0";
		    if ( jn<10 ) zn = "00";
		    jyn = 0; 
		    goodthis = 0; 
		    onegood = 0; 
		    while ( jyn < 100 ){
		    begwhile:  
			nyd.str("");
			nyd << jyn;
			if ( jyn<10 ){
			    zyn = "0";
			} else {
			    zyn = "";
			};
			stringstream nfile;
			nfile.str("");
			if ( upper ){
			    nfile << "DW_" << zy;
			} else {
			    nfile << "dw_" << zy;
			};
			nfile << nyear.str().c_str();
			nfile << zm;
			nfile << nmonth.str().c_str();
			nfile << zd;
			nfile << nday.str().c_str() << "_";
			nfile << zn;
			nfile << nfno.str().c_str();
			nfile << zyn;
			nfile << nyd.str().c_str();
			//
			fgfile = "";
			stringcopy(fgfile,ffile,0,posiz);
			stringappend(fgfile,nfile.str().c_str());
			const char *tosee = fgfile;
			if ( debug ) printf("file: %s retval %i retvar %i \n",tosee,retval,retvar);
			calheadFile = getEmiFile(fgfile, "CalibCalPed", "Header");
			calcalibFile = getEmiFile(fgfile, "CalibCalPed","Event");    
			//
			if ( !calheadFile || !calcalibFile ) {
			    if ( calheadFile ) calheadFile->Close();
			    if ( calcalibFile ) calcalibFile->Close();
			    if ( jyn > 0 && onegood ) {
				goodthis = 1;
				jyn--;
				goto begwhile;
			    };			    
			    if ( jyn == 0 ) {
				onegood = 0;
				jyn = 100;
				goto salta;
			    };
			    if ( jyn > 0 && !onegood ) {
				onegood = 0;
				jyn = 100;
				goto salta;
			    };
			    goto salta;
			} else {
			    onegood = 1;
			    if ( !goodthis ){
				calheadFile->Close();
				calcalibFile->Close();
				goto salta;
			    };
			};
			tr = (TTree*)calheadFile->Get("Pscu");
			tr->AddFriend("CalibCalPed", calcalibFile);						
			ncalibs = tr->GetEntries();
			if ( ncalibs == 0 ) {
			    jyn = 100;
			    onegood = 0;
			    goodthis = 0;
			    calheadFile->Close();
			    calcalibFile->Close();
			    goto salta;
			};
			calheadFile->Close();
			calcalibFile->Close();
			if ( retval == retvar ) return(fgfile);	       
		    salta:
			jyn++;
		    };
		    jn--;
		};
		jn = 999;
		jd--;
		gday++;
		if ( gday>7 ) return("wrong");
	    };
	    jm--;
	    if ( jm==4 || jm==6 || jm==9 || jm==11 ){
		jd = 30;
	    } else {
		jd = 31;
	    };
	    if ( jm==2 ) jd = 29;
	};	    
	jm = 12;
	jy--;
    };
    return("wrong");
}

char *getLEVname(TString filename, TString detc, TString numb, TString frame="root"){
    //    char *file;
    stringstream file;
    file.str("");
    const char *num = numb;
    const char *det = detc;
    const char *fra = frame;
    string fil = (const char *)filename;
    Int_t posiz = fil.find("dw_");
    if ( posiz == -1 ) posiz = fil.find("DW_");
    if ( posiz == -1 ) return(0);
    Int_t spos = posiz;
    Int_t epos = posiz+15;
    TString tmptempf;
    stringcopy(tmptempf,filename,spos,epos);
    const char *tempf = tmptempf;
    file << tempf << ".Physics.Level";
    file << num << ".";
    file << det << ".Event.";
    file << fra;
    const char *rfile = file.str().c_str();
    return((char*)rfile);
};

char *getfileLEVname(TString filename, TString detc, TString numb, TString frame="root"){
    //    char *file;
    stringstream file;
    const char *num = numb;
    const char *det = detc;
    const char *fra = frame;
    string fil = (const char *)filename;
    Int_t posiz = fil.find("dw_");
    if ( posiz == -1 ) posiz = fil.find("DW_");
    if ( posiz == -1 ) return(0);
    Int_t spos = posiz;
    Int_t epos = posiz+13;
    TString tmptempf;
    stringcopy(tmptempf,filename,spos,epos);
    const char *tempf = tmptempf;
    file.str("");
    file << tempf << "00.Physics.Level";
    //    file << "00.Physics.Level";
    file << num << ".";
    file << det << ".Event.";
    file << fra;
    const char *rfile = file.str().c_str();
    return((char*)rfile);
};

void OLDCaloFindCalibs(TString &filename, Calib & calib){
    for (Int_t s = 0; s < 4; s++){
	for (Int_t d = 1; d<50; d++){
	    calib.ttime[s][d] = 0;
	    if ( d < 49 ) calib.time[s][d] = 0;
	};
    };
    TFile *calheadFile = 0;
    calheadFile = getEmiFile(filename, "CalibCalPed", "Header");
    TFile *calcalibFile = 0;
    calcalibFile = getEmiFile(filename, "CalibCalPed","Event");    
    //
    if ( !calheadFile || !calcalibFile ) {
	if ( calheadFile ) calheadFile->Close();
	if ( calcalibFile ) calcalibFile->Close();
	printf("No calibration header file! Exiting...\n");	
	return;
    };
    //Takes the tree of the header file
    TTree *tr = (TTree*)calheadFile->Get("Pscu");
    tr->AddFriend("CalibCalPed", calcalibFile);
    pamela::CalibCalPedEvent *ce = 0;   
    pamela::PscuHeader *cph = 0;
    pamela::EventHeader *ceh = 0; 
    tr->SetBranchAddress("Header", &ceh);
    tr->SetBranchAddress("Event", &ce);
    Long64_t ncalibs    = tr->GetEntries();
    Int_t inter;    	
    for (Int_t s = 0; s < 4; s++){
	for (Int_t d = 1; d<50; d++){
	    calib.ttime[s][d] = 0;
	    if ( d < 49 ) calib.time[s][d] = 0;
	};
    };
    for (Int_t s = 0; s < 4; s++){
	inter = 0;
	for (Int_t c = 0; c < ncalibs; c++){
	    tr->GetEntry(c);
	    cph = ceh->GetPscuHeader(); 
	    calib.ttime[s][inter] = 0;
	    if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
		calib.ttime[s][inter] = cph->GetOrbitalTime();
		inter++;
	    } else {
		if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) {
		    printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]);
		};
	    };
	};
	if ( inter == 0 ){
	    printf(" ERROR: no suitable calibration for section %i in this file!\n",s);
	};
	for (Int_t d = 1; d<50; d++){
	    if ( calib.ttime[s][d] != 0 ) {
		calib.time[s][d-1] = calib.ttime[s][d-1] + (int)((calib.ttime[s][d] - calib.ttime[s][d-1])/2.);
	    } else {
		if ( d == 1 ) {
		    calib.time[s][d-1] = 0;
		    };
	    };
	};
    };
    calheadFile->Close();
    calcalibFile->Close();
}

void CaloFindCalibs(TString &filename, TString &calibfile, Int_t &wused, Calib & calib){
    //
    // First of all check if the file has a monotone growing OBT time function. 
    //
    Int_t debug = 0;
    Int_t obtjump = 0;
    Int_t evjump[50];
    Int_t OBT = 0;
    Int_t OOBT = 0;
    Int_t OBT1 = 0;
    Int_t LOBT = 0;
    pamela::PscuHeader *ph = 0;
    pamela::EventHeader *eh = 0;
    TFile *headerF = 0;
    headerF = emigetFile(filename, "Physics", "Header");
    if ( !headerF ) {
	printf(" No physics file! \n ");
	return;
    };
    TTree *ctr = (TTree*)headerF->Get("Pscu");
    Long64_t nevents = ctr->GetEntries();
    if ( nevents == 0 ) {
	printf(" The file is empty! \n ");
	return;
    };
    for (Int_t i = 0; i < nevents; i++){
	    ctr->SetBranchAddress("Header", &eh);
	    ctr->GetEntry(i);
	    ph = eh->GetPscuHeader();
	    OBT = ph->GetOrbitalTime();
	    if ( !i ) OBT1 = OBT;
	    if ( i == nevents-1 ) LOBT = OBT;
	    if ( OOBT > OBT ) {
		evjump[obtjump] = i;
		obtjump++;
	    };
	    OOBT = OBT;
    };
    printf("\n Obtjump = %i - FIRST OBT %i - LAST OBT %i \n\n",obtjump,OBT1,LOBT);
    //
    TTree *tr;
    pamela::CalibCalPedEvent *ce = 0;   
    pamela::PscuHeader *cph = 0;
    pamela::EventHeader *ceh = 0; 
    Long64_t ncalibs;
    //
    for (Int_t s = 0; s < 4; s++){
	for (Int_t d = 1; d<50; d++){
	    calib.ttime[s][d] = 0;
	    if ( d < 49 ) calib.time[s][d] = 0;
	};
    };	
    for (Int_t s = 0; s < 4; s++){
	//
	Int_t firstlook = 0;
    repeatsearch:	
	//
	Int_t inter = 0;    	
	TFile *calheadFile = 0;
	calheadFile =  getEmiFile(filename, "CalibCalPed", "Header");
	TFile *calcalibFile = 0;
	calcalibFile =  getEmiFile(filename, "CalibCalPed","Event");    
	//
	if ( !calheadFile || !calcalibFile ) {
	    if ( calheadFile ) calheadFile->Close();
	    if ( calcalibFile ) calcalibFile->Close();
	    printf(" No calibration in this file! \n");
	    if ( !firstlook ){
		wused = 2;
		filename = calibfile;
		firstlook = 1;
		goto repeatsearch;
	    } else {
		wused = 0;
		goto jte;
	    };
	};
	wused = 1;
	//Takes the tree of the header file
	tr = (TTree*)calheadFile->Get("Pscu");
	tr->AddFriend("CalibCalPed", calcalibFile);
	tr->SetBranchAddress("Header", &ceh);
	tr->SetBranchAddress("Event", &ce);
	//
	ncalibs = tr->GetEntries();
	calib.obtjump = 0;
	inter = 0;
	for (Int_t c = 0; c < ncalibs; c++){
	    tr->GetEntry(c);
	    cph = ceh->GetPscuHeader(); 
	    calib.ttime[s][inter] = 0;
	    if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
		calib.ttime[s][inter] = cph->GetOrbitalTime();
		inter++;
	    } else {
		if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) {
		    printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]);
		};
	    };
	};
    jte: 
	if ( inter == 0 ){
	    printf(" WARNING: no suitable calibration for section %i in this file!\n",s);
	    printf(" I WILL SEARCH IN PREVIOUS FILES\n");
	    if ( !firstlook ){
		wused = 2;
		filename = calibfile;
		firstlook = 1;
		goto repeatsearch;
	    };
	};
	if ( inter > 50 ){
	    printf(" WARNING: cannot handle more than 50 calibrations for file!\n");
	    printf(" I WILL SEARCH IN PREVIOUS FILES\n");
	    inter = 0;
	    if ( !firstlook ){
		wused = 2;
		filename = calibfile;
		firstlook = 1;
		goto repeatsearch;
	    };
	};
	if ( obtjump ){
	    calib.obtjump = 1;
	    printf(" ERROR: jump in OBT! CPU restarded? \n");
	    printf(" WARNING: I will use only the first calibration (if any)\n");
	};
	//
	Int_t lastcal = 0;
	Int_t of = 0;
	if ( OBT1 < calib.ttime[s][0] || !inter ) {
	    //
	    // here we must look for a calibration in the previous file...
	    //
	    TString ffile = "";
	    stringcopy(ffile,filename);
	    Int_t fcodep = fetchpreviousfile(ffile,s);
	    if ( fcodep == 1 ) {
		if ( inter ){
		    of = -3;
		    printf(" WARNING: Problems to find a good calibration for section %i \n",s);
		} else {
		    of = -2;
		    printf(" ERROR: I was not able to find any good calibration for section %i \n",s);
		};
	    } else {
		stringcopy(ffile,filename);
		TString pfile = whatnamewith(ffile,fcodep);		    
		struct Calib tcal;
		OLDCaloFindCalibs(pfile,tcal);
		for ( Int_t d = 0; d<51; d++ ){
		    if ( tcal.ttime[s][d] == 0 ){ 
			lastcal = d-1; 
			break;
		    };
		};
		if ( !inter ){
		    calib.ttime[s][0] = tcal.ttime[s][lastcal];			
		    calib.time[s][0] = OBT1;
		    calib.time[s][1] = LOBT;
		    calib.fcode[s][0] = fcodep;
		    if ( debug ) printf("boh! %i \n",fcodep);
		    if ( debug ) getchar();
		    of = -1;
		} else {			
		    if ( obtjump ){
			if ( debug ) printf("boh2! %i \n",fcodep);
			if ( debug ) getchar();
			inter = 1;
			calib.time[s][0] = OBT1;
			calib.time[s][1] = LOBT;
			calib.fcode[s][0] = 10;
			of = -5;
			for (Int_t d = 1; d<50; d++){
			    calib.ttime[s][d] = 0;
			    if ( d>1 && d<49 ) calib.time[s][d] = 0;
			};	
		    } else {
			if ( lastcal >= 0 && tcal.ttime[s][lastcal]<OBT1 && (OBT1-tcal.ttime[s][lastcal])<7200000 && (OBT1-tcal.ttime[s][lastcal])>0 ){
			    for (Int_t i = inter+1; i>0; i--){
				calib.ttime[s][i] = calib.ttime[s][i-1];
			    };
			    calib.time[s][0] = OBT1;
			    calib.ttime[s][0] = tcal.ttime[s][lastcal];
			    calib.fcode[s][0] = fcodep;
			    of = 2;
			} else {
			    calib.time[s][0] = OBT1;
			    calib.fcode[s][0] = 10;
			    of = 1;
			};
		    };
		};
	    };
	} else {
	    of = 0;
	};
	if ( debug ) printf("boh3! %i \n",of);
	Int_t of1 = 0;
	if ( inter && of != 0 && of != 2 && of != -2 && of != -5) {
	    for (Int_t i = inter+1; i>0; i--){
		calib.ttime[s][i] = calib.ttime[s][i-1];
	    };
	    of1 = 0;
	};
	if ( of == -2 ) of = -1;
	if ( of == -3 ) {
	    of = 1;
	    calib.fcode[s][0] = 10;
	};
	//
	if ( of != -1 && of != -5 ){
	    if ( of == 2 ) of = 1;
	    for (Int_t d = 0; d<inter+1; d++){
		calib.time[s][d+of] = calib.ttime[s][d+of];
		calib.fcode[s][d+of] = 10;
		if ( d == inter ){
		    if ( debug ) printf("boh4! %i \n",of);
		    calib.time[s][d+of1+of] = LOBT;
		    calib.fcode[s][d+of] = 0;
		};		
	    };
	};
	if ( calheadFile ) calheadFile->Close();
	if ( calcalibFile ) calcalibFile->Close();
    };   
}

void Calo1stcalib(TString filename, Calib & calib, Int_t b[4]){
    Float_t estrip[2][22][96];
    //    
    // 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<50; d++){
	    calib.ttime[s][d] = 0 ;
	    if ( d < 49 ) calib.time[s][d] = 0 ;
	};
    };
    for (Int_t m = 0; m < 2 ; m++ ){
	for (Int_t k = 0; k < 22; k++ ){
	    for (Int_t l = 0; l < 96; l++ ){
		calib.calped[m][k][l] = 0. ;
		estrip[m][k][l] = 0.;
	    };
	};
    }
    //
    // first of all find the calibrations in the file 
    //
    OLDCaloFindCalibs(filename, calib);    
    //
    // print on the screen the results:   
    //
    printf(" ---------------------------------------------------------- \n \n");	
    Int_t calibex = 0;
    for (Int_t s=0; s<4;s++){
	Int_t stop = 0;
	for (Int_t d = 0; d<48; d++){
	    if ( calib.ttime[s][d] != 0 ) calibex++;
	    if ( calib.time[s][0] != 0 ){
		if ( d == 0 )  printf(" Section %i from time 0 to time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d]);
		if ( calib.time[s][d+1] != 0 ) {
		    printf(" Section %i from time %i to time %i use calibration at time %i \n",s,calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d+1]);
		} else {
		    if ( !stop ){
			printf(" Section %i from time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d+1]);		
			stop = 1;
		    };
		};
	    } else {
		if ( calib.ttime[s][d] != 0 ) printf(" Section %i from time 0 use calibration at time %i \n",s,calib.ttime[s][d]);
	    };
	};
	printf("\n");	
    };
    printf(" ---------------------------------------------------------- \n");	
    if ( calibex < 1 ) {
	printf("No full calibration data in this file, sorry!\n");
    } else {
	//
	// calibrate before starting
	//
	for (Int_t s = 0; s < 4; s++){
	    b[s]=0;
	    CaloPede(filename,s,calib.ttime[s][b[s]],calib);
	};	    
    };
}


void ColorMIP(Float_t mip, int& colo){
    //    printf("mip = %f \n",mip);
    if ( colo > 0 ){
	colo = 10;
	if ( mip > 0.7  ) colo = 38;
	if ( mip > 2.   ) colo = 4;    
	if ( mip > 10.  ) colo = 3;
	if ( mip > 100. ) colo = 2;
	if ( mip > 500. ) colo = 6;
    } else {
	colo = 10;
	if ( mip > 0.7  ) colo = 17;
	if ( mip > 2.   ) colo = 15;    
	if ( mip > 10.  ) colo = 14;
	if ( mip > 100. ) colo = 13;
	if ( mip > 500. ) colo = 12;
    };
}		

void ColorTOFMIP(Float_t mip, int& colo){
    //    printf("mip = %f \n",mip);
    if ( colo > 0 ){
	colo = 10;
	if ( mip > 0.  ) colo = 38;
	if ( mip > 2.   ) colo = 4;    
	if ( mip > 10.  ) colo = 3;
	if ( mip > 100. ) colo = 2;
	if ( mip > 500. ) colo = 6;
    } else {
	colo = 10;
	if ( mip > 0.  ) colo = 17;
	if ( mip > 2.   ) colo = 15;    
	if ( mip > 10.  ) colo = 14;
	if ( mip > 100. ) colo = 13;
	if ( mip > 500. ) colo = 12;
    };
}

Int_t WhatToDo(int& i, int& doflag, int& jumpto, Long64_t nevents, Float_t lastevno, Float_t firstevno, const char *file, TString outDir, TString figty, TCanvas &figure) {
    char *bw;
    if ( jumpto == -1 ){
	bw = "bw";
    } else {
	bw = "";
    };
    jumpto = 0;
    char tellme[256];
    stringstream input;
    char tellme2[256];
    stringstream input2;
    char tellme3[256];
    stringstream input3;
    stringstream out;
    out.str("x");
    input2.str("zzzzzzzzzzzzzz");
    input3.str("z");
    input << out.str().c_str();    
    stringstream stc;
    stringstream stc2;
    while ( !strcmp(input.str().c_str(),out.str().c_str()) ) {
	printf("\nPress <enter> to continue, b<enter> to go backward, j<enter> to jump to a\nparticular event number, p<enter> to save the figure, q<enter> to quit: \n");
	cin.getline(tellme,256);
	//
	input.str("");
	input << tellme;
	out << "y";
	//
	stc.str("b");
	//
	if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    if ( i > 0  ) {
		printf("WARNING: going backward!\n\n");
		doflag = 2;
	    } else {
		printf("This is the first event, you can't go backward! \n");
		out.str("");
		out << input.str().c_str();
	    };
	};
	stc.str("j");
	if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    printf("\n Do you want to jump to a progressive number [p] or to an event number [e]?  ");
	    cin.getline(tellme3,256);
	    input3.str("");
	    input3 << tellme3;
	    stc.str("p");
	    if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) {
		printf("\n Enter the progressive number you want to jump to: ");
		cin.getline(tellme2,256);
		input2.str("");
		input2 << tellme2;
		Int_t j;
		j = atoi(input2.str().c_str());
		if ( j < 1 || j > nevents+1 ) {
		    printf("\n You can choose between 1 and %i \n",(int)nevents+1);
		    out.str("");
		    out << input.str().c_str();
		} else { 
		    printf("\n Jumping to progressive number %i\n\n",j);
		    i = j-2;
		};
	    };
	    stc.str("e");
	    if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) {
		printf("\n Enter the event number you want to jump to: ");
		cin.getline(tellme2,256);
		input2.str("");
		input2 << tellme2;
		Int_t j;
		j = atoi(input2.str().c_str());
		if ( j < firstevno || j > lastevno ) {
		    printf("\n You can choose between %i and %i \n",(int)firstevno,(int)lastevno);
		    out.str("");
		    out << input.str().c_str();
		} else { 
		    printf("\n Jumping to event number %i\n\n",j);
		    jumpto = j;
		    i = -1;		    
		};
	    };
	    stc.str("e");
	    stc2.str("p");
	    if ( strcmp(input3.str().c_str(),stc.str().c_str()) && strcmp(input3.str().c_str(),stc2.str().c_str()) ) {
		printf(" You must type or \"p\" or \"e\"\n");
		out.str("");
		out << input.str().c_str();
	    };
	};
	//
	stc.str("q");
	if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    printf("Exiting...\n");
	    return(1);
	};	    
	//
	stc.str("p");
	if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n");
	    cin.getline(tellme2,256);
	    input2.str("");
	    input2 << tellme2;
	    out.str("");
	    out << input.str().c_str();
	    //
	    TString filename = file;
	    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 *ty = figty;
	    figsave.str("");
	    figsave << outdir << "/";
	    figsave << ty << "_";
	    figsave << (i+1) << "_";
	    figsave << figrec;
	    figsave << bw << ".";
	    figsave << input2.str().c_str();
	    figure.SaveAs(figsave.str().c_str());
	    printf("\n");		    
	};	    
    };	    
    return(0);
}

void PrintFigure(TString filename, TString outDir, TString figty, TString saveas, TCanvas& figure) {
    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;
    stringstream figsave;
    figsave.str("");
    figsave << outDir.Data() << "/";
    figsave << figrec << "_";
    figsave << figty.Data() << ".";
    figsave << saveas.Data();
    figure.SaveAs(figsave.str().c_str());
    return;
}

Double_t langaufun(Double_t *x, Double_t *par) {
      // Numeric constants
      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
      Double_t mpshift  = -0.22278298;       // Landau maximum location
      // Control constants
      Double_t np = 100.0;      // number of convolution steps
      Double_t sc =   3.;      // convolution extends to +-sc Gaussian sigmas
      // Variables
      Double_t xx;
      Double_t mpc;
      Double_t fland;
      Double_t sum = 0.0;
      Double_t xlow,xupp;
      Double_t step;
      Double_t i;
      // MP shift correction
      mpc = par[1] - mpshift * par[0]; 
      // Range of convolution integral
      xlow = x[0] - sc * par[3];
      xupp = x[0] + sc * par[3];
      step = (xupp-xlow) / np;
      // Convolution integral of Landau and Gaussian by sum
      for(i=1.0; i<=np/2; i++) {
         xx = xlow + (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);
         xx = xupp - (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);
      };
      return (par[2] * step * sum * invsq2pi / par[3]);
}



TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF, TString drw="QR0"){
    // Once again, here are the Landau * Gaussian parameters:
   //   par[0]=Width (scale) parameter of Landau density
   //   par[1]=Most Probable (MP, location) parameter of Landau density
   //   par[2]=Total area (integral -inf to inf, normalization constant)
   //   par[3]=Width (sigma) of convoluted Gaussian function
   //
   // Variables for langaufit call:
   //   his             histogram to fit
   //   fitrange[2]     lo and hi boundaries of fit range
   //   startvalues[4]  reasonable start values for the fit
   //   parlimitslo[4]  lower parameter limits
   //   parlimitshi[4]  upper parameter limits
   //   fitparams[4]    returns the final fit parameters
   //   fiterrors[4]    returns the final fit errors
   //   ChiSqr          returns the chi square
   //   NDF             returns ndf
   Int_t i;
   Char_t FunName[100];
   sprintf(FunName,"Fitfcn_%s",his->GetName());
   //
   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
   if (ffitold) delete ffitold;
   //
   TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
   ffit->SetParameters(startvalues);
   ffit->SetParNames("Width","MP","Area","GSigma");
   //
   for (i=0; i<4; i++) {
      ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
   };
   his->Fit(FunName,drw);   // fit within specified range, use ParLimits, do not plot
   //
   ffit->GetParameters(fitparams);    // obtain fit parameters
   for (i=0; i<4; i++) {
      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
   }
   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
   NDF[0] = ffit->GetNDF();           // obtain ndf
   //
   return (ffit);              // return fit function
}


//Double_t langaupro(Double_t *params, Double_t &maxx) {
Double_t langaupro(Double_t *params) {
   // Seaches for the location (x value) at the maximum of the 
   // Landau-Gaussian convolute and its full width at half-maximum.
   //
   Double_t p,x;
   Double_t step;
   Double_t l,lold;
   Int_t i = 0;
   Int_t MAXCALLS = 5000;
   // Search for maximum
   p = params[1] - 0.1 * params[0];
   step = 0.05 * params[0];
   lold = -2.0;
   l    = -1.0;
   //   while ( (l != lold) && (i < MAXCALLS) ) {
   while ( fabs(l-lold)>10e-12 && (i < MAXCALLS) ) {
      i++;
      lold = l;
      x = p + step;
      l = langaufun(&x,params);
      if ( l < lold ) step = -step/10.;
      p += step;
      //      printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step);
      //      printf(" l %.16g \n",l);
   };
   //   printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step);
   //printf(" l %.16g \n",l);
   if (i == MAXCALLS) return (-1.);
   //
   return(x);
}


void delay(int selftrig, int &del){
  int i;
  i = 0;
  del = 300;
  //
  while ((i < 16) && ((selftrig & (0x8000>>i)) == 0)) { 
    del += 50;
    i++;
  };
  if (i > 15) del = 1100;
}

Double_t fitraw(Double_t *x, Double_t *par){
    //
    Double_t fitval = par[3] + par[0]*TMath::Exp(-par[2]*(x[0]-par[1]));
    //
    return fitval;
}



int checkmerging(Chmerge & chmerge){
    for (Int_t i = chmerge.fromflag; i<1000; i++){
	if ( chmerge.pscuOBT == 0 && chmerge.pscucount == 0 ) return(1);
	if ( chmerge.pscuOBT == chmerge.lastOBT[i] && chmerge.pscucount == chmerge.lastcount[i] ){
	    chmerge.fromflag = i;
	    printf("\n WARNING! %i event number %i at OBT %i already processed, skipped!\n",i,chmerge.pscucount,chmerge.pscuOBT);
	    return(1);
	};
    };
    return(0);
}

int fillmerging(Chmerge & chmerge){
    chmerge.lastOBT[chmerge.scanflag] = chmerge.pscuOBT;
    chmerge.lastcount[chmerge.scanflag] = chmerge.pscucount;
    chmerge.scanflag++;
    if ( chmerge.scanflag > 999 ) return(0);
    chmerge.lastOBT[chmerge.scanflag] = 0;
    chmerge.lastcount[chmerge.scanflag] = 0;
    return(0);
}

//short int checkifpretyodaversion(TString filename){
//   string fil = (const char *)filename;
//  Int_t posiz = fil.find("dw_");
//  if ( posiz == -1 ) posiz = fil.find("DW_");
//  if ( posiz == -1 ) return (-1);    
//  TString year;
//  stringcopy(year,filename,posiz+3,posiz+5);   
//  const char *ye = year;
//  Int_t y = atoi(ye);
//  TString month;
//  stringcopy(month,filename,posiz+5,posiz+7);   
//  const char *mo = month;
//  Int_t m = atoi(mo);
//  TString day;
//  stringcopy(day,filename,posiz+7,posiz+9);   
//  const char *da = day;
//  Int_t d = atoi(da);
//  TString fno;
//  stringcopy(fno,filename,posiz+10,posiz+13);   
//  const char *fn = fno;
//  Int_t f = atoi(fn);
//  if ( y >= 5 && m >= 5 && d < 15 ) return(1); 
//  if ( y == 5 && m == 5 && d == 15 && f < 7 ) return(1); 
//  return(0);
//}


