//
//- Emiliano Mocchiutti
//
//   FCaloADC2MIP.c      version 1.02  (2006-03-31)
//
//   The only input needed is 
//
//   Changelog:
//
//   1.01 - 1.02 (2006-03-31): Added FCaloALIG2DAT and FCaloREADALIG to create and read the calorimeter alignment data file.
//
//   1.00 - 1.01 (2006-03-20): First flight version, read single yoda file.
//
//   0.00 - 1.00 (2006-03-20): Clone of CaloADC2MIP v4r01.
//
//
#include <fstream>
#include <sstream>
#include <string>
#include <iostream>
#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 <TGraph.h>
#include <TMultiGraph.h>
#include <TLine.h>
#include <TStyle.h>
#include <TF1.h>
//
#include <fcalostructs.h>
//
extern const char *pathtocalibration();
extern TString getFilename(const 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 CaloFindBaseRawNC(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);
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 <PamelaRun.h>
#include <physics/calorimeter/CalorimeterEvent.h>
#include <physics/trigger/TriggerEvent.h>
#include <CalibCalPedEvent.h>
//
#include <caloclassesfun.h>
#include <FCaloADC2MIPfun.h>
using namespace std;

void FCaloADC2MIP(TString Filename, TString calcalibfile = "", TString OutDir="", TString flist=""){
  const char *pam_calib = pathtocalibration();
  if ( !strcmp(OutDir.Data(),"") ){
    OutDir = pam_calib;
  };
  //
  TString calcalib;
  struct Evento evento;
  struct Calib calib;
  Int_t debug = 0;
  Int_t b[4];
  Int_t etime,evno;
  Float_t sdexy[2][22][96],sdexyc[2][22][96];
  Float_t  edelta[2][22][96];
  Float_t estrip[2][22][96];
  Float_t estripnull[2][22][96];
  Float_t delta;
  delta = 1.;
  evento.emin = 0.7;
  Int_t UPYES = 0;
  Int_t se =5;
  Int_t pre;
  Int_t rdone=0;
  Int_t done;
  bool isCOMP = 0;
  bool isFULL = 0;
  bool isRAW = 0;
  TH1F *pfmip[2][22][96];
  TF1 *pfitsnr[2][22][96];
  stringstream Fmip;
  for (Int_t i=0; i<2;i++){
    for (Int_t j=0; j<22;j++){
      for (Int_t m=0; m<96;m++){
	Fmip.str("");
	Fmip << "Fmip " << i;
	Fmip << " " << j;
	Fmip << " " << m;		   
	pfmip[i][j][m] = new TH1F(Fmip.str().c_str(),"",56,0.,55.);
	pfitsnr[i][j][m] = new TF1;
      };
    };	
  };       
  //
  //  the file which contains the calorimeter ADC to MIP conversion values is file2 = path to the script dir / CaloADC2MIP.root
  //
  stringstream file2;
  file2.str("");
  file2 << OutDir.Data() << "/CaloADC2MIP.root";
  //
  // this the figures and fit (4224 x 2 objects!)
  //
  stringstream file3;
  file3.str("");
  file3 << OutDir.Data() << "/CaloADC2MIPf.root";
  //
  stringstream file4;
  file4.str("");
  file4 << OutDir.Data() << "/CaloADC2MIPf";
  //
  stringstream file5;
  file5.str("");
  file5 << OutDir.Data() << "/CaloADC2MIPdata.root";
  //
  Int_t nfile = 0;
  Int_t fileup = 0;
  if ( flist != "" ){
    printf("\n You have given an input file list \n\n Files: \n");
    ifstream in;
    in.open(flist, ios::in);
    char ctime[20];
    while (1) {       
      in >> ctime;
      if (!in.good()) break;
      printf(" File number %i = %s \n",nfile+1,ctime);	    
      nfile++;
    };
    in.close();
    if ( !nfile ) {
      printf("\n\n The list is empty, exiting... \n\n");	    
      return;
    };
  } else {
    nfile = 1;
  };
  printf("\n");
  //
  // First check if the calibration file exists:
  //
  Int_t EFILE = 0;
  printf(" Check the existence of CaloADC2MIP.root file: \n\n");
  if ( !existfile(file2.str().c_str()) ){
    printf(" %s :no such file or directory \n",file2.str().c_str());
    printf("\n OK, I will create it!\n\n");
    EFILE = 0;
  } else {
    printf(" The file already exists! Check if an update is needed \n\n");
    EFILE = 1;
  };
  //
  // the check for an update is done once even if we have a list.
  //
  TTree *tree;
  CalorimeterCalibration *calo = new CalorimeterCalibration();
  if ( EFILE == 1) {
    //
    // The file exists, open it (READ ONLY) and check if any update is needed:
    //
    TFile *hfile = 0;
    hfile = new TFile(file2.str().c_str(),"READ","Calorimeter CALIBRATION data");
    tree = (TTree*)hfile->Get("CaloADC");  
    calo  = new CalorimeterCalibration();
    tree->SetBranchAddress("Event", &calo);    
    //
    Long64_t nevents = tree->GetEntries();
    printf("\n Number of updates: %i\n",(int)nevents);
    tree->GetEntry(nevents-1);
    //
    // Read the status variable of the last event: if it is 0 then we need an update, if it is one print out informations and exit
    //
    if ( calo->status == 0 ) {
      printf("\n An update is needed: \n");
      Int_t upstrip = 0;
      for (Int_t i = 0; i<2; i++){
	for (Int_t j = 0; j<22; j++){
	  for (Int_t l = 0; l<96; l++){
	    if ( calo->mask[i][j][l] == 0. ) upstrip++;
	  };
	};
      };
      printf(" %i strip out of 4224 have no statistic enough.\n\n",upstrip);
    } else {
      printf(" No update is needed! \n The data files used to obtain the ADC to MIP conversion are: \n");
      for (Int_t i = 0; i<nevents; i++){	
	tree->GetEntry(i);
	const char *fout = calo->fname;
	printf(" - %s \n",fout);
      };	
      printf("\n\n");
      return;
    };	
    hfile->Close();
  };
  //
  const char *filename2 = Filename;
  //
  TString lfile;
  TString sfile;
  //
  Int_t e = 0;
  ifstream in;
  in.open(flist, ios::in);
  char ctime[20];
  TString file2f;
  TString pfile;
  const char *cali;
  const char *ffile;    
  Int_t calibex = 0;
  //
  TString filename;
  const char *nome;
  TFile *hfile = 0;
  calo = new CalorimeterCalibration(); 
  Long64_t nevents;
  Int_t used = 0;
  const char *compare;
  Int_t hf3 = 0;
  Int_t ufile = 0;
  TFile *headerFile = 0;
  file2f = "";
  pfile = "";
  stringstream fmippa;
  const char *porca;
  Int_t S4;
  Int_t upnn;
  Int_t notgood = 0;
  Int_t notgood1 = 0;
  stringstream files;
  TFile *hfileb = 0;
  //
  Int_t pl1, pl2, strmin, strmax;
  Int_t spre=0;
  Int_t surstrip = 0;
  Float_t estrippa = 0.;
  Int_t surmam = 0;
  Int_t surman = 0;
  Int_t surmatrix[3][3];
  //
  TString fififile;
  const char *file;
  //
  stringstream nofi;
  nofi.str("");
  TFile *hfile3 = 0;
  //
  while ( e < nfile ){
    if ( flist != "" ){
      in >> ctime;
      nofi.str("");
      nofi << filename2 << "/";
      nofi << ctime;
      lfile = TString(nofi.str().c_str());
      fififile = getFilename(lfile);
      file = fififile;
      nofi.str("");
      nofi << file;
      sfile = TString(nofi.str().c_str());
    } else {
      nofi.str("");
      nofi << filename2;
      lfile = TString(nofi.str().c_str());
      fififile = getFilename(lfile);
      file = fififile;
      nofi.str("");
      nofi << file;
      sfile = TString(nofi.str().c_str());
    };
    //
    filename = lfile;	
    nome = filename;
    printf("\n Processing file number %i: %s \n\n",e,nome);
    nome = sfile;
    //
    // this must be done for each file.
    //
    const char *calname;
    calname = "";
    Int_t wused = 0;
    TTree *otr;
    pamela::calorimeter::CalorimeterEvent *de = 0;   
    pamela::trigger::TriggerEvent *trigger = 0;   
    pamela::PscuHeader *ph = 0;
    pamela::EventHeader *eh = 0;
    //
    if ( EFILE == 1 ){
      hfile = new TFile(file2.str().c_str(),"UPDATE","Calorimeter CALIBRATION data");	
      tree = (TTree*)hfile->Get("CaloADC");  
      calo = new CalorimeterCalibration();
      tree->SetBranchAddress("Event", &calo);		
      //
      // Check if the input file has already been used
      //	
      printf(" Check if file %s has already been used... \n",nome);
      nevents = tree->GetEntries();
      used = 0;
      for (Int_t i = 0; i<nevents; i++){	
	tree->GetEntry(i);
	compare= calo->fname;
	if ( !strcmp(compare,nome) ) {
	  used = 1;
	  break;
	};
      };		
      if ( used ) {
	printf(" The file %s has already been used!\n You cannot run twice this program with that file!\n\n",nome);
	hfile->Close();		    
	goto jumpfile;
      } else {
	printf(" Good, the file has not been used yet. \n\n");
	nevents    = tree->GetEntries() -1;
	tree->GetEntry(nevents);
      };
    } else {
      //
      // If the file doesn't exist create it
      //
      calo = new CalorimeterCalibration();
      hfile = new TFile(file2.str().c_str(),"NEW","Calorimeter CALIBRATION data");
      if ( !existfile(file2.str().c_str()) ){
	printf(" %s :no such file or directory \n",filename.Data());
	printf(" ERROR: Cannot create file \n");
	return;
      };
      //	    if ( hfile->IsZombie() ){
      //		printf(" ERROR: Cannot create file \n");
      //	return;
      //};
      tree = new TTree("CaloADC","Calorimeter calibration data");
      tree->Branch("Event","CalorimeterCalibration",&calo);	    
    };
    //
    //
    //
    hf3 = 0;
    if ( e == 0 ){
      if ( !existfile(file3.str().c_str()) ){
	ufile = 0;
      } else {
	printf(" Update figures file! \n\n");
	ufile = 1;
      };
    };
    if ( (EFILE == 1 || ufile == 1) && e == 0 ){
      printf(" Reading the old MIP figures to be updated...\n ");
      hfile3 = new TFile(file3.str().c_str());
      hf3 = 1;
      //
      stringstream fmino;
      fmino.str("");
      for (Int_t l = 0; l < 2; l++){
	for (Int_t m = 0; m < 22; m++){
	  for (Int_t n =0 ; n < 96; n++){		    
	    fmino.str("");
	    fmino << "fmip " << l;
	    fmino << " " << m;
	    fmino << " " << n;
	    gDirectory->Delete(fmino.str().c_str());
	    TH1F *temp  = (TH1F*)hfile3->Get(fmino.str().c_str());
	    pfmip[l][m][n] = (TH1F*)temp->Clone();			
	  };
	};
      };
      printf(" ...done!\n");
    };
    EFILE = 1;
    //
    // Check if the input calibration file contains any calibration data
    //
    done = 0;
    printf(" Check for calorimeter calibration: \n");
    calcalib = filename;
  findcalibs: printf(" Calibration file: %s \n",calname); 
    calname = calcalib;	
    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. ;
	  evento.dexy[m][k][l] = 0. ;
	  evento.dexyc[m][k][l] = 0. ;
	  edelta[m][k][l] = 0.;
	  estrip[m][k][l] = 0.;
	  estripnull[m][k][l] = 0.;
	};
      };
    }
    wused = 0;
    CaloFindCalibs(calcalib, calcalib, wused, calib);        
    //
    // print on the screen the results:   
    //
    ffile = filename;		
    printf(" ------ %s ------- \n \n",ffile);	
    calibex = 0;	
    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 ){
	    file2f = "";
	    stringcopy(file2f,calcalib,0,0);
	    pfile = (TString)whatnamewith(file2f,calib.fcode[s][d]);		    
	  } else {
	    pfile = (TString)calcalib;
	  };
	  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);
	};
      };
      printf("\n");	
    };
    printf(" ----------------------------------------------------------------------- \n \n");
    //
    if ( calibex < 1 ) {
      printf("No calibration data in this file!\n");
      calcalib = calcalibfile;
      if ( !done && calcalibfile != "" ) {
	cali = calcalibfile;
	printf(" Search in the calibration file %s \n",cali);
	done = 1;
	goto findcalibs;
      } else {
	printf(" Cannot find any calibration for this file! \n");		    
	goto jumpfile;
      };
    };
    if ( calibex < 4 ) {
      printf(" WARNING! No full calibration data in this file! \n\n");
    };		
    //
    // upgrade data...
    //    
    if ( !existfile(filename) ){
      printf(" %s :no such file or directory \n",filename.Data());
      goto jumpfile;
    };
    headerFile = new TFile(filename.Data());
    otr = (TTree*)headerFile->Get("Physics");
    otr->SetBranchAddress("Header", &eh);
    otr->SetBranchAddress("Calorimeter", &de);
    otr->SetBranchAddress("Trigger", &trigger);
    //
    nevents    = otr->GetEntries();
    //
    // calibrate before starting
    //
    for (Int_t s = 0; s < 4; s++){
      b[s]=0;
      if ( calib.fcode[s][b[s]] != 10 ){
	stringcopy(file2f,calcalib,0,0);
	pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);		    
      } else {
	pfile = (TString)calcalib;
      };
      porca = pfile;
      if ( debug ) printf(" porca miseria %s \n",porca);
      CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
    };
    //
    // run over all the events
    //
    printf("\n Processed events: \n\n");
    //
    for (Int_t i = 0; i < nevents; i++){
      evno = i;
      if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);	    
      //      printf(" %i \n",i);
      //
      // read from the header of the event the absolute time at which it was recorded
      //
      otr->GetEntry(i);
      //
      S4 = trigger->patterntrig[1] & (1<<0);
      //
      if ( !S4 || i == 0 ) {
	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 ){
		  stringcopy(file2f,calcalib,0,0);
		  pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);		    
		} else {
		  pfile = (TString)calcalib;
		};			 
		porca = pfile;
		if ( debug ) printf(" porca miseria %s \n",porca);
		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
	    //
	    if (l == 0 && m%2 == 0) se = 3;
	    if (l == 0 && m%2 != 0) se = 2;
	    if (l == 1 && m%2 == 0) se = 0;
	    if (l == 1 && m%2 != 0) se = 1;
	    //
	    // 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
	    //			
	    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++;
		rdone = 0;
		done = 0;
	      };	
	      if ( calo->mask[l][m][n] == 0. && ( isRAW || ( de->dexyc[l][m][n] > 0. &&  de->dexyc[l][m][n] < 32000.)) ){			    
		// 
		// baseline check and calculation
		//
		if ( isRAW ) {
		  if ( !rdone ){
		    CaloFindBaseRaw(calib,evento,l,m,pre);
		    rdone = 1;			
		  };
		} else {
		  if ( !done ){
		    evento.base[l][m][pre] = de->base[l][m][pre] ;   
		    evento.dexyc[l][m][n] = de->dexyc[l][m][n] ; 
		  };				    
		};
		//
		// no suitable new baseline, use old ones and compress data!
		//
		if ( !done && (evento.base[l][m][pre] > 30000. || evento.base[l][m][pre] == 0.) ){
		  evento.base[l][m][pre] = calib.sbase[l][m][pre];
		  upnn = n+16;
		  if ( upnn > 96 ) upnn = 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;			
		};
		//
		// CALIBRATION ALGORITHM
		//
		pl1 = -1;
		pl2 = -1;
		pl1 = m - 1;
		pl2 = m + 1;				    
		if ( pl1 < 0 ) {
		  pl1 = pl2;
		  pl2 = m + 2;
		};
		if ( pl2 > 21 ) {
		  pl2 = pl1;
		  pl1 = m - 2;
		};
		strmin = n - 1;
		if ( strmin < 0 ) strmin = 0;
		strmax = n + 1;
		if ( strmax > 95 ) strmax = 95;
		surstrip = 0;
		surmam = 0;
		surman = 0;
		for (Int_t im = 0; im<3 ;im++){
		  for (Int_t jm = 0; jm<3 ;jm++){
		    surmatrix[im][jm] = 0;
		  };
		};
		for (Int_t spl = pl1; spl<pl2+1; spl++){
		  if ( spl != m ){
		    surmam = spl - pl1;
		    for (Int_t sst = strmin; sst <= strmax ; sst++){					    
		      surman = sst - strmin;
		      if ( sst < 16 ) spre = 0;
		      if ( 15 < sst && sst < 32 ) spre = 1;
		      if ( 31 < sst && sst < 48 ) spre = 2;
		      if ( 47 < sst && sst < 64 ) spre = 3;
		      if ( 63 < sst && sst < 80 ) spre = 4;
		      if ( 79 < sst ) spre = 5;
		      if ( isRAW ){
			estrippa = ( evento.dexy[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
		      } else {
			estrippa = ( evento.dexyc[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
		      };
		      if ( estrippa > 0.7 && estrippa < 10. && calib.calgood[l][spl][sst] == 0 ){
			surstrip++;
			surmatrix[surmam][surman]++;
			//printf("surstrip = %i estrip = %f \n",surstrip,estrip);
		      };
		    };
		  };				    
		};
		if ( (surmatrix[0][0] && surmatrix[2][0]) ) surstrip = 0;
		if ( (surmatrix[0][2] && surmatrix[2][2]) ) surstrip = 0;
		//
		//
		estrip[l][m][n] = ( evento.dexyc[l][m][n] - calib.calped[l][m][n] - evento.base[l][m][pre]) ;
		if ( estrip[l][m][n] >= 0. && estrip[l][m][n] <= 56. && surstrip >=2 ){
		  pfmip[l][m][n]->Fill(estrip[l][m][n]);					    
		  UPYES = 1;
		};								
		calib.sbase[l][m][pre] = evento.base[l][m][pre];				
	      };
	    };
	  };
	};
      };
    };	
    //
    // save filename
    //
    calo->fname = sfile;	
    tree->Fill();
    hfile->Write();
    fileup++;
    hfile->Close();	
    //
    // create a backup copy!
    //
	
    if ( nfile > 1 && ( e%10 == 0 || e == nfile-1 ) && e > 0 ){
      files.str("");
      files << file4.str().c_str() << (e+1);
      files << ".bak";
      printf(" Save the backup file with figures:\n %s \n ",files.str().c_str());
      hfileb = new TFile(files.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");	
      if ( !existfile(files.str().c_str()) ){
	printf(" ERROR: Cannot create file!!! \n");
	goto oltre;
      };
      TH1F *bfmip[2][22][96];
      stringstream bmippa;
      bmippa.str("");
      for (Int_t l = 0; l < 2; l++){
	for (Int_t m = 0; m < 22; m++){		    
	  for (Int_t n =0 ; n < 96; n++){
	    bmippa.str("");
	    bmippa << "bmip " << l;
	    bmippa << " " << m;
	    bmippa << " " << n;
	    bfmip[l][m][n] = new TH1F(bmippa.str().c_str(),"",56,0.,55.);
	    bfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(bmippa.str().c_str());
	    bfmip[l][m][n]->Write();
	  };
	};
      };
      hfileb->Close();	
    };
    //
    //
    //
    headerFile->Close();
    //	caloFile->Close();
    //triggerFile->Close();
    goto oltre;
  jumpfile: printf(" File not processed! \n ");
    hfile->Close();	
  oltre:
    //	gObjectTable->Print();
    e++;
  };
  //
  if ( !UPYES ) goto end;
  printf("\n OK, now I will fit the MIP distribution for each strip \n\n");
  //
  hfile = new TFile(file2.str().c_str(),"UPDATE","Calorimeter CALIBRATION data");	
  calo = new CalorimeterCalibration();   
  tree = (TTree*)hfile->Get("CaloADC");  
  tree->SetBranchAddress("Event", &calo);		
  nevents = tree->GetEntries();
  if ( (nevents-1-fileup) > 0 ){
    tree->GetEntry(nevents-1-fileup);
  };
  //
  // Fit all 4224 distributions and save parameters, fit, figures (all basically).
  //    
  notgood = 0;
  notgood1 = 0;
  //
  //
  for (Int_t l = 0; l < 2; l++){
    for (Int_t m = 0; m < 22; m++){
      for (Int_t n =0 ; n < 96; n++){		    
	Double_t SNRPeak = 0.;		    
	if ( calo->mask[l][m][n] == 0. ){		      
	  //
	  // Setting fit range and start values
	  //
	  Double_t fr[2];
	  Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
	  if ( SNRPeak > 16. && SNRPeak < 35. ){
	    fr[0] = (Float_t)SNRPeak - 7.;
	    sv[0]= fp[0]; 
	    sv[1]= fp[1]; 
	    sv[2]= fp[2]; 
	    sv[3]= fp[3];		    
	  } else {
	    fr[0] = 19.;
	    sv[0] = 2.8; 
	    sv[1] = 21.0;
	    sv[2] = 1000.0; 
	    sv[3] = 0.1;		    
	  };
	  Int_t doitagain = 0;
	fitting: printf("Fitting strip %i %i %i \n",l,m,n);
	  fr[1] = 60.;
	  pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
	  plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
	  Double_t chisqr;
	  Int_t    ndf;
	  //
	  // fit the figure pfmip and put fit results in pfitsnr
	  //
	  pfitsnr[l][m][n] = langaufit(pfmip[l][m][n],fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf,"QR0");
	  //
	  // search for the maximum of the distribution, it will be the conversion factor between ADC channels and the MIP
	  //
	  Double_t SNRPeak = langaupro(fp);
	  printf("\n Conversion factor: %f \n\n",SNRPeak);
	  //
	  //
	  //
	  if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){
	    printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);
	    doitagain++;
	    if ( chisqr > 100. ) {
	      fr[0] = fr[0] + 5.;
	      sv[0] = fp[0]; 
	      sv[1] = fp[1]; 
	      sv[2] = fp[2]; 
	      sv[3] = fp[3];		    
	      goto fitting;
	    };
	    if ( doitagain == 2 ) {
	      fr[0] = 19.;
	      sv[0] = 2.8; 
	      sv[1] = 21.0;
	      sv[2] = 1000.0; 
	      sv[3] = 0.1;
	      goto fitting;
	    };
	    if ( doitagain == 1 ) {
	      fr[0] = 22.;
	      sv[0] = 2.8; 
	      sv[1] = 25.0;
	      sv[2] = 1000.0; 
	      sv[3] = 0.1;
	      goto fitting;
	    };
	    if ( doitagain == 3 ) {
	      fr[0] = 12.;
	      sv[0] = 2.8; 
	      sv[1] = 15.0;
	      sv[2] = 1000.0; 
	      sv[3] = 0.1;
	      goto fitting;
	    };
	  };

	  //
	  // save data
	  //
	  calo->mip[l][m][n] = (float)SNRPeak;
	  calo->ermip[l][m][n] = (float)fpe[1];
	  for (Int_t a = 0; a < 4 ; a++){
	    calo->fp[a][l][m][n] = (float)fp[a];
	    calo->fpe[a][l][m][n] = (float)fpe[a];
	  };
	  calo->ndf[l][m][n] = (float)ndf;
	  calo->chi2[l][m][n] = (float)chisqr;
	  if ( calo->fp[1][l][m][n] > 16. && calo->fp[1][l][m][n] < 35. ) notgood1++;
	  //
	  if ( ndf!=0 && chisqr/ndf < 2. && fpe[1] < 0.15 && fp[1] > 15. && fp[1] < 40. && SNRPeak > 0. ) {	
	    printf(" THIS IS A GOOD FIT, SAVED! \n\n");
	    calo->mask[l][m][n] = 1.;			
	  } else {
	    notgood++;
	    calo->mask[l][m][n] = 0.;			
	  };
	  printf("Fitting done, strip %i %i %i \n\n\n",l,m,n);
	  //
	  // draw the figure and fit on the screen
	  //		   
	  //c1->Clear();
	  //c1->cd();
	  //pfmip[l][m][n]->DrawCopy();	    
	  //pfitsnr[l][m][n]->DrawCopy("lsame");
	  //c1->Modified();
	  //c1->Update();
	};
      };
    };
  };
  //
  calo->status = 0;
  if ( !notgood ) calo->status = 1;
  //
  // save fit data
  //
  tree->Fill();
  hfile->Write();
  hfile->Close();
  //
  // Save histograms and fit figures
  //
  TFile *hfile2;
  hfile2 = new TFile(file3.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");	
  if ( !existfile(file3.str().c_str()) ){
    printf(" ERROR: Cannot create file \n");
    return;
  };
  TH1F *sfmip[2][22][96];
  TF1 *sfitsnr[2][22][96];
  fmippa.str("");
  for (Int_t l = 0; l < 2; l++){
    for (Int_t m = 0; m < 22; m++){
      for (Int_t n =0 ; n < 96; n++){
	fmippa.str("");
	fmippa << "fmip " << l;
	fmippa << " " << m;
	fmippa << " " << n;
	sfmip[l][m][n] = new TH1F(fmippa.str().c_str(),"",56,0.,55.);
	sfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(fmippa.str().c_str());
	fmippa.str("");
	fmippa << "fit " << l;
	fmippa << " " << m;
	fmippa << " " << n;
	sfitsnr[l][m][n] = (TF1*)pfitsnr[l][m][n]->Clone(fmippa.str().c_str());
	sfmip[l][m][n]->Write();
	sfitsnr[l][m][n]->Write();
      };
    };
  };
  hfile2->Close();	
  if ( hf3 ) hfile3->Close();		
  //
  if ( notgood ) {
    printf("\n An update will be necessary:\n %i strip with too low statistic (%i of them gave a result anyway)\n\n",notgood,notgood1);
  };
  //
 end: printf("\n");
  //
}

