FUTILITIES/macros/FCaloADC2MIP.cxx

//
//- Emiliano Mocchiutti
//
//   FCaloADC2MIP.c      version 1.01  (2006-03-20)
//
//   The only input needed is 
//
//   Changelog:
//
//   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);
}