CalorimeterLevel2/src/CaloProcessing.cpp

/**
 * \file src/CaloProcessing.cpp
 * \author Emiliano Mocchiutti
**/
//
// C/C++ headers
//
#include <sstream>
#include <fstream>
//
// ROOT headers
//
#include <TTree.h>
#include <TBranch.h>
#include <TFile.h>
#include <TObject.h>
//
// YODA headers
//
#include <PamelaRun.h>
#include <physics/calorimeter/CalorimeterEvent.h>
#include <CalibCalPedEvent.h>
//
//
//
#include <GLTables.h>
//
// this package headers
//
#include <delay.h>
#include <CaloProcessing.h>
//
//
// Declaration of the core fortran routines
//
#define calol2cm calol2cm_
extern "C" int calol2cm();
#define calol2tr calol2tr_
extern "C" int calol2tr();
//
using namespace std;
//
//
// Public methods
//

CaloProcessing::~CaloProcessing(){
  delete de;
  delete this;
}

CaloProcessing::CaloProcessing(){
  //
  extern struct FlCaLevel1 clevel1_;
  extern struct FlCaLevel2 clevel2_;
  clevel1 = &clevel1_;
  clevel2 = &clevel2_;
  //
  trkseqno = 0;
  ClearStructs();
  //
  memset(dexy, 0, 2*22*96*sizeof(Float_t));
  memset(dexyc, 0, 2*22*96*sizeof(Float_t));
  memset(mip, 0, 2*22*96*sizeof(Float_t));
  memset(base, 0, 2*22*6*sizeof(Float_t));
  memset(sbase, 0, 2*22*6*sizeof(Float_t));
  calopar1 = true;
  calopar2 = true;
  calopar3 = true;
  ftcalopar1 = 0;
  ttcalopar1 = 0;
  ftcalopar2 = 0;
  ttcalopar2 = 0;
  ftcalopar3 = 0;
  ttcalopar3 = 0;
}

/**
 * Initialize CaloProcessing object
**/
void CaloProcessing::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
  //
  debug = isdeb;
  verbose = isverb;
  //
  l0tr=(TTree*)l0tree;
  de = new pamela::calorimeter::CalorimeterEvent();
  l0calo = (TBranch*)l0tr->GetBranch("Calorimeter");
  l0tr->SetBranchAddress("Calorimeter", &de);
  //
  trkseqno = 0;
  ClearStructs();
  //
  GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
  //
  sgnl = 0;
  UInt_t uptime = 0;
  //
  for (Int_t s = 0; s < 4; s++){
    idcalib[s] = 0;
    fromtime[s] = 0;
    totime[s] = 0;
    calibno[s] = 0;
    ClearCalibVals(s);
    //
    sgnl = glcalo->Query_GL_CALO_CALIB(hs,uptime,s,dbc);
    if ( sgnl < 0 ){
        if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
        return;
    };
    //  
    idcalib[s] = glcalo->ID_ROOT_L0;
    fromtime[s] = glcalo->FROM_TIME;
    if ( glcalo->TO_TIME < hs ){ // calibration is corrupted and we are using the one that preceed the good one
      totime[s] = uptime;
    } else {
      totime[s] = glcalo->TO_TIME;
    };
    calibno[s] = glcalo->EV_ROOT;
    //
    if ( totime[s] == 0 ){
      if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
      ClearCalibVals(s);
      sgnl = 100;      
    };
  };
  //
  // determine path and name and entry of the calibration file
  //
  GL_ROOT *glroot = new GL_ROOT();  
  if ( verbose ) printf("\n");
  for (Int_t s = 0; s < 4; s++){
    if ( verbose ) printf(" ** SECTION %i **\n",s);
    if ( totime[s] > 0 ){
      //
      sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
      if ( sgnl < 0 ){
        if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
        return;
      };
      //        
      stringstream name;
      name.str("");
      name << glroot->PATH.Data() << "/";
      name << glroot->NAME.Data();
      //
      fcalname[s] = (TString)name.str().c_str();
      if ( verbose ) printf(" - runheader at time %u. From time %u to time %u \n   use file %s \n   calibration at entry %i \n\n",hs,fromtime[s],totime[s],fcalname[s].Data(),calibno[s]);
    } else {
      if ( verbose ) printf(" - runheader at time %u. NO CALIBRATION INCLUDE THE RUNHEADER! ",hs);
    };
    sgnl = LoadCalib(s);
    if ( sgnl ) break;
  };    
  //
  delete glcalo;
  delete glroot;
  //
  return;
  //
}