void FCaloRAWADC2MIPPLOT(TString Filename, TString calcalibfile = "", TString OutDir="", TString flist=""){
  const char *pam_calib = pathtocalibration();
  if ( !strcmp(OutDir.Data(),"") ){
    OutDir = pam_calib;
  };
  //
  TString calcalib;
  struct Evento evento;
  struct Calib calib;
  Int_t debug = 0;
  Int_t b[4];
  Int_t etime,evno;
  Float_t sdexy[2][22][96],sdexyc[2][22][96];
  Float_t  edelta[2][22][96];
  Float_t estrip[2][22][96];
  Float_t estripnull[2][22][96];
  Float_t delta;
  CalorimeterCalibration *calo = new CalorimeterCalibration();   
  delta = 1.;
  evento.emin = 0.7;
  Int_t se =5;
  Int_t pre;
  Int_t rdone = 0;
  Int_t done;
  bool isCOMP = 0;
  bool isFULL = 0;
  bool isRAW = 0;
  TH1F *pfmip[2][22][96];
  TF1 *pfitsnr[2][22][96];
  stringstream Fmip;
  Fmip.str("");
  for (Int_t i=0; i<2;i++){
    for (Int_t j=0; j<22;j++){
      for (Int_t m=0; m<96;m++){
	Fmip.str("");
	Fmip << "Fmip " << i;
	Fmip << " " << j;
	Fmip << " " << m;
	pfmip[i][j][m] = new TH1F(Fmip.str().c_str(),"",70,-15.,55.);
	pfitsnr[i][j][m] = new TF1;
      };
    };	
  };       
  //
  //  the file which contains the calorimeter ADC to MIP conversion values is file2 = path to the script dir / CaloADC2MIP.root
  //
  stringstream file2;
  file2.str("");
  file2 << OutDir.Data() << "/CaloADC2MIP.raw";
  //
  // this the figures and fit (4224 x 2 objects!)
  //
  stringstream file3;
  file3 << OutDir.Data() << "/CaloADC2MIPfr.raw";
  //
  stringstream file4;
  file4.str("");
  file4 << OutDir.Data() << "/CaloADC2MIPfr";
  //
  stringstream file5;
  file5.str("");
  file5 << OutDir.Data() << "/CaloADC2MIPdata.raw";
  //
  Int_t nfile = 0;
  Int_t fileup = 0;
  if ( flist != "" ){
    printf("\n You have given an input file list \n\n Files: \n");
    ifstream in;
    in.open(flist, ios::in);
    char ctime[20];
    while (1) {       
      in >> ctime;
      if (!in.good()) break;
      printf(" File number %i = %s \n",nfile+1,ctime);	    
      nfile++;
    };
    in.close();
    if ( !nfile ) {
      printf("\n\n The list is empty, exiting... \n\n");	    
      return;
    };
  } else {
    nfile = 1;
  };
  printf("\n");
  //
  // First check if the calibration file exists:
  //
  Int_t EFILE = 0;
  printf(" Check the existence of CaloADC2MIP.raw file: \n\n");

  if ( !existfile(file2.str().c_str()) ){
    printf(" %s :no such file or directory \n",file2.str().c_str());
    printf("\n OK, I will create it!\n\n");
    EFILE = 0;
  } else {
    printf(" The file already exists! Check if an update is needed \n\n");
    EFILE = 1;
  };
  //
  // the check for an update is done once even if we have a list.
  //
  TTree *tree = 0;
  if ( EFILE == 1) {
    //
    // The file exists, open it (READ ONLY) and check if any update is needed:
    //
    TFile *hfile = 0;
    hfile = new TFile(file2.str().c_str(),"READ","Calorimeter CALIBRATION data");
    calo = new CalorimeterCalibration();   
    tree = (TTree*)hfile->Get("CaloADC");  
    tree->SetBranchAddress("Event", &calo);    
    //
    Long64_t nevents = tree->GetEntries();
    printf("\n Number of updates: %i\n",(int)nevents);
    tree->GetEntry(nevents-1);
    //
    // Read the status variable of the last event: if it is 0 then we need an update, if it is one print out informations and exit
    //
    if ( calo->status == 0 ) {
      printf("\n An update is needed: \n");
      Int_t upstrip = 0;
      for (Int_t i = 0; i<2; i++){
	for (Int_t j = 0; j<22; j++){
	  for (Int_t l = 0; l<96; l++){
	    if ( calo->mask[i][j][l] == 0. ) upstrip++;
	  };
	};
      };
      printf(" %i strip out of 4224 have no statistic enough.\n\n",upstrip);
    } else {
      printf(" No update is needed! \n The data files used to obtain the ADC to MIP conversion are: \n");
      for (Int_t i = 0; i<nevents; i++){	
	tree->GetEntry(i);
	const char *fnout = calo->fname;	       
	printf(" - %s \n",fnout);
      };	
      printf("\n\n");
      return;
    };	
    hfile->Close();
  };
  //
  const char *filename2 = Filename;
  //
  TString lfile;
  TString sfile;
  //
  Int_t e = 0;
  ifstream in;
  in.open(flist, ios::in);
  char ctime[20];
  TString pfile;
  const char *cali;
  const char *ffile;    
  Int_t calibex = 0;
  //
  TString filename;
  const char *nome = 0;
  TFile *hfile = 0;
  //  CalorimeterCalibration *calo = 0;   
  Long64_t nevents;
  Int_t used = 0;
  const char *compare = 0;
  Int_t hf3 = 0;
  TFile *headerFile = 0;
  TString file2f = "";
  pfile = "";
  const char *porca = 0;
  Int_t S4;
  Int_t upnn;
  stringstream files;
  TFile *hfileb = 0;
  //
  Int_t pl1, pl2, strmin, strmax;
  Int_t spre=0;
  Int_t surstrip = 0;
  Float_t estrippa = 0.;
  //
  TFile *hfile3 = 0;
  TString fififile;
  const char *file;
  //
  calo = new CalorimeterCalibration();   
  stringstream finome;
  finome.str("");
  while ( e < nfile ){
    if ( flist != "" ){
      in >> ctime;
      finome.str("");
      finome << filename2 << "/";
      finome << ctime;
      lfile = TString(finome.str().c_str());
      fififile = getFilename(lfile);
      file = fififile;
      finome.str("");
      finome << file;
      sfile = TString(finome.str().c_str());
    } else {
      finome.str("");
      finome << filename2;
      lfile = TString(finome.str().c_str());
      fififile = getFilename(lfile);
      file = fififile;
      finome.str("");
      finome << file;
      sfile = TString(finome.str().c_str());
    };
    //
    filename = lfile;	
    nome = filename;
    printf("\n Processing file number %i: %s \n\n",e,nome);
    nome = sfile;
    //
    // this must be done for each file.
    //
    TTree *otr;
    pamela::calorimeter::CalorimeterEvent *de = 0;   
    pamela::trigger::TriggerEvent *trigger = 0;   
    pamela::PscuHeader *ph = 0;
    pamela::EventHeader *eh = 0;
    const char *calname;
    Int_t wused = 0;
    //
    //    TTree *tree;
    if ( EFILE == 1 ){
      hfile = new TFile(file2.str().c_str(),"UPDATE","Calorimeter CALIBRATION data");	
      tree = (TTree*)hfile->Get("CaloADC");  
      tree->SetBranchAddress("Event", &calo);		
      //
      // Check if the input file has already been used
      //	
      printf(" Check if file %s has already been used... \n",nome);
      nevents = tree->GetEntries();
      used = 0;
      for (Int_t i = 0; i<nevents; i++){	
	tree->GetEntry(i);
	compare= calo->fname;
	if ( !strcmp(compare,nome) ) {
	  used = 1;
	  break;
	};
      };		
      if ( used ) {
	printf(" The file %s has already been used!\n You cannot run twice this program with that file!\n\n",nome);
	hfile->Close();		    
	goto jumpfile;
      } else {
	printf(" Good, the file has not been used yet. \n\n");
	nevents    = tree->GetEntries() -1;
	tree->GetEntry(nevents);
      };
    } else {
      //
      // If the file doesn't exist create it
      //
      calo = new CalorimeterCalibration();   
      hfile = new TFile(file2.str().c_str(),"NEW","Calorimeter CALIBRATION data");
      if ( !existfile(file2.str().c_str()) ){
	printf(" %s :no such file or directory \n",filename.Data());
	printf(" ERROR: Cannot create file \n");
	return;
      };
      tree = new TTree("CaloADC","Calorimeter calibration data");
      tree->Branch("Event","CalorimeterCalibration",&calo);	    
    };
    //
    //
    //
    hf3 = 0;
    if ( EFILE == 1 && e == 0 ){
      printf(" Reading the old MIP figures to be updated...\n ");
      hfile3 = new TFile(file3.str().c_str());
      hf3 = 1;
      //
      stringstream fmippa;
      fmippa.str("");
      for (Int_t l = 0; l < 2; l++){
	for (Int_t m = 0; m < 22; m++){
	  for (Int_t n =0 ; n < 96; n++){		    
	    fmippa.str("");
	    fmippa << "fmip " << l;		       
	    fmippa << " " << m;
	    fmippa << " " << n;
	    TH1F *temp  = (TH1F*)hfile3->Get(fmippa.str().c_str());
	    pfmip[l][m][n] = (TH1F*)temp->Clone();			
	  };
	};
      };
      printf(" ...done!\n");
    };
    EFILE = 1;
    //
    // Check if the input calibration file contains any calibration data
    //
    done = 0;
    printf(" Check for calorimeter calibration: \n");
    calcalib = filename;	
  findcalibs: calname = calcalib; 
    printf(" Calibration file: %s \n",calname);
    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. ;
	  evento.dexy[m][k][l] = 0. ;
	  evento.dexyc[m][k][l] = 0. ;
	  edelta[m][k][l] = 0.;
	  estrip[m][k][l] = 0.;
	  estripnull[m][k][l] = 0.;
	};
      };
    }
    wused = 0;
    CaloFindCalibs(calcalib, calcalib, wused, calib);   
    //
    // print on the screen the results:   
    //
    ffile = filename;		
    printf(" ------ %s ------- \n \n",ffile);	
    calibex = 0;
    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 ){
	    file2f = "";
	    stringcopy(file2f,calcalib,0,0);
	    pfile = (TString)whatnamewith(file2f,calib.fcode[s][d]);		    
	  } else {
	    pfile = (TString)calcalib;
	  };
	  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);
	};
      };
      printf("\n");	
    };
    printf(" ----------------------------------------------------------------------- \n \n");
    //
    if ( calibex < 1 ) {
      printf("No calibration data in this file!\n");
      calcalib = calcalibfile;
      if ( !done && calcalibfile != "" ) {
	cali = calcalibfile;
	printf(" Search in the calibration file %s \n",cali);
	done = 1;
	goto findcalibs;
      } else {
	printf(" Cannot find any calibration for this file! \n");		    
	goto jumpfile;
      };
    };
    if ( calibex < 4 ) {
      printf(" WARNING! No full calibration data in this file! \n\n");
    };		
    //
    // upgrade data...
    //    
    if ( !existfile(filename) ){
      printf(" %s :no such file or directory \n",filename.Data());
      goto jumpfile;
    };
    headerFile = new TFile(filename.Data());
    otr = (TTree*)headerFile->Get("Physics");
    otr->SetBranchAddress("Header", &eh);
    otr->SetBranchAddress("Calorimeter", &de);
    otr->SetBranchAddress("Trigger", &trigger);
    //
    //
    nevents    = otr->GetEntries();
    //
    // calibrate before starting
    //
    for (Int_t s = 0; s < 4; s++){
      b[s]=0;
      if ( calib.fcode[s][b[s]] != 10 ){
	stringcopy(file2f,calcalib,0,0);
	pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);		    
      } else {
	pfile = (TString)calcalib;
      };
      porca = pfile;
      if ( debug ) printf(" porca miseria %s \n",porca);
      CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
    };
    //
    // run over all the events
    //
    printf("\n Processed events: \n\n");
    //
    for (Int_t i = 0; i < nevents; i++){
      evno = 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);
      //
      S4 = trigger->patterntrig[1] & (1<<0);
      //
      if ( !S4 || i == 0 ) {
	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 ){
		  stringcopy(file2f,calcalib,0,0);
		  pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);		    
		} else {
		  pfile = (TString)calcalib;
		};			 
		porca = pfile;
		if ( debug ) printf(" porca miseria %s \n",porca);
		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
	    //
	    if (l == 0 && m%2 == 0) se = 3;
	    if (l == 0 && m%2 != 0) se = 2;
	    if (l == 1 && m%2 == 0) se = 0;
	    if (l == 1 && m%2 != 0) se = 1;
	    //
	    // 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
	    //			
	    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++;
		  rdone = 0;
		  done = 0;
		};	
		if ( calo->mask[l][m][n] == 0. && ( isRAW || ( de->dexyc[l][m][n] > 0. &&  de->dexyc[l][m][n] < 32000.)) ){			    
		  // 
		  // baseline check and calculation
		  //
		  if ( !rdone ){
		    CaloFindBaseRawNC(calib,evento,l,m,pre);
		    rdone = 1;			
		  };
		  //
		  // no suitable new baseline, use old ones and compress data!
		  //
		  if ( !done && (evento.base[l][m][pre] > 30000. || evento.base[l][m][pre] == 0.) ){
		    evento.base[l][m][pre] = calib.sbase[l][m][pre];
		    upnn = n+16;
		    if ( upnn > 96 ) upnn = 96;
		    for ( Int_t nn = n; nn<upnn; nn++ ){
		      evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
		    };
		    done = 1;			
		  };
		  //
		  // CALIBRATION ALGORITHM
		  //
		  pl1 = -1;
		  pl2 = -1;
		  pl1 = m - 1;
		  pl2 = m + 1;				    
		  if ( pl1 < 0 ) {
		    pl1 = pl2;
		    pl2 = m + 2;
		  };
		  if ( pl2 > 21 ) {
		    pl2 = pl1;
		    pl1 = m - 2;
		  };
		  strmin = n - 2;
		  if ( strmin < 0 ) strmin = 0;
		  strmax = n + 2;
		  if ( strmax > 95 ) strmax = 95;
		  surstrip = 0;
		  for (Int_t spl = pl1; spl<pl2+1; spl++){
		    if ( spl != m ){
		      for (Int_t sst = strmin; sst<strmax+1; sst++){
			if ( sst < 16 ) spre = 0;
			if ( 15 < sst && sst < 32 ) spre = 2;
			if ( 31 < sst && sst < 48 ) spre = 3;
			if ( 47 < sst && sst < 64 ) spre = 4;
			if ( 63 < sst && sst < 80 ) spre = 5;
			if ( 79 < sst ) spre = 6;
			if ( isRAW ){
			  estrippa = ( evento.dexy[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
			} else {
			  estrippa = ( evento.dexyc[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
			};
			if ( estrippa > 0.7 && estrippa < 10. && calib.calgood[l][spl][sst] == 0 ){
			  surstrip++;
			  //printf("surstrip = %i estrip = %f \n",surstrip,estrip);
			};
		      };
		    };
		  };
		  //
		  //
		  estrip[l][m][n] = ( evento.dexyc[l][m][n] - calib.calped[l][m][n] - evento.base[l][m][pre]) ;
		  pfmip[l][m][n]->Fill(estrip[l][m][n]);					    
		  calib.sbase[l][m][pre] = evento.base[l][m][pre];				
		};
	      };
	    };
	  };
	};	    
      };
    };	
    //
    // save filename
    //
    calo->fname = sfile;	
    tree->Fill();
    hfile->Write();
    fileup++;
    hfile->Close();	
    //
    // create a backup copy!
    //
    if ( nfile > 1 && ( e%10 == 0 || e == nfile-1 ) && e > 0 ){
      files.str("");
      files << file4.str().c_str() << (e+1);
      files << ".bak";
      printf(" Save the backup file with figures:\n %s \n ",files.str().c_str());
      hfileb = new TFile(files.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");	
      if ( !existfile(files.str().c_str()) ){
	printf(" ERROR: Cannot create file!!! \n");
	goto oltre;
      };
      TH1F *bfmip[2][22][96];
      stringstream bmippa;
      for (Int_t l = 0; l < 2; l++){
	for (Int_t m = 0; m < 22; m++){
	  for (Int_t n =0 ; n < 96; n++){
	    bmippa.str("");
	    bmippa << "bmip " << l;
	    bmippa << " " << m;
	    bmippa << " " << n;
	    bfmip[l][m][n] = new TH1F(bmippa.str().c_str(),"",70,-15.,55.);
	    bfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(bmippa.str().c_str());
	    bfmip[l][m][n]->Write();
	  };
	};
      };
      hfileb->Close();	
    };
    //
    //
    //
    headerFile->Close();
    goto oltre;
  jumpfile: printf(" File not processed! \n ");
    hfile->Close();	
  oltre:
    //	gObjectTable->Print();
    e++;
  };
  //
  // Save histograms and fit figures
  //
  TFile *hfile2;
  hfile2 = new TFile(file3.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");	
  if ( !existfile(file3.str().c_str()) ){
    printf(" ERROR: Cannot create file!!! \n");
    return;
  };
  TH1F *sfmip[2][22][96];    
  stringstream fmippa;
  for (Int_t l = 0; l < 2; l++){
    for (Int_t m = 0; m < 22; m++){
      for (Int_t n =0 ; n < 96; n++){
	fmippa.str("");
	fmippa << "fmip " << l;
	fmippa << " " << m;
	fmippa << " " << n;
	sfmip[l][m][n] = new TH1F(fmippa.str().c_str(),"",72,-15.,55.);
	sfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(fmippa.str().c_str());
	sfmip[l][m][n]->Write();
      };
    };
  };
  hfile2->Close();	
  if ( hf3 ) hfile3->Close();	
}