Int_t CaloProcessing::ChkCalib(TSQLServer *dbc, UInt_t atime){
  Int_t sgnl = 0;
  for ( Int_t s = 0; s < 4; s++){
    if ( atime > totime[s] ){
      if ( !dbc->IsConnected() ) throw -116;
      sgnl = Update(dbc,atime,s);
      if ( sgnl < 0 ) return(sgnl); 
    };
  };
  return(sgnl);
}

Int_t CaloProcessing::ChkParam(TSQLServer *dbc, UInt_t runheader){
  stringstream calfile;
  stringstream bmfile;
  stringstream aligfile;
  Int_t error = 0;
  FILE *f = 0;
  ifstream badfile;
  GL_PARAM *glparam = new GL_PARAM();
  //
  if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){
    calopar1 = false;
    //
    // determine where I can find calorimeter ADC to MIP conversion file  
    //
    if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n");
    //
    error = 0;
    error = glparam->Query_GL_PARAM(runheader,101,dbc);
    if ( error < 0 ) return(error);
    //
    calfile.str("");
    calfile << glparam->PATH.Data() << "/";
    calfile << glparam->NAME.Data();
    ftcalopar1 = glparam->FROM_TIME;
    ttcalopar1 = glparam->TO_TIME;
    //
    if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
    f = fopen(calfile.str().c_str(),"rb");
    if ( !f ){
      if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n");
      return(-105);
    };
    //
    for (Int_t m = 0; m < 2 ; m++ ){
      for (Int_t k = 0; k < 22; k++ ){
        for (Int_t l = 0; l < 96; l++ ){
          fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);
        };
      };
    };
    fclose(f);
  };
  //
  if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){
    calopar2 = false;
    //
    // determine where I can find calorimeter alignment file  
    //
    //
    error = 0;
    error = glparam->Query_GL_PARAM(runheader,102,dbc);
    if ( error < 0 ) return(error);
    //
    aligfile.str("");
    aligfile << glparam->PATH.Data() << "/";
    aligfile << glparam->NAME.Data();
    ftcalopar2 = glparam->FROM_TIME;
    ttcalopar2 = glparam->TO_TIME;
    //
    if ( verbose ) printf("\n Using alignment file: \n %s \n",aligfile.str().c_str());
    f = fopen(aligfile.str().c_str(),"rb");
    if ( !f ){
      if ( verbose ) printf(" CALORIMETER - ERROR: no alignement file!\n");
      return(-106);
    };
    //
    fread(&clevel1->xalig,sizeof(clevel1->xalig),1,f);
    if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
    fread(&clevel1->yalig,sizeof(clevel1->yalig),1,f);
    if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
    fread(&clevel1->zalig,sizeof(clevel1->zalig),1,f);
    if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
    fread(&clevel1->emin,sizeof(clevel1->emin),1,f);  
    if ( debug ) printf(" signal threshold = %f \n",clevel1->emin);
    //
    fclose(f);
  };
  //
  // Load offline bad strip mask 
  //
  if ( calopar3 || ( ttcalopar3 != 0 && ttcalopar3 < runheader ) ){
    calopar3 = false;
    //
    // determine where I can find calorimeter alignment file  
    //
    //
    error = 0;
    error = glparam->Query_GL_PARAM(runheader,103,dbc);
    if ( error < 0 ) return(error);
    //
    bmfile.str("");
    bmfile << glparam->PATH.Data() << "/";
    bmfile << glparam->NAME.Data();
    ftcalopar3 = glparam->FROM_TIME;
    ttcalopar3 = glparam->TO_TIME;
    //
    if ( verbose ) printf("\n Using bad strip offline mask file: \n %s \n\n",bmfile.str().c_str());
    badfile.open(bmfile.str().c_str());
    if ( !badfile ){
      if ( verbose ) printf(" CALORIMETER - ERROR: no bad strip offline mask file!\n");
      return(-115);
    };
    //
    Bool_t isdone = false;
    Int_t bad = 0;
    Int_t view = 1;
    Int_t strip = 0;
    Int_t plane = 21;
    while ( !isdone ) {
      badfile >> bad;
      obadmask[view][plane][strip] = bad;
      if ( debug && bad ) printf(" SETTING view %i plane %i strip %i  BAD = %i \n",view,plane,strip,bad);
      strip++;
      if ( strip > 95 ){
        strip = 0;
        plane--;
        if ( plane < 0 ){
          plane = 21;
          view--;
        };
        if ( view < 0 ) isdone = true;
      };
    };
    //
    badfile.close();
  };
  //
  delete glparam;
  //
  return(0);
}


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

Int_t CaloProcessing::Calibrate(Int_t ei){
  //
  // get entry ei
  //
  l0calo->GetEntry(ei);
  //
  // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3.
  //
  Int_t val = 0;
  Int_t del = 1100;
  if ( clevel2->trigty != 2. ){
    for (Int_t sec = 0; sec < 4; sec++){
      val = (Int_t)de->calselftrig[sec][6];
      del = delay(val);
      if ( del < 1100 ){
        clevel2->trigty = 3.;
        break;
      };
    };
  };
  //
  Int_t se = 5;
  Int_t done = 0;
  Int_t pre = -1; 
  Bool_t isCOMP = false;
  Bool_t isFULL = false;
  Bool_t isRAW = false;
  Float_t ener;
  Int_t doneb = 0;
  Int_t donec = 0;
  Int_t ck = 0;
  Int_t ipre = 0;
  Int_t ip[3] = {0};
  Float_t base0, base1, base2;
  base0 = 0.;
  base1 = 0.;
  base2 = 0.;
  Float_t qpre[6] = {0.,0.,0.,0.,0.,0.};
  Float_t ene[96];
  Int_t chdone[4] = {0,0,0,0};
  Int_t pe = 0;
  //
  Float_t ener0 = 0.;
  Float_t cbase0 = 0.;
  Bool_t pproblem = false;
  //
  Float_t tim = 0.;
  Int_t plo = 0;
  Int_t fbi = 0;
  Int_t cle = 0;
  //
  // run over views and planes
  //
  for (Int_t l = 0; l < 2; l++){
    for (Int_t m = 0; m < 22; m++){
      //
      // determine the section number
      //
      se = 5;
      if (l == 0 && m%2 == 0) se = 3;
      if (l == 0 && m%2 != 0) se = 2;
      if (l == 1 && m%2 == 0) se = 1;
      if (l == 1 && m%2 != 0) se = 0;           
      //
      // determine what kind of event we are going to analyze
      //
      isCOMP = false;
      isFULL = false;
      isRAW = false;
      if ( de->stwerr[se] & (1 << 16) ) isCOMP = true; 
      if ( de->stwerr[se] & (1 << 17) ) isFULL = true; 
      if ( de->stwerr[se] & (1 << 3) ) isRAW = true; 
      if ( !chdone[se] ){
        //
        // check for any error in the event
        //
        clevel2->crc[se] = 0;
        if ( de->perror[se] == 132 ){
          clevel2->crc[se] = 1;
          pe++;
        };
        clevel2->perr[se] = 0;
        if ( de->perror[se] != 0 ){
          clevel2->perr[se] = 1;
          pe++;
        };
        clevel2->swerr[se] = 0;
        for (Int_t j = 0; j < 7 ; j++){
          if ( (j != 3) && (de->stwerr[se] & (1 << j)) ){
            clevel2->swerr[se] = 1; 
            pe++;
          };
        };          
        chdone[se] = 1;
      };
      if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){
        pre = -1; 
        //
        for (Int_t nn = 0; nn < 96; nn++){                  
          ene[nn] = 0.;
          dexy[l][m][nn] = de->dexy[l][m][nn] ; 
          dexyc[l][m][nn] = de->dexyc[l][m][nn] ; 
        };
        //
        // run over preamplifiers
        //
        pre = -1;
        cbase0 = 0.;
        for (Int_t i = 0; i < 3; i++){
          for (Int_t j = 0; j < 2; j++){
            pre = j + i*2;
            // 
            // baseline check and calculation
            //
            if ( !isRAW ) {
              base[l][m][pre] = de->base[l][m][pre] ;   
              cbase0 += base[l][m][pre];
            } else {
              //
              // if it is a raw event and we haven't checked 
              // yet, calculate the baseline. 
              //
              FindBaseRaw(l,m,pre);
              cbase0 += base[l][m][pre];
            };
          };
        };
        //
        // run over strips
        //
        pre = -1;
        ener0 = 0.;
        for (Int_t i = 0 ; i < 3 ; i++){
          ip[i] = 0;
          for (Int_t n = i*32 ; n < (i+1)*32 ; n++){                
            if (n%16 == 0) {
              ck = 0;
              done = 0;
              doneb = 0;
              donec = 0;
              pre++;
              qpre[pre] = 0.;
            };
            // 
            // baseline check and calculation
            //
            // no suitable new baseline, use old ones!
            //
            if ( !done ){
              if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){
                ck = 1;
                if (pre%2 == 0) {
                  ip[i] = pre + 1;
                } else { 
                  ip[i] = pre - 1;
                };
                if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0.) ){
                  //
                  ck = 2;
                  if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) {
                    ck = 3;
                  };
                };
                done = 1;
              };
            };                  
            //
            // CALIBRATION ALGORITHM
            // 
            if ( !doneb ){
              switch (ck) {
              case 0: 
                base0 = base[l][m][pre];
                base2 = calbase[l][m][pre];
                break;
              case 1: 
                base0 = base[l][m][ip[i]];
                base2 = calbase[l][m][ip[i]];
                break;
              case 2: 
                base0 = sbase[l][m][pre];
                base2 = calbase[l][m][pre];
                break;
              case 3: 
                base0 = calbase[l][m][pre];
                base2 = calbase[l][m][pre];
                break;
              };
              base1 = calbase[l][m][pre];
              doneb = 1;
            };
            ener = dexyc[l][m][n];
            ener0 += ener;
            clevel1->estrip[n][m][l] = 0.;
            if ( base0>0 && base0 < 30000. ){
              if ( !donec && (base0 - base1 + base2) != 0. ){
                sbase[l][m][pre] = base0 - base1 + base2;
                donec = 1;
              };
              if ( ener > 0. ){
                clevel1->estrip[n][m][l] = (ener - calped[l][m][n] - base0 - base1 + base2)/mip[l][m][n] ;
                //
                // OK, now in estrip we have the energy deposit in MIP of all the strips for this event (at the end of loops of course)
                //
                qpre[pre] += clevel1->estrip[n][m][l];
              };
            };
          };
          if (ck == 1){
            if (ip[i]%2 == 0) {
              ipre = ip[i] + 1;
            } else {
              ipre = ip[i] - 1;
            };
            for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
              clevel1->estrip[j][m][l] += (qpre[ipre] - qpre[ip[i]]) * 0.00478;
            };
          };
          if (ck == 2){
            for (Int_t j = i*32 ; j < (i+1)*32 ; j++){
              ipre = j/16 + 1;
              clevel1->estrip[j][m][l] +=  qpre[ipre] * 0.00478;
            };
          };
        };
        //
        if ( ener0 == 0. && cbase0 == 0. && !pproblem ){
          if ( verbose ) printf(" Calorimeter power problems! event marked as bad \n");
          pproblem = true;
          pe++;
        };
        //
        Int_t j4 = -4;
        Int_t jjj = -3;
        Int_t jj = -2;
        for (Int_t j = 0 ; j < 100 ; j++){
          jj++;
          jjj++;
          j4++;
          if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l];
          if ( jj >= 0 && jj < 96 ){
            if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581;
            if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581;                  
          };
          if ( jjj >= 0 && jjj < 96 ){
            if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] +=  -ene[jjj] * 0.01581;
            if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] +=  -ene[jjj] * 0.01581;
          };
          if ( j4 >= 0 && j4 < 96 ){
            //
            // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================>
            //
            if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || calrms[l][m][j4] > 26 ){
              clevel1->estrip[j4][m][l] = 0.;
            };
            //
            // code and save the energy for each strip in svstrip
            //
            if ( clevel1->estrip[j4][m][l] > clevel1->emin ){
              //
              tim = 100000.;
              plo = m;
              fbi = 0;
              if ( clevel1->estrip[j4][m][l] > 1. ){
                tim = 10000.;
                plo = m;
                fbi = 1;
              };
              if ( clevel1->estrip[j4][m][l] > 10. ){
                tim = 1000.;
                plo = 22 + m;
                fbi = 1;
              };
              if ( clevel1->estrip[j4][m][l] > 100. ){
                tim = 100.;
                plo = 22 + m;
                fbi = 0;
              };
              if ( clevel1->estrip[j4][m][l] > 1000. ){
                tim = 10.;
                plo = 44 + m;
                fbi = 0;
              };
              if ( clevel1->estrip[j4][m][l] > 10000. ){
                tim = 1.;
                plo = 66 + m;
                fbi = 0;
              };
              //
              cle = (Int_t)lroundf(tim*clevel1->estrip[j4][m][l]); 
              //
              if ( l == 0 ){
                //
                // +-PPSSmmmm.mmmm
                //
                svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;
              } else {
                svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);
              };
              //              if ( ei == 770 ) printf(" j %i l %i m %i estrip %f \n",j4,l,m,clevel1->estrip[j4][m][l]);
              //              if ( ei == 770 ) printf(" num lim %i fbi %i tim %f plo %i cle %i \n",numeric_limits<Int_t>::max(),fbi,tim,plo,cle);
              //              if ( ei == 770 ) printf(" svstrip %i \n",svstrip[istrip]);
              //
              istrip++;
            };
          };
        };
        //
      } else {
        for (Int_t nn = 0; nn < 96; nn++){                  
          clevel1->estrip[nn][m][l] = 0.; 
        };
      };
    };
  };
  if ( !pe  ){
    clevel2->good = 1;
  } else {
    clevel2->good = 0;
  };
  return(0);
}