void FCaloRAWADC2MIPDATA(TString Filename, TString calcalibfile = "", TString OutDir="", TString flist=""){
  const char *pam_calib = pathtocalibration();
  if ( !strcmp(OutDir.Data(),"") ){
    OutDir = pam_calib;
  };
  TString calcalib;
  struct Evento evento;
  struct Calib calib;
  Int_t debug = 0;
  Int_t b[4];
  Int_t etime,evno;
  Float_t sdexy[2][22][96],sdexyc[2][22][96];
  Float_t  edelta[2][22][96];
  Float_t estrip[2][22][96];
  Float_t estripnull[2][22][96];
  Float_t delta;
  delta = 1.;
  evento.emin = 0.7;
  Int_t se =5;
  Int_t pre;
  Int_t rdone = 0;
  Int_t  done;
  bool isCOMP = 0;
  bool isFULL = 0;
  bool isRAW = 0;
  //
  //  the file which contains the calorimeter ADC to MIP conversion values is file2 = path to the script dir / CaloADC2MIP.root
  //
  stringstream file5;
  file5.str("");
  file5 << OutDir.Data() << "/CaloADC2MIPdata.raw";
  //
  Int_t nfile = 0;
  Int_t fileup = 0;
  if ( flist != "" ){
    printf("\n You have given an input file list \n\n Files: \n");
    ifstream in;
    in.open(flist, ios::in);
    char ctime[20];
    while (1) {       
      in >> ctime;
      if (!in.good()) break;
      printf(" File number %i = %s \n",nfile+1,ctime);	    
      nfile++;
    };
    in.close();
    if ( !nfile ) {
      printf("\n\n The list is empty, exiting... \n\n");	    
      return;
    };
  } else {
    nfile = 1;
  };
  printf("\n");
  //
  // First check if the calibration file exists:
  //
  //
  const char *filename2 = Filename;
  //
  TString lfile;
  TString sfile;
  //
  Int_t e = 0;
  ifstream in;
  in.open(flist, ios::in);
  char ctime[20];
  TString file2f;
  TString pfile;
  const char *cali;
  const char *ffile;    
  Int_t calibex = 0;
  //
  TString filename;
  const char *nome = 0;
  Long64_t nevents;
  TFile *headerFile = 0;
  file2f = "";
  pfile = "";
  const char *porca = 0;
  Int_t upnn;
  //
  //
  TFile *hfl1 = 0;
  hfl1 = new TFile(file5.str().c_str(),"RECREATE","Calorimeter LEVEL1 data");
  // 
  // class CalorimeterLevel1
  //
  CalorimeterADCRAW *clev1 = new CalorimeterADCRAW();   
  TTree *trl1 = 0;
  trl1 = new TTree("CaloADCRAW","PAMELA Level1 calorimeter data");
  trl1->Branch("Rawdata","CalorimeterADCRAW",&clev1);
  //
  Int_t glbc = 0;
  //
  TString fififile;
  const char *file;
  //
  stringstream lfillo;
  while ( e < nfile ){
    if ( flist != "" ){
      in >> ctime;
      lfillo.str("");
      lfillo << filename2 << "/";
      lfillo << ctime;
      lfile = TString(lfillo.str().c_str());
      fififile = getFilename(lfile);
      file = fififile;
      lfillo.str("");
      lfillo << file;
      sfile = TString(lfillo.str().c_str());
    } else {
      lfillo.str("");
      lfillo << filename2;
      lfile = TString(lfillo.str().c_str());
      fififile = getFilename(lfile);
      file = fififile;
      lfillo.str("");
      lfillo << file;
      sfile = TString(lfillo.str().c_str());
    };
    //
    filename = lfile;	
    nome = filename;
    printf("\n Processing file number %i: %s \n\n",e,nome);
    nome = sfile;
    //
    // Check if the input calibration file contains any calibration data
    //
    TTree *otr;
    pamela::calorimeter::CalorimeterEvent *de = 0;   
    pamela::PscuHeader *ph = 0;
    pamela::EventHeader *eh = 0;
    //
    done = 0;
    printf(" Check for calorimeter calibration: \n");
    calcalib = filename;
  findcalibs: const char *calname = calcalib;
    printf(" Calibration file: %s \n",calname);
    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. ;
	  evento.dexy[m][k][l] = 0. ;
	  evento.dexyc[m][k][l] = 0. ;
	  edelta[m][k][l] = 0.;
	  estrip[m][k][l] = 0.;
	  estripnull[m][k][l] = 0.;
	};
      };
    }
    Int_t wused = 0;
    CaloFindCalibs(filename, calcalibfile, wused, calib);        
    if ( wused == 1 ) calcalibfile = filename;
    //
    // print on the screen the results:   
    //
    ffile = filename;		
    printf(" ------ %s ------- \n \n",ffile);	
    calibex = 0;
    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 ){
	    file2f = "";
	    stringcopy(file2f,calcalib,0,0);
	    pfile = (TString)whatnamewith(file2f,calib.fcode[s][d]);		    
	  } else {
	    pfile = (TString)calcalib;
	  };
	  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);
	};
      };
      printf("\n");	
    };
    printf(" ----------------------------------------------------------------------- \n \n");
    //
    if ( calibex < 1 ) {
      printf("No calibration data in this file!\n");
      calcalib = calcalibfile;
      if ( !done && calcalibfile != "" ) {
	cali = calcalibfile;
	printf(" Search in the calibration file %s \n",cali);
	done = 1;
	goto findcalibs;
      } else {
	printf(" Cannot find any calibration for this file! \n");		    
	goto jumpfile;
      };
    };
    if ( calibex < 4 ) {
      printf(" WARNING! No full calibration data in this file! \n\n");
    };		
    //
    // upgrade data...
    //    
    if ( !existfile(filename) ){
      printf(" %s :no such file or directory \n",filename.Data());
      goto jumpfile;
    };
    headerFile = new TFile(filename.Data());
    otr = (TTree*)headerFile->Get("Physics");
    otr->SetBranchAddress("Header", &eh);
    otr->SetBranchAddress("Calorimeter", &de);
    //
    nevents    = otr->GetEntries();
    //
    // calibrate before starting
    //
    for (Int_t s = 0; s < 4; s++){
      b[s]=0;
      if ( calib.fcode[s][b[s]] != 10 ){
	stringcopy(file2f,calcalib,0,0);
	pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);		    
      } else {
	pfile = (TString)calcalib;
      };
      porca = pfile;
      if ( debug ) printf(" porca miseria %s \n",porca);
      CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
    };
    //
    // run over all the events
    //
    printf("\n Processed events: \n\n");
    //
    Int_t updata;
    for (Int_t i = 0; i < nevents; i++){
      evno = 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);
      //
      //
      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 ){
		stringcopy(file2f,calcalib,0,0);
		pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);		    
	      } else {
		pfile = (TString)calcalib;
	      };			 
	      porca = pfile;
	      if ( debug ) printf(" porca miseria %s \n",porca);
	      CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
	    };
	  };
	};
      };
      //
      // do whatever you want with data
      //
      evento.iev = de->iev;
      //
      // run over views and planes
      //
      updata = 0;
      glbc++;
      for (Int_t l = 0; l < 2; l++){		    
	for (Int_t m = 0; m < 22; m++){
	  //
	  // determine the section number
	  //
	  if (l == 0 && m%2 == 0) se = 3;
	  if (l == 0 && m%2 != 0) se = 2;
	  if (l == 1 && m%2 == 0) se = 0;
	  if (l == 1 && m%2 != 0) se = 1;
	  //
	  // 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
	  //			
	  pre = -1; 
	  if ( isRAW || isFULL ){
	    updata = 1;
	    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->dexy[l][m][nn] ; 
	      // !!!
	    };
	    //
	    // run over strips
	    //
	    pre = -1;
	    for (Int_t n =0 ; n < 96; n++){		    
	      if ( n%16 == 0 ) {
		pre++;
		rdone = 0;
		done = 0;
	      };	
	      if ( isRAW || ( de->dexyc[l][m][n] > 0. &&  de->dexyc[l][m][n] < 32000.) || isFULL ){
		if ( isFULL && !rdone ) rdone = 1;					
		// 
		// baseline check and calculation
		//				    
		if ( !rdone ){
		  CaloFindBaseRawNC(calib,evento,l,m,pre);
		  rdone = 1;			
		};
		//
		// no suitable new baseline, use old ones and compress data!
		//
		if ( !done && (evento.base[l][m][pre] > 30000. || evento.base[l][m][pre] == 0.) ){
		  evento.base[l][m][pre] = calib.sbase[l][m][pre];
		  upnn = n+16;
		  if ( upnn > 96 ) upnn = 96;
		  for ( Int_t nn = n; nn<upnn; nn++ ){
		    evento.dexyc[l][m][nn] = de->dexy[l][m][nn] ;
		  };
		  done = 1;			
		};
		//
		// CALIBRATION ALGORITHM
		//
		estrip[l][m][n] = ( evento.dexyc[l][m][n] - calib.calped[l][m][n] - evento.base[l][m][pre]) ;
		calib.sbase[l][m][pre] = evento.base[l][m][pre];				
		clev1->diffbas[l][m][pre] = evento.base[l][m][pre];
		clev1->estrip[l][m][n] = estrip[l][m][n];
	      };
	    };
	  };
	};
      };	    
      //		if ( updata ){
      clev1->evno = glbc;
      clev1->etime = etime;
      trl1->Fill();
      for (Int_t l =0 ; l < 2; l++){
	for (Int_t m =0 ; m < 22; m++){
	  pre = -1;
	  for (Int_t n =0 ; n < 96; n++){
	    //			    memcpy(estrip, estripnull, sizeof(estrip));
	    if ( n%16 == 0 ) pre++;
	    estrip[l][m][n] = 0.;
	    evento.base[l][m][pre] = 0.;
	    clev1->diffbas[l][m][pre] = evento.base[l][m][pre];
	    clev1->estrip[l][m][n] = estrip[l][m][n];
	  };
	};
      };
      //gObjectTable->Print();
      //};
    };
    //
    // save filename
    //
    fileup++;
    //
    headerFile->Close();
    goto oltre;
  jumpfile: printf(" File not processed! \n ");
  oltre:
    e++;
  };
  hfl1->Write();
  hfl1->Close();
}