void CaloProcessing::GetTrkVar(){
  calol2tr();
}

void CaloProcessing::FillTrkVar(CaloLevel2 *ca, Int_t nutrk){
  //
  CaloTrkVar *t_ca = new CaloTrkVar();
  //
  t_ca->trkseqno = trkseqno; 
  t_ca->ncore = (Int_t)clevel2->ncore; 
  t_ca->qcore = clevel2->qcore;
  t_ca->noint = (Int_t)clevel2->noint;
  t_ca->ncyl = (Int_t)clevel2->ncyl;
  t_ca->qcyl = clevel2->qcyl;
  t_ca->qtrack = clevel2->qtrack;
  t_ca->qtrackx = clevel2->qtrackx;
  t_ca->qtracky = clevel2->qtracky;
  t_ca->dxtrack = clevel2->dxtrack;
  t_ca->dytrack = clevel2->dytrack;
  t_ca->qlast = clevel2->qlast;
  t_ca->nlast = (Int_t)clevel2->nlast;
  t_ca->qpre = clevel2->qpre;
  t_ca->npre = (Int_t)clevel2->npre;
  t_ca->qpresh = clevel2->qpresh;
  t_ca->npresh = (Int_t)clevel2->npresh;
  t_ca->qtr = clevel2->qtr;
  t_ca->ntr = (Int_t)clevel2->ntr;
  t_ca->planetot = (Int_t)clevel2->planetot;
  t_ca->qmean = clevel2->qmean;
  t_ca->dX0l = clevel2->dX0l;
  t_ca->qlow = clevel2->qlow;
  t_ca->nlow = (Int_t)clevel2->nlow;
  //
  memcpy(t_ca->tibar,clevel2->tibar,sizeof(clevel2->tibar));
  memcpy(t_ca->tbar,clevel2->tbar,sizeof(clevel2->tbar));
  //
  if ( trkseqno == -1 ){
    ca->impx = clevel2->impx;
    ca->impy = clevel2->impy;
    ca->tanx = clevel2->tanx;
    ca->tany = clevel2->tany;
    ca->elen = clevel2->elen;
    ca->selen = clevel2->selen;
    memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
    memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
    memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
    memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
    memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
  };
  // 
  TClonesArray &t = *ca->CaloTrk;
  new(t[nutrk]) CaloTrkVar(*t_ca);
  //
  delete t_ca;
  //
  ClearTrkVar();
}