void FCalo4224STATUS(TString filename = "",  TString style = ""){
  const char *pam_calib = pathtocalibration();
  //
  // this routine shows the calorimeter calibration status with a figure. You can choose style between "error" (or "") to show the error on the 
  // peak maximum determined and "status" to show which strips have converged to a certain value.
  //
  if ( filename == "" ){
    stringstream filenome;
    filenome.str("");
    filenome << pam_calib << "/CaloADC2MIP.root";
    filename = (TString)filenome.str().c_str();
  };
  if ( style == "" ) style = "error";
  const char *file2 = filename;
  TFile *hfile;
  hfile = new TFile(file2,"READ","Calorimeter CALIBRATION data");
  CalorimeterCalibration *calo = 0;
  TTree *tree;
  tree = (TTree*)hfile->Get("CaloADC");  
  tree->SetBranchAddress("Event", &calo);    
  //
  Long64_t nevents = tree->GetEntries();
  tree->GetEntry(nevents-1);
  //
  // Book the histograms:
  //
  //
  TH2F *Xview = new TH2F("x-view event ","",96,-0.5,95.5,22,-0.5,21.5);
  TH2F *Yview = new TH2F("y-view event ","",96,-0.5,95.5,22,-0.5,21.5);	
  //
  // figures:
  //
  TCanvas *figura2 = new TCanvas("Calorimeter:_strip_RMS", "Calorimeter:_strip_RMS", 750, 650);
  figura2->SetFillColor(10);
  figura2->Range(0,0,100,100);
  Int_t bgcolor = 10;
  TPad *pd1 = new TPad("pd1","This is pad1",0.02,0.05,0.88,0.49,bgcolor);
  TPad *pd2 = new TPad("pd2","This is pad2",0.02,0.51,0.88,0.95,bgcolor);
  TPad *palette = new TPad("palette","This is palette",0.90,0.05,0.98,0.90,bgcolor);
  figura2->cd();
  TLatex *t=new TLatex();
  t->SetTextFont(32);
  t->SetTextColor(1);
  t->SetTextSize(0.03);
  t->SetTextAlign(12);
  t->DrawLatex(90.,92.5,"error");
  pd1->Range(0,0,100,100);
  pd2->Range(0,0,100,100);
  palette->Range(0,0,101,100);
  pd1->SetTicks();
  pd2->SetTicks();
  pd1->Draw();
  pd2->Draw();
  palette->Draw();
  palette->cd();
  const char *style2 = style;
  if (  !strcmp(style2,"error") ){
    // palette
    TPaveLabel *box1 = new TPaveLabel(1,2,99,14,"no fit","");
    box1->SetTextFont(32);
    box1->SetTextColor(1);	
    box1->SetTextSize(0.25);	
    box1->SetFillColor(10);
    box1->Draw();
    TPaveLabel *box2 = new TPaveLabel(1,16,99,28,"< 0.15","");
    box2->SetTextFont(32);
    box2->SetTextColor(1);		
    box2->SetTextSize(0.25);	
    box2->SetFillColor(38);
    box2->Draw();
    TPaveLabel *box3 = new TPaveLabel(1,30,99,42,"0.15 - 0.25","");
    box3->SetTextFont(32);
    box3->SetTextColor(1);	
    box3->SetTextSize(0.2);	
    box3->SetFillColor(4);
    box3->Draw();
    TPaveLabel *box4 = new TPaveLabel(1,44,99,56,"0.25 - 0.55","");
    box4->SetTextFont(32);
    box4->SetTextColor(1);	
    box4->SetTextSize(0.2);	
    box4->SetFillColor(3);
    box4->Draw();
    TPaveLabel *box5 = new TPaveLabel(1,58,99,70,"0.55 - 1","");
    box5->SetTextFont(32);
    box5->SetTextColor(1);	
    box5->SetTextSize(0.25);	
    box5->SetFillColor(2);
    box5->Draw();
    TPaveLabel *box6 = new TPaveLabel(1,72,99,84,"1 - 5","");
    box6->SetTextFont(32);
    box6->SetTextColor(1);	
    box6->SetTextSize(0.25);	
    box6->SetFillColor(6);
    box6->Draw();
    TPaveLabel *box7 = new TPaveLabel(1,86,99,98,"> 5","");
    box7->SetTextFont(32);
    box7->SetTextColor(10);	
    box7->SetTextSize(0.25);	
    box7->SetFillColor(1);
    box7->Draw();
  } else {
    TPaveLabel *box1 = new TPaveLabel(1,2,99,14,"not done","");
    box1->SetTextFont(32);
    box1->SetTextColor(1);	
    box1->SetTextSize(0.25);	
    box1->SetFillColor(10);
    box1->Draw();
    TPaveLabel *box2 = new TPaveLabel(1,16,99,28,"done","");
    box2->SetTextFont(32);
    box2->SetTextColor(1);		
    box2->SetTextSize(0.25);	
    box2->SetFillColor(38);
    box2->Draw();
  };
  //
  pd1->cd();	   
  gStyle->SetOptStat(0);
  Xview->SetXTitle("strip");
  Xview->SetYTitle("X - plane");
  Xview->GetYaxis()->SetTitleOffset(0.5);
  Xview->SetFillColor(bgcolor);	    
  Xview->Fill(1.,1.,1.);	    
  Xview->Draw("box");
  pd1->Update();
  pd2->cd();
  gStyle->SetOptStat(0);
  Yview->SetXTitle("strip");
  Yview->SetYTitle("Y - plane");
  Yview->GetYaxis()->SetTitleOffset(0.5);	
  Yview->SetFillColor(bgcolor);	    
  Yview->Fill(1.,1.,1.);	    
  Yview->Draw("box");
  pd2->Update();
  //
  // run over views and planes
  //
  stringstream xview;
  stringstream yview;
  for (Int_t m = 0; m < 22; m++){
    for (Int_t l = 0; l < 2; l++){
      for (Int_t n = 0; n < 96; n++){		    		    
	xview.str("");
	xview << "x-view event " << n;
	xview << " " << m;
	xview << " " << l;		
	yview.str("");
	yview << "y-view event " << n;
	yview << " " << m;
	yview << " " << l;		
	gDirectory->Delete(xview.str().c_str());
	gDirectory->Delete(yview.str().c_str());
	TH2F *Xview = new TH2F(xview.str().c_str(),"",96,-0.5,95.5,22,-0.5,21.5);
	TH2F *Yview = new TH2F(yview.str().c_str(),"",96,-0.5,95.5,22,-0.5,21.5);			   
	if (  !strcmp(style2,"error") ){
	  if ( calo->fpe[1][l][m][n] < 0.15 ){
	    Xview->SetFillColor(38);
	    Yview->SetFillColor(38);
	  };
	  if ( calo->fpe[1][l][m][n] >= 0.15 ){
	    Xview->SetFillColor(4);
	    Yview->SetFillColor(4);
	  };
	  if ( calo->fpe[1][l][m][n] >= 0.25){
	    Xview->SetFillColor(3);
	    Yview->SetFillColor(3);
	  };
	  if ( calo->fpe[1][l][m][n] >= 0.55 ){
	    Xview->SetFillColor(2);
	    Yview->SetFillColor(2);
	  };
	  if ( calo->fpe[1][l][m][n] >= 1. ){
	    Xview->SetFillColor(6);
	    Yview->SetFillColor(6);
	  };
	  if ( calo->fpe[1][l][m][n] >= 5. ){
	    Xview->SetFillColor(1);
	    Yview->SetFillColor(1);
	  };
	  if ( calo->mip[l][m][n] < 15. || calo->mip[l][m][n] > 40. ){
	    Xview->SetFillColor(10);
	    Yview->SetFillColor(10);
	  };
	} else {
	  if ( calo->mask[l][m][n] == 1 ){
	    Xview->SetFillColor(38);
	    Yview->SetFillColor(38);
	  } else {
	    Xview->SetFillColor(10);
	    Yview->SetFillColor(10);
	  };
	};
	if ( l == 0 ) {
	  Xview->Fill(n,m,1.);
	  pd1->cd();				    
	  Xview->Draw("box same");
	};				    
	if ( l == 1 ) {
	  Yview->Fill(n,m,1.);
	  pd2->cd();
	  Yview->Draw("box same");
	};		    
      };
    };
  };
  figura2->cd();
  gStyle->SetOptDate(1);
  //
  t=new TLatex();
  t->SetTextFont(32);
  t->SetTextColor(8);
  t->SetTextAlign(12);
  t->SetTextSize(0.02);
  figura2->Update();
  pd1->Update();
  pd2->Update();
  //
}

void FCalo4224FIT(TString filename = "", TString OutDir="", TString filevalue = "", TString type=""){
  const char *pam_calib = pathtocalibration();
  if ( !strcmp(OutDir.Data(),"") ){
    OutDir=pam_calib;
  };
  //
  // this routine shows the 4224 figures with fit from the main program
  //
  const char* startingdir = gSystem->WorkingDirectory();
  stringstream filenome;
  filenome.str("");
  if ( filename == "" ){
    filenome << OutDir.Data() << "/CaloADC2MIPf.root";
    filename = (TString)filenome.str().c_str();
    filenome.str("");
  };
  if ( filevalue == "" ){
    filenome.str("");
    filenome << OutDir.Data() << "/CaloADC2MIP.root";
  } else {
    filenome << filevalue.Data();
  };
  // EM
  const char *file2;
  file2 = filenome.str().c_str();
  TFile *hfile = 0;
  hfile = new TFile(file2,"UPDATE","Calorimeter CALIBRATION data");
  CalorimeterCalibration *calo = 0;
  TTree *tree;
  tree = (TTree*)hfile->Get("CaloADC");  
  tree->SetBranchAddress("Event", &calo);    
  Long64_t nevents = tree->GetEntries();
  tree->GetEntry(nevents-1);
  // EM end
  //
  const char *file = filename;
  stringstream file3;
  file3.str("");
  //  file3 << OutDir.Data() << "/";
  file3 << file;
  TFile hfile3(file3.str().c_str());
  Int_t l = 0;
  TCanvas *c1;
  Int_t ty = 1;
  if ( type != "" ) ty = 0;
  c1 = new TCanvas("canvas","canvas",800,800);
  c1->Draw();
  gStyle->SetOptFit(111);
  gPad->SetLogy();
  gPad->SetTickx();
  gPad->SetTicky();
  gPad->SetGrid();
  Int_t j = -1;
  Int_t k = -1;
  Int_t h = -1;;
  Int_t changed = 0;
 beginning:    
  if ( j >= 0 ) l = j;
  stringstream fmippa;
  TF1 *tempf = 0;
  TH1F *temp;
  while ( l < 2){
    Int_t m = 0;
    if ( k >= 0 ) m = k;
    while ( m < 22 ){	    
      Int_t n =0 ;
      if ( h >= 0 ) n = h;
      while ( n < 96 ){		    
	c1->Clear();
	c1->cd();		
	fmippa.str("");
	fmippa << "fmip " << l;
	fmippa << " " << m;
	fmippa << " " << n;
	temp  = (TH1F*)hfile3.Get(fmippa.str().c_str());
	if ( ty ) {
	  fmippa.str("");
	  fmippa << "fit " << l;
	  fmippa << " " << m;
	  fmippa << " " << n;
	  tempf  = (TF1*)hfile3.Get(fmippa.str().c_str());
	};
	j = -1;
	k = -1;
	h = -1;
	temp->DrawCopy();	    
	if ( ty ) tempf->DrawCopy("lsame");	    
	c1->Update();
	//
	//
	//
	printf(" **** STRIP: %i %i %i **** \n",l,m,n);
	printf(" MIP: %f \n",calo->mip[l][m][n]);
	printf(" ERMIP: %f \n",calo->ermip[l][m][n]);
	printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);
	printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);
	printf(" CHI2: %f \n",calo->chi2[l][m][n]);
	printf(" NDF: %f \n",calo->ndf[l][m][n]);
	printf(" MASK: %f \n",calo->mask[l][m][n]);		
	// 
	// what to do?
	//
	char    tellme[256];
	char    tellme2[256];
	stringstream input;
	stringstream input2;
	stringstream out;
	stringstream stc;
	input.str("a");
	out.str("a");    
	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\nfigure, p<enter> to save the figure, f<enter> to do the fit, q<enter> to quit: \n");
	  cin.getline(tellme,256);
	  input.str("");
	  input << tellme;
	  out.str("1");
	  //
	  stc.str("f");
	  if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    printf("\n Do you want to give a number by hand, h<enter>, or do you want to try a better fit, f<enter>?");
	    cin.getline(tellme,256);
	    input.str("");
	    input << tellme;
	    out.str("f");
	    if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {			    
	      Double_t fr[2],vtemp;
	      Double_t sv[4], pllo[4], plhi[4], ffp[4], ffpe[4];
	      for ( Int_t fk = 0; fk < 4; fk++){
		sv[fk] = (int)calo->fp[fk][l][m][n];
		printf(" Give input parameter [fp[%i]=%f] \n",fk,sv[fk]);
		cin.getline(tellme2,256);
		input2.str("");
		input2 << tellme2;
		vtemp = (Double_t)atoi(input2.str().c_str());
		if ( vtemp != 0. ) sv[fk] = vtemp;				
		printf(" -> fp[%i] = %f \n\n",fk,sv[fk]);			    
	      };
	      printf(" Give lower limit of fitting range fr[0] \n");
	      cin.getline(tellme2,256);
	      input2.str("");
	      input2 << tellme2;				
	      fr[0] = (Double_t)atoi(input2.str().c_str());
	      printf(" -> fr[0] = %f \n\n",fr[0]);
	      printf(" Give upper limit of fitting range fr[1] \n");
	      cin.getline(tellme2,256);
	      input2.str("");
	      input2 << tellme2;	
	      fr[1] = (Double_t)atoi(input2.str().c_str());
	      printf(" -> fr[1] = %f \n",fr[1]);
	      pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
	      plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
	      Double_t chisqr;
	      Int_t    fndf;
	      //
	      // fit the figure pfmip and put fit results in pfitsnr
	      //
	      tempf = langaufit(temp,fr,sv,pllo,plhi,ffp,ffpe,&chisqr,&fndf,"R0");
	      Double_t SNRPeak = langaupro(ffp);
	      printf("\n Conversion factor: %f \n",SNRPeak);
	      tempf->DrawCopy("lsame");	    
	      c1->Update();
	      c1->Modified();
	      c1->Update();
	      printf(" Do you want to save this result y/n<enter> ?\n");
	      cin.getline(tellme,256);
	      input.str("");
	      input << tellme;
	      out.str("y");
	      if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {			
		calo->mip[l][m][n] = (float)SNRPeak;
		calo->ermip[l][m][n] = (float)ffpe[1];
		for (Int_t a = 0; a < 4 ; a++){
		  calo->fp[a][l][m][n] = (float)ffp[a];
		  calo->fpe[a][l][m][n] = (float)ffpe[a];
		};
		calo->ndf[l][m][n] = (float)fndf;
		calo->chi2[l][m][n] = (float)chisqr;
		//
		if ( fndf!=0 && (float)chisqr/(float)fndf < 2. && ffpe[1] < 0.15 && ffp[1] > 15. && ffp[1] < 40. ) {					  
		  calo->mask[l][m][n] = 1.;			
		} else {
		  calo->mask[l][m][n] = 0.;			
		};
		printf(" Saved! \n");
		changed = 1;
	      };
	    } else {
	      printf(" Give the MIP conversion value \n");
	      cin.getline(tellme2,256);
	      input2.str("");
	      input2 << tellme2;	
	      calo->mip[l][m][n] = (Float_t)atoi(input2.str().c_str());
	      printf(" -> mip = %f \n",calo->mip[l][m][n]);
	      printf(" Give the mask value \n");
	      cin.getline(tellme2,256);
	      input2.str("");
	      input2 << tellme2;	
	      calo->mask[l][m][n] = (Float_t)atoi(input2.str().c_str());
	      printf(" -> mask = %f \n",calo->mask[l][m][n]);
	      changed = 1;
	    };
	    out.str("");
	    out << input.str().c_str();
	  };
	  stc.str("b");
	  if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    if ( n > 0  ) {
	      printf("WARNING: going one figure backward!\n\n");
	      n -= 2;	   
	    } else {
	      printf("This is the first strip of plane, 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 Enter the view number you want to jump to (0 = x, 1 = y): ");
	    cin.getline(tellme2,256);
	    input2.str("");
	    input2 << tellme2;	
	    j = atoi(input2.str().c_str());
	    if ( j < 0 || j > 1 ) {
	      printf("\n You can choose between 0 and 1 \n");
	      out.str("");
	      out << input.str().c_str();			   
	    } else { 
	      printf("\n Enter the plane number you want to jump to (0 to 21): ");
	      cin.getline(tellme2,256);
	      input2.str("");
	      input2 << tellme2;
	      k = atoi(input2.str().c_str());
	      if ( k < 0 || k > 21 ) {
		printf("\n You can choose between 0 and 21 \n");
		out.str("");
		out << input.str().c_str();	
	      } else {
		printf("\n Enter the strip number you want to jump to (0 to 95): ");
		cin.getline(tellme2,256);
		input2.str("");
		input2 << tellme2;
		h = atoi(input2.str().c_str());
		if ( h < 0 || h > 95 ) {
		  printf("\n You can choose between 0 and 95 \n");
		  out.str("");
		  out << input.str().c_str();
		} else {				   
		  printf("\n Jumping to strip %i %i %i \n\n",j,k,h);				    
		  goto beginning;
		};
	      };
	    };
	  };
	  //
	  stc.str("q");
	  if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    printf("Exiting...\n");
	    goto end;
	  };	    
	  //
	  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();
	    stringstream figsave;
	    figsave.str("");
	    figsave << startingdir << "/mip2adc_";
	    figsave << l << "_";
	    figsave << m << "_";
	    figsave << n << ".";
	    figsave << input2.str().c_str();
	    c1->SaveAs(figsave.str().c_str());
	    printf("\n");		    
	  };	    
	};	    
	n++;
      };
      m++;
    };
    l++;
  };
 end: printf("\n");
  if ( changed ){
    tree->Fill();
    hfile->Write();
  };
  hfile->Close();	
}