void CaloProcessing::GetCommonVar(){
  calol2cm();
}

void CaloProcessing::FillCommonVar(CaloLevel2 *ca){
  //
  ca->good = clevel2->good;
  if ( clevel2->trigty == 2. ){
    ca->selftrigger = 1;
  } else {
    ca->selftrigger = 0;
  };
  memcpy(ca->perr,clevel2->perr,sizeof(clevel2->perr));
  memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr));
  memcpy(ca->crc,clevel2->crc,sizeof(clevel2->crc));
  ca->nstrip = (Int_t)clevel2->nstrip;
  ca->qtot = clevel2->qtot;
  ca->impx = clevel2->impx;
  ca->impy = clevel2->impy;
  ca->tanx = clevel2->tanx;
  ca->tany = clevel2->tany;
  ca->nx22 = (Int_t)clevel2->nx22;
  ca->qx22 = clevel2->qx22;
  ca->qmax = clevel2->qmax;
  ca->elen = clevel2->elen;
  ca->selen = clevel2->selen;
  ca->estrip = TArrayI(ca->nstrip,svstrip);
  memcpy(ca->qq,clevel2->qq,sizeof(clevel2->qq));
  memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
  memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
  memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
  memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
  memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
  //
}

void CaloProcessing::ClearStructs(){
  ClearTrkVar();
  ClearCommonVar();
}