void FCalo4224MIPVALUES(TString filevalue = "", TString type=""){
  const char* startingdir = gSystem->WorkingDirectory();
  const char *pam_calib = pathtocalibration();
  //
  // this routine shows the 4224 figures with fit from the main program
  //
  stringstream filenome;
  filenome.str("");
  if ( filevalue == "" ){
    filenome << pam_calib << "/CaloADC2MIP.root";
  } else {
    filenome << filevalue.Data();
  };
  // EM
  const char *file2;
  file2 = filenome.str().c_str();
  TFile *hfile = 0;
  hfile = new TFile(file2,"READ","Calorimeter CALIBRATION data");
  CalorimeterCalibration *calo = 0;
  TTree *tree;
  tree = (TTree*)hfile->Get("CaloADC");  
  tree->SetBranchAddress("Event", &calo);    
  Long64_t nevents = tree->GetEntries();
  tree->GetEntry(nevents-1);
  // EM end
  //
  TCanvas *c1;
  Int_t ty = 1;
  if ( type != "" ) ty = 0;
  c1 = new TCanvas("canvas","canvas",1200,800);
  c1->Draw();
  gStyle->SetOptFit(111);
  gPad->SetTickx();
  gPad->SetTicky();
  Int_t k = -1;
  Int_t changed = 0;
  Float_t mipf[1]={0};
  Float_t val[1]={0};
  TLatex *text=new TLatex();
  TLine *linea;
  stringstream testo;
 beginning:    
  Int_t m = 0;
  if ( k >= 0 ) m = k;
  while ( m < 22 ){	    
    Int_t l = 0;
    c1->Clear();
    c1->cd();		
    TMultiGraph *mg = new TMultiGraph();
    while ( l < 2){
      Int_t n = 0;
      while ( n < 96 ){
	val[0] = (float)n + 100.*(float)l ;
	mipf[0] = (float)calo->mip[l][m][n];
	TGraph *graph = new TGraph(1, val, mipf);
	graph->SetMarkerColor(2);
	graph->SetMarkerStyle(22+l);
	graph->SetMarkerSize(0.7);
	mg->Add(graph);
	printf(" **** STRIP: %i %i %i **** \n",l,m,n);
	printf(" MIP: %f \n",calo->mip[l][m][n]);
	printf(" ERMIP: %f \n",calo->ermip[l][m][n]);
	printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);
	printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);
	printf(" CHI2: %f \n",calo->chi2[l][m][n]);
	printf(" NDF: %f \n",calo->ndf[l][m][n]);
	printf(" MASK: %f \n",calo->mask[l][m][n]);		
	n++;
      };
      l++;
    };
    mg->SetMaximum(50.);
    mg->SetMinimum(0.);
    mg->Draw("AP");
    for (Int_t e = 0; e<2 ; e++){
      for (Int_t ep = 0; ep<96 ; ep++){
	if ( ep%16 == 0 ) {
	  linea = new TLine((float)ep+100.*(float)e,0.,(float)ep+100.*(float)e,50.);
	  linea->SetLineColor(15);
	  linea->SetLineStyle(2);
	  linea->SetLineWidth(1);
	  linea->Draw("lsame");
	  if ( ep > 0 ){
	    linea = new TLine((float)ep-1.+100.*(float)e,0.,(float)ep-1.+100.*(float)e,50.);
	    linea->SetLineColor(15);
	    linea->SetLineStyle(2);
	    linea->SetLineWidth(1);
	    linea->Draw("lsame");
	  };
	};
	if ( ep == 95 ){
	  linea = new TLine(95.+100.*(float)e,0.,95.+100.*(float)e,50.);
	  linea->SetLineColor(15);
	  linea->SetLineStyle(2);
	  linea->SetLineWidth(1);
	  linea->Draw("lsame");
	};
      };
    };
    text->SetTextAngle(0);
    text->SetTextFont(32);
    text->SetTextColor(1);
    text->SetTextSize(0.050); 
    text->SetTextAlign(12);
    testo.str("");
    testo << "Plane " << m;
    text->DrawLatex(85.,52.,testo.str().c_str());
    c1->Update();
    // 
    // what to do?
    //
	
    char    tellme[256];
    char    tellme2[256];
    stringstream input;
    stringstream input2;
    stringstream out;
    stringstream stc;
    input.str("a");
    out.str("a");    
    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\nfigure, p<enter> to save the figure, q<enter> to quit: \n");
      cin.getline(tellme,256);
      input.str("");
      input << tellme;
      out.str("1");
      //
      stc.str("b");
      if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	if ( m > 0  ) {
	  printf("WARNING: going one figure backward!\n\n");
	  m -= 2;	   
	} else {
	  printf("This is the first plane, 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 Enter the plane number you want to jump to (0 to 21): ");
	cin.getline(tellme2,256);
	input2.str("");
	input2 << tellme2;
	k = atoi(input2.str().c_str());
	if ( k < 0 || k > 21 ) {
	  printf("\n You can choose between 0 and 21 \n");
	  out.str("");
	  out << input.str().c_str();	
	} else {
	  printf("\n Jumping to plane %i \n\n",k);				    
	  m = k;
	  goto beginning;
	};
      };
      //
      stc.str("q");
      if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	printf("Exiting...\n");
	goto end;
      };	    
      //
      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();
	stringstream figsave;
	figsave.str("");
	figsave << startingdir << "/mip2adc_plane_";
	figsave << m << ".";
	figsave << input2.str().c_str();
	c1->SaveAs(figsave.str().c_str());
	printf("\n");		    
      };	    
    };	    
    m++;
  };
 end: printf("\n");
  if ( changed ){
    tree->Fill();
    hfile->Write();
  };
  hfile->Close();	
}

void FCalo4224BAK(TString filename){
  //
  // this shows the 4224 figure from a backup file
  //
  const char* startingdir = gSystem->WorkingDirectory();
  const char *file3 = filename;
  TFile hfile3(file3);
  Int_t l = 0;
  TCanvas *c1;
  c1 = new TCanvas("canvas","canvas",800,800);
  c1->Draw();
  gStyle->SetOptFit(111);
  gPad->SetLogy();
  Int_t j = 0;
  Int_t k = 0;
  Int_t h = 0;
 beginning:
  if ( j ) l = j;
  stringstream bmippa;
  while ( l < 2){
    Int_t m = 0;
    if ( k ) m = k;
    while ( m < 22 ){	    
      Int_t n =0 ;
      if ( h ) n = h;
      while ( n < 96 ){		    
	c1->Clear();
	c1->cd();		
	bmippa.str("");
	bmippa << "bmip " << l;
	bmippa << " " << m;
	bmippa << " " << n;
	TH1F *temp  = (TH1F*)hfile3.Get(bmippa.str().c_str());
	j = 0;
	k = 0;
	h = 0;
	temp->DrawCopy();	    
	c1->Update();
	// 
	// what to do?
	//
	char    tellme[256];
	char    tellme2[256];
	stringstream input;
	stringstream input2;
	stringstream out;
	stringstream stc;
	input.str("a");
	out.str("a");    
	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\nfigure, p<enter> to save the figure, q<enter> to quit: \n");
	  cin.getline(tellme,256);
	  input.str("");
	  input << tellme;
	  out.str("1");
	  //
	  stc.str("b");
	  if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    if ( n > 0  ) {
	      printf("WARNING: going one figure backward!\n\n");
	      n -= 2;	   
	    } else {
	      printf("This is the first strip of plane, 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 Enter the view number you want to jump to (0 = x, 1 = y): ");
	    cin.getline(tellme2,256);
	    input2.str("");
	    input2 << tellme2;	
	    j = atoi(input2.str().c_str());
	    if ( j < 0 || j > 1 ) {
	      printf("\n You can choose between 0 and 1 \n");
	      out.str("");
	      out << input.str().c_str();			   
	    } else { 
	      printf("\n Enter the plane number you want to jump to (0 to 21): ");
	      cin.getline(tellme2,256);
	      input2.str("");
	      input2 << tellme2;
	      k = atoi(input2.str().c_str());
	      if ( k < 0 || k > 21 ) {
		printf("\n You can choose between 0 and 21 \n");
		out.str("");
		out << input.str().c_str();	
	      } else {
		printf("\n Enter the strip number you want to jump to (0 to 95): ");
		cin.getline(tellme2,256);
		input2.str("");
		input2 << tellme2;
		h = atoi(input2.str().c_str());
		if ( h < 0 || h > 95 ) {
		  printf("\n You can choose between 0 and 95 \n");
		  out.str("");
		  out << input.str().c_str();
		} else {				   
		  printf("\n Jumping to strip %i %i %i \n\n",j,k,h);				    
		  goto beginning;
		};
	      };
	    };
	  };
	  //
	  stc.str("q");
	  if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
	    printf("Exiting...\n");
	    goto end;
	  };	    
	  //
	  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();
	    stringstream figsave;
	    figsave.str("");
	    figsave << startingdir << "/mip2adcbak_";
	    figsave << l << "_";
	    figsave << m << "_";
	    figsave << n << ".";
	    figsave << input2.str().c_str();
	    c1->SaveAs(figsave.str().c_str());
	    printf("\n");		    
	  };	    
	};	    
	n++;
      };
      m++;
    };
    l++;
  };
 end: printf("\n");
}