void CaloProcessing::RunClose(){
  l0tr->Delete();
  ClearStructs();
  //
  memset(dexy, 0, 2*22*96*sizeof(Float_t));
  memset(dexyc, 0, 2*22*96*sizeof(Float_t));
  memset(base, 0, 2*22*6*sizeof(Float_t));
  memset(sbase, 0, 2*22*6*sizeof(Float_t));
  //
}

//
// Private methods
//

void CaloProcessing::ClearTrkVar(){
  clevel2->ncore = 0; 
  clevel2->qcore = 0.;
  clevel2->noint = 0.;
  clevel2->ncyl = 0.;
  clevel2->qcyl = 0.;
  clevel2->qtrack = 0.;
  clevel2->qtrackx = 0.;
  clevel2->qtracky = 0.;
  clevel2->dxtrack = 0.;
  clevel2->dytrack = 0.;
  clevel2->qlast = 0.;
  clevel2->nlast = 0.;
  clevel2->qpre = 0.;
  clevel2->npre = 0.;
  clevel2->qpresh = 0.;
  clevel2->npresh = 0.;
  clevel2->qlow = 0.;
  clevel2->nlow = 0.;
  clevel2->qtr = 0.;
  clevel2->ntr = 0.;
  clevel2->planetot = 0.;
  clevel2->qmean = 0.;
  clevel2->dX0l = 0.;
  clevel2->elen = 0.;
  clevel2->selen = 0.;
  memset(clevel1->al_p, 0, 5*2*sizeof(Double_t));
  memset(clevel2->tibar, 0, 2*22*sizeof(Int_t));
  memset(clevel2->tbar, 0, 2*22*sizeof(Float_t));
}