void FCaloBAKFIT(TString filename, TString OutDir=""){
  const char *pam_calib = pathtocalibration();
  if ( !strcmp(OutDir.Data(),"") ){
    OutDir = pam_calib;
  };
  //
  // this routine perform the 4224 fit on a backup figure file and save - with extension .bak - the two output file
  // in the correct format, ready to be updated etc. etc. 
  // NB: you will LOSE information about used files
  //
  stringstream file2;
  file2.str("");
  file2 << OutDir.Data() << "/CaloADC2MIP.bak";
  stringstream file4;
  file4.str("");
  file4 << OutDir.Data() << "/CaloADC2MIPf.bak";
  //
  Int_t notgood = 0;
  //
  TFile *hfile = 0;
  CalorimeterCalibration *calo = new CalorimeterCalibration();   
  hfile = new TFile(file2.str().c_str(),"RECREATE","Calorimeter CALIBRATION data");
  if ( !existfile(file2.str().c_str()) ){
    printf(" ERROR: Cannot create file!!! \n");
    return;
  };
  TTree *tree = 0;
  tree = new TTree("CaloADC","Calorimeter calibration data");
  tree->Branch("Event","CalorimeterCalibration",&calo);	    
  //
  const char *file3 = filename;
  TFile hfile3(file3);
  TCanvas *c1;
  c1 = new TCanvas("canvas","canvas",800,800);
  c1->Draw();
  gStyle->SetOptFit(111);
  TH1F *pfmip[2][22][96];
  TF1 *pfitsnr[2][22][96];
  stringstream Fmip;
  for (Int_t i=0; i<2;i++){
    for (Int_t j=0; j<22;j++){
      for (Int_t m=0; m<96;m++){
	Fmip.str("");
	Fmip << "Fmip " << i;
	Fmip << " " << j;
	Fmip << " " << m;
	pfmip[i][j][m] = new TH1F(Fmip.str().c_str(),"",56,0.,55.);
	pfitsnr[i][j][m] = new TF1;
      };
    };	
  };       
  gStyle->SetOptFit(111);
  gPad->SetLogy();
  //
  Int_t l = 0;
  stringstream bmippa;
  while ( l < 2){
    Int_t m = 0;
    while ( m < 22 ){	    
      Int_t n =0 ;
      Double_t SNRPeak = 0.;
      while ( n < 96 ){		    
	c1->Clear();
	c1->cd();		
	bmippa.str("");
	bmippa << "bmip " << l;
	bmippa << " " << m;
	bmippa << " " << n;		
	TH1F *temp  = (TH1F*)hfile3.Get(bmippa.str().c_str());
	bmippa.str("");
	bmippa << "fmip " << l;
	bmippa << " " << m;
	bmippa << " " << n;		
	pfmip[l][m][n] = (TH1F*)temp->Clone(bmippa.str().c_str());
	Int_t doitagain = 0;
	Double_t fr[2];
	Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
	if ( SNRPeak > 16. && SNRPeak < 35. ){
	  fr[0] = (Float_t)SNRPeak - 7.;
	  sv[0]= fp[0]; 
	  sv[1]= fp[1]; 
	  sv[2]= fp[2]; 
	  sv[3]= fp[3];		    
	} else {
	  fr[0] = 19.;
	  sv[0] = 2.8; 
	  sv[1] = 21.0;
	  sv[2] = 1000.0; 
	  sv[3] = 0.1;		    
	};
      fitting: printf("Fitting strip %i %i %i \n",l,m,n);
	fr[1] = 60.;
	pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
	plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
	Double_t chisqr;
	Int_t    ndf;
	pfitsnr[l][m][n] = langaufit(pfmip[l][m][n],fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf,"QR0");
	//
	Double_t SNRPeak = langaupro(fp);
	printf("\n Conversion factor: %f \n\n",SNRPeak);
	//
	//
	//
	if ( (SNRPeak<0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){
	  printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);
	  doitagain++;
	  if ( chisqr > 100. ) {
	    fr[0] = fr[0] + 5.;
	    sv[0] = fp[0]; 
	    sv[1] = fp[1]; 
	    sv[2] = fp[2]; 
	    sv[3] = fp[3];		    
	    goto fitting;
	  };
	  if ( doitagain == 1 ) {
	    fr[0] = 19.;
	    sv[0] = 2.8; 
	    sv[1] = 21.0;
	    sv[2] = 1000.0; 
	    sv[3] = 0.1;
	    goto fitting;
	  };
	  if ( doitagain == 2 ) {
	    fr[0] = 22.;
	    sv[0] = 2.8; 
	    sv[1] = 25.0;
	    sv[2] = 1000.0; 
	    sv[3] = 0.1;
	    goto fitting;
	  };
	  if ( doitagain == 3 ) {
	    fr[0] = 12.;
	    sv[0] = 2.8; 
	    sv[1] = 15.0;
	    sv[2] = 1000.0; 
	    sv[3] = 0.1;
	    goto fitting;
	  };
	};
	//
	calo->mip[l][m][n] = (float)SNRPeak;
	calo->ermip[l][m][n] = (float)fpe[1];
	for (Int_t a = 0; a < 4 ; a++){
	  calo->fp[a][l][m][n] = (float)fp[a];
	  calo->fpe[a][l][m][n] = (float)fpe[a];
	};
	calo->ndf[l][m][n] = (float)ndf;
	calo->chi2[l][m][n] = (float)chisqr;
	//
	if ( ndf!=0 && chisqr/ndf < 2. && fpe[1] < 0.15 && fp[1] > 15. && fp[1] < 40. && SNRPeak > 0. ) {	
	  printf("\n THIS IS A GOOD FIT, SAVED! \n\n");
	  calo->mask[l][m][n] = 1.;			
	} else {
	  notgood++;
	  calo->mask[l][m][n] = 0.;			
	};
	printf("Fitting done, strip %i %i %i \n\n\n",l,m,n);
	c1->cd();		    
	pfmip[l][m][n]->DrawCopy();	    
	pfitsnr[l][m][n]->DrawCopy("lsame");
	c1->Modified();
	c1->Update();
	n++;
      };
      m++;
    };
    l++;
  };
  //
  calo->status = 0;
  if ( !notgood ) calo->status = 1;
  //
  tree->Fill();
  hfile->Write();
  hfile->Close();
  //
  // Save histograms
  //
  TFile *hfile2;
  hfile2 = new TFile(file4.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");	
  if ( !existfile(file4.str().c_str()) ){
    printf(" ERROR: Cannot create file!!! \n");
    return;
  };
  TH1F *sfmip[2][22][96];
  TF1 *sfitsnr[2][22][96];
  stringstream fmipfi;
  for (Int_t l = 0; l < 2; l++){
    for (Int_t m = 0; m < 22; m++){
      for (Int_t n =0 ; n < 96; n++){
	fmipfi.str("");
	fmipfi << "fmip " << l;
	fmipfi << " " << m;
	fmipfi << " " << n;
	sfmip[l][m][n] = new TH1F(fmipfi.str().c_str(),"",56,0.,55.);
	sfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(fmipfi.str().c_str());
	fmipfi.str("");
	fmipfi << "fit " << l;
	fmipfi << " " << m;
	fmipfi << " " << n;
	sfitsnr[l][m][n] = (TF1*)pfitsnr[l][m][n]->Clone(fmipfi.str().c_str());
	sfmip[l][m][n]->Write();
	sfitsnr[l][m][n]->Write();
      };
    };
  };
  hfile2->Close();	
  printf("\n");
  hfile3.Close();
  FCalo4224STATUS("CaloADC2MIP.bak","");
}

void FCaloROOT2TXT(TString rootple, TString txtple){
  FILE *f;
  f = fopen(txtple.Data(),"wb");
  TTree *ctree = 0;
  TFile *chfile;
  CalorimeterCalibration *ccalo = new CalorimeterCalibration();
  chfile = new TFile(rootple.Data(),"READ","Calorimeter CALIBRATION data");
  if ( chfile->IsZombie() ){
    printf(" ERROR: no calorimeter calibration file! \n");
  } else {
    ctree = (TTree*)chfile->Get("CaloADC");  
    ctree->SetBranchAddress("Event", &ccalo);    
    //
    Long64_t cnevents = ctree->GetEntries();
    ctree->GetEntry(cnevents-1);
  };
  //
  Float_t mip = 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++ ){
	mip = 0.;
	if ( (ccalo->fp[1][m][k][l] > 20. && ccalo->fp[1][m][k][l] < 32.) || ccalo->mask[m][k][l] == 1. ) {
	  mip = ccalo->mip[m][k][l];
	} else {
	  mip = 26. ;
	};
	if ( mip == 0 ) mip = 26. ;
	//	if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);
	fwrite(&mip,sizeof(mip),1,f);
      };
    };
  };
  fclose(f);
}

void FCaloREADTXT(TString txtple){
  FILE *f;
  f = fopen(txtple.Data(),"rb");
  Float_t mip = 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++ ){
	mip = 0.;
	fread(&mip,sizeof(mip),1,f);
	//	if ( m == 1 && k == 17 ) printf(" %i mip = %f \n",l,mip);
      };
    };
  };
  fclose(f);
}

void FCaloALIG2DAT(TString txtple){
//   clevel1_.xalig = 121.1;
//   clevel1_.yalig = 119.8;
//   clevel1_.zalig = -263.1;
  char value[256];
  Int_t ixalig;
  Int_t iyalig;
  Int_t izalig;
  Float_t xalig;
  Float_t yalig;
  Float_t zalig;
  //
  FILE *f;
  f = fopen(txtple.Data(),"wb");
  //
  printf("\n Insert four cifres parameters, the last one will become the decimal (1234 will be transformed into 123.4) \n");
  printf("\n Insert the x alignment parameter: \n");
  cin.getline(value,256);
  ixalig = atoi(value);
  xalig = (Float_t)ixalig/10.;
  fwrite(&xalig,sizeof(xalig),1,f);
  printf(" xalig = %f \n",xalig);
  //
  printf("\n Insert the y alignment parameter: \n");
  cin.getline(value,256);
  iyalig = atoi(value);
  yalig = (Float_t)iyalig/10.;
  fwrite(&yalig,sizeof(yalig),1,f);
  printf(" yalig = %f \n",yalig);
  //
  printf("\n Insert the z alignment parameter: \n");
  cin.getline(value,256);
  izalig = atoi(value);
  zalig = (Float_t)izalig/10.;
  fwrite(&zalig,sizeof(zalig),1,f);
  printf(" zalig = %f \n",zalig);
  //
  fclose(f);
}

void FCaloREADALIG(TString txtple){
  FILE *f;
  f = fopen(txtple.Data(),"rb");
  Float_t mip = 0.;
  fread(&mip,sizeof(mip),1,f);
  printf(" xalig = %f \n",mip);
  mip = 0.;
  fread(&mip,sizeof(mip),1,f);
  printf(" yalig = %f \n",mip);
  mip = 0.;
  fread(&mip,sizeof(mip),1,f);
  printf(" zalig = %f \n",mip);
  fclose(f);
}