void CaloProcessing::ClearCommonVar(){
  istrip = 0;
  clevel2->trigty = -1.;
  clevel2->good = 0;
  clevel2->nstrip = 0.;
  clevel2->qtot = 0.;
  clevel2->impx = 0.;
  clevel2->impy = 0.;
  clevel2->tanx = 0.;
  clevel2->tany = 0.;
  clevel2->qmax = 0.;
  clevel2->nx22 = 0.;
  clevel2->qx22 = 0.;
  memset(clevel2->perr, 0, 4*sizeof(Int_t));
  memset(clevel2->swerr, 0, 4*sizeof(Int_t));
  memset(clevel2->crc, 0, 4*sizeof(Int_t));
  memset(clevel2->qq, 0, 4*sizeof(Int_t));
  memset(clevel2->varcfit, 0, 2*sizeof(Float_t));
  memset(clevel2->npcfit, 0, 2*sizeof(Int_t));
  memset(clevel2->planemax, 0, 2*sizeof(Int_t));
  memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));
  memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));
}

void CaloProcessing::ClearCalibVals(Int_t s){
  //
  for ( Int_t d=0 ; d<11 ;d++  ){
    Int_t pre = -1;
    for ( Int_t j=0; j<96 ;j++){
      if ( j%16 == 0 ) pre++;
      if ( s == 2 ){
        calped[0][2*d+1][j] = 0.;
        cstwerr[3] = 0.;
        cperror[3] = 0.;
        calgood[0][2*d+1][j] = 0.;
        calthr[0][2*d+1][pre] = 0.;
        calrms[0][2*d+1][j] = 0.;
        calbase[0][2*d+1][pre] = 0.;
        calvar[0][2*d+1][pre] = 0.;
      };
      if ( s == 3 ){
        calped[0][2*d][j] = 0.;
        cstwerr[1] = 0.;
        cperror[1] = 0.;
        calgood[0][2*d][j] = 0.;
        calthr[0][2*d][pre] = 0.;
        calrms[0][2*d][j] = 0.;
        calbase[0][2*d][pre] = 0.;
        calvar[0][2*d][pre] = 0.;
      };
      if ( s == 0 ){
        calped[1][2*d][j] = 0.;
        cstwerr[0] = 0.;
        cperror[0] = 0.;
        calgood[1][2*d][j] = 0.;
        calthr[1][2*d][pre] = 0.;
        calrms[1][2*d][j] = 0.;
        calbase[1][2*d][pre] = 0.;
        calvar[1][2*d][pre] = 0.;
      };
      if ( s == 1 ){
        calped[1][2*d+1][j] = 0.;
        cstwerr[2] = 0.;
        cperror[2] = 0.;
        calgood[1][2*d+1][j] = 0.;
        calthr[1][2*d+1][pre] = 0.;
        calrms[1][2*d+1][j] = 0.;
        calbase[1][2*d+1][pre] = 0.;
        calvar[1][2*d+1][pre] = 0.;
      };
    };
  };
  return;
}

Int_t CaloProcessing::Update(TSQLServer *dbc, UInt_t atime, Int_t s){
  //
  Int_t sgnl = 0;
  //
  GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
  //
  sgnl = 0;
  //
  idcalib[s] = 0;
  fromtime[s] = 0;
  totime[s] = 0;
  calibno[s] = 0;
  ClearCalibVals(s);
  //
  UInt_t uptime = 0;
  //
  sgnl = glcalo->Query_GL_CALO_CALIB(atime,uptime,s,dbc);
  if ( sgnl < 0 ){
    if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
    return(sgnl);
  };
  //  
  idcalib[s] = glcalo->ID_ROOT_L0;
  fromtime[s] = glcalo->FROM_TIME;
  if ( glcalo->TO_TIME < atime ){ // calibration is corrupted and we are using the one that preceed the good one
    totime[s] = uptime;
  } else {
    totime[s] = glcalo->TO_TIME;
  };
  //  totime[s] = glcalo->TO_TIME;
  calibno[s] = glcalo->EV_ROOT;
  //
  if ( totime[s] == 0 ){
    if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
    ClearCalibVals(s);
    sgnl = 100;      
  };
  //
  // determine path and name and entry of the calibration file
  //
  GL_ROOT *glroot = new GL_ROOT();  
  if ( verbose ) printf("\n");
  if ( verbose ) printf(" ** SECTION %i **\n",s);
  //
  sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
  if ( sgnl < 0 ){
    if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
    return(sgnl);
  };
  //    
  stringstream name;
  name.str("");
  name << glroot->PATH.Data() << "/";
  name << glroot->NAME.Data();
  //
  fcalname[s] = (TString)name.str().c_str();
  if ( verbose ) printf(" - event at time %u. From time %u to time %u \n   use file %s \n   calibration at entry %i \n\n",atime,fromtime[s],totime[s],fcalname[s].Data(),calibno[s]);
  //
  sgnl = LoadCalib(s);
  //
  if ( sgnl != 0 ) return(sgnl);
  delete glcalo;
  delete glroot;
  //
  return(0);
  //
}

Int_t CaloProcessing::LoadCalib(Int_t s){
  //
  ifstream myfile;
  myfile.open(fcalname[s].Data());
  if ( !myfile ){    
    return(-107);
  };
  myfile.close();
  //
  TFile *File = new TFile(fcalname[s].Data());
  if ( !File ) return(-108);
  TTree *tr = (TTree*)File->Get("CalibCalPed");
  if ( !tr ) return(-109);
  //
  TBranch *calo = tr->GetBranch("CalibCalPed");
  //
  pamela::CalibCalPedEvent *ce = 0;
  tr->SetBranchAddress("CalibCalPed", &ce);
  //
  Long64_t ncalibs = calo->GetEntries();
  //
  if ( !ncalibs ) return(-110); 
  //
  calo->GetEntry(calibno[s]);
  //
  if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
    for ( Int_t d=0 ; d<11 ;d++  ){
      Int_t pre = -1;
      for ( Int_t j=0; j<96 ;j++){
        if ( j%16 == 0 ) pre++;
        if ( s == 2 ){
          calped[0][2*d+1][j] = ce->calped[3][d][j];
          cstwerr[3] = ce->cstwerr[3];
          cperror[3] = ce->cperror[3];
          calgood[0][2*d+1][j] = ce->calgood[3][d][j];
          calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
          calrms[0][2*d+1][j] = ce->calrms[3][d][j];
          calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
          calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
        };
        if ( s == 3 ){
          calped[0][2*d][j] = ce->calped[1][d][j];
          cstwerr[1] = ce->cstwerr[1];
          cperror[1] = ce->cperror[1];
          calgood[0][2*d][j] = ce->calgood[1][d][j];
          calthr[0][2*d][pre] = ce->calthr[1][d][pre];
          calrms[0][2*d][j] = ce->calrms[1][d][j];
          calbase[0][2*d][pre] = ce->calbase[1][d][pre];
          calvar[0][2*d][pre] = ce->calvar[1][d][pre];
        };
        if ( s == 0 ){
          calped[1][2*d][j] = ce->calped[0][d][j];
          cstwerr[0] = ce->cstwerr[0];
          cperror[0] = ce->cperror[0];
          calgood[1][2*d][j] = ce->calgood[0][d][j];
          calthr[1][2*d][pre] = ce->calthr[0][d][pre];
          calrms[1][2*d][j] = ce->calrms[0][d][j];
          calbase[1][2*d][pre] = ce->calbase[0][d][pre];
          calvar[1][2*d][pre] = ce->calvar[0][d][pre];
        };
        if ( s == 1 ){
          calped[1][2*d+1][j] = ce->calped[2][d][j];
          cstwerr[2] = ce->cstwerr[2];
          cperror[2] = ce->cperror[2];
          calgood[1][2*d+1][j] = ce->calgood[2][d][j];
          calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
          calrms[1][2*d+1][j] = ce->calrms[2][d][j];
          calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
          calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
        };
      };
    };
  } else {
    if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
    return(-111);
  };
  File->Close();
  return(0);
}