CaloEnergy/src/CaloEnergy.cpp

#include <CaloEnergy.h>
#include <PamLevel2.h>

//--------------------------------------
/**
 * Default constructor 
 */
CaloEnergy::CaloEnergy(){
  Clear();
}

CaloEnergy::CaloEnergy(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  fSimu = false;
  this->Set();
  //
}

CaloEnergy::CaloEnergy(PamLevel2 *l2p, Bool_t simulation){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  fSimu = simulation;
  this->Set();
  //
}

void CaloEnergy::UseCaloPreSampler(){
  //
  // use the presampler setting forcefitmode to 1000 means to force the DV routine to find the track inside the calorimeter using the "shower" approach developed for electrons
  //
  if ( !cp ) cp = new CaloPreSampler(L2);
  cp->SplitInto(0,22);
  cp->SetForceFitMode(1000);
//   cp->UseTracker(false);
//   cp->ForceCaloFit();
//   cp->SetDebug(true);
//   cp->Process();
  if ( clong ) clong->SetCaloLevel2Pointer(cp->GetLevel2Pointer());
}


void CaloEnergy::UseLongFit(){ 
  fPl = 0; 
  fLong = true;
  if ( !clong ){
    clong = new CaloLong(L2);
    if ( cp ) clong->SetCaloLevel2Pointer(cp->GetLevel2Pointer());
    clong->SplitInto(0,22);
  };
  //
}

void CaloEnergy::Set(){
  //
  // set default values, NB default conversion factor for energy is just very approximated!
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  sntr = "start";
  //
  AOBT = 0;
  APKT = 0;
  aatime = 0;
  asntr = "start";
  //
  debug = false;
  //
  indep = false;
  //
  fLong = false;
  fPl = 1;
  fRad = -1;
  cp = NULL;
  clong = NULL;
  x0max = -1.;
  //
  this->DefineGeometry();
  fXosel =true;
  fXesel = true;
  fYosel = true;
  fYesel = true;
  fConv_rxe = 44.4;
  fConv_rxo = 44.4;
  fConv_ryo = 44.4;
  fConv_rye = 44.4;
  fXomin = 1000;
  fXemin = 1000;
  fYomin = 1000;
  fYemin = 1000;
}

void CaloEnergy::SetMinimumContainment(TString section, Int_t plane){
  section.ToUpper();
  if ( section.Contains("XO") ) fXomin = plane;
  if ( section.Contains("XE") ) fXemin = plane;
  if ( section.Contains("YO") ) fYomin = plane;
  if ( section.Contains("YE") ) fYemin = plane;
}

void CaloEnergy::SetMinimumContainment(Int_t plane){
  this->SetMinimumContainment("XEXOYEYO",plane);
}

Int_t CaloEnergy::GetMinimumContainment(TString section){
  section.ToUpper();
  if ( section.Contains("XO") ) return(fXomin);
  if ( section.Contains("XE") ) return(fXemin);
  if ( section.Contains("YE") ) return(fYemin);
  if ( section.Contains("YO") ) return(fYomin);
  printf(" ERROR: section not recognized \n");
  return(-1000);
}

void CaloEnergy::SetConversionFactor(TString section, Float_t conv){
  section.ToUpper();
  if ( section.Contains("XO") ) fConv_rxo = conv;
  if ( section.Contains("XE") ) fConv_rxe = conv;
  if ( section.Contains("YO") ) fConv_ryo = conv;
  if ( section.Contains("YE") ) fConv_rye = conv;
}

void CaloEnergy::SetConversionFactor(Float_t conv){
  this->SetConversionFactor("XEXOYEYO",conv);
}

Float_t CaloEnergy::GetConversionFactor(TString section){
  section.ToUpper();
  if ( section.Contains("XO") ) return(fConv_rxo);
  if ( section.Contains("XE") ) return(fConv_rxe);
  if ( section.Contains("YO") ) return(fConv_ryo);
  if ( section.Contains("YE") ) return(fConv_rye);
  printf(" ERROR: section not recognized \n");
  return(-1000.);
}

Int_t CaloEnergy::GetMaxplane(TString section){
  section.ToUpper();
  if ( section.Contains("XO") ) return fMax_planexo;
  if ( section.Contains("XE") ) return fMax_planexe;
  if ( section.Contains("YO") ) return fMax_planeyo;
  if ( section.Contains("YE") ) return fMax_planeye;
  return(-1);
}

Float_t CaloEnergy::GetEnergyAtMaxplane(TString section){
  section.ToUpper();
  if ( section.Contains("XO") ) return xomax_en;
  if ( section.Contains("XE") ) return xemax_en;
  if ( section.Contains("YO") ) return yomax_en;
  if ( section.Contains("YE") ) return yemax_en;
  return(-1);
}

Float_t CaloEnergy::GetMaxEnergy(TString section){
  section.ToUpper();
  if ( fLong ){ 
    this->Process(section);    
    return fXOen_maxplane;
  } else {
    if ( section.Contains("XO") ) return fXOen_maxplane;
    if ( section.Contains("XE") ) return fXEen_maxplane;
    if ( section.Contains("YO") ) return fYOen_maxplane;
    if ( section.Contains("YE") ) return fYEen_maxplane;
  };
  return(-1);
}

Float_t CaloEnergy::GetMaxEnergy(){   
  if ( fLong ){
    if ( debug ) printf(" oh! call process! with asntr %s and sntr %s \n",asntr.Data(),sntr.Data());
    this->Process(asntr);
  };
  return((fXEen_maxplane+fYOen_maxplane+fYEen_maxplane+fXOen_maxplane));
}

void CaloEnergy::DefineGeometry(){
  //
  // Use CaloStrip to determine once the position of border strips for each section
  //
  //
  //  fM = 2. + 0.096; // real position from cbar BUG the 0.096 is already taken into account in the border calculation made by Giovanna
  fM = 2. ; // real position from cbar 
  //  fM1 = 2. - 0.122 - 0.096; // due to calculation of xe1 etc. BUG! this way we count from the silicon border not from the silicon sensitive area
  fM1 = 2. - 0.122 - 0.096 + 0.096; // due to calculation of xe1 etc.
  if ( fM1 < 0. ) fM1 = 0.;
  //
  CaloStrip *cs = new CaloStrip(fSimu);
  //
  // view y plane 0 strip 0
  cs->Set(1,0,0);
  xe1= cs->GetY();
  //  z1= cs->GetZ();
  // view y plane 0 strip 31
  cs->Set(1,0,31);
  xe2= cs->GetY();
  // view y plane 0 strip 32
  cs->Set(1,0,32);
  xe3= cs->GetY();
  // view y plane 0 strip 63
  cs->Set(1,0,63);
  xe4= cs->GetY();
  // view y plane 0 strip 64
  cs->Set(1,0,64);
  xe5= cs->GetY();
  // view y plane 0 strip 95
  cs->Set(1,0,95);
  xe6= cs->GetY();

  //  printf(" xe1 %f xe2 %f xe3 %f xe4 %f xe5 %f xe6 %f \n",xe1,xe2,xe3,xe4,xe5,xe6);

  // view x plane 0 strip 0
  cs->Set(0,0,0);
  yo1= cs->GetX();
  //  z2= cs->GetZ();
  // view x plane 0 strip 31
  cs->Set(0,0,31);
  yo2= cs->GetX();
  // view x plane 0 strip 32
  cs->Set(0,0,32);
  yo3= cs->GetX();
  // view x plane 0 strip 63
  cs->Set(0,0,63);
  yo4= cs->GetX();
  // view x plane 0 strip 64
  cs->Set(0,0,64);
  yo5= cs->GetX();
  // view x plane 0 strip 95
  cs->Set(0,0,95);
  yo6= cs->GetX();

  // view y plane 1 strip 0
  cs->Set(1,1,0);
  xo1= cs->GetY();
  //  z3= cs->GetZ();
  // view y plane 1 strip 31
  cs->Set(1,1,31);
  xo2= cs->GetY();
  // view y plane 1 strip 32
  cs->Set(1,1,32);
  xo3= cs->GetY();
  // view y plane 1 strip 63
  cs->Set(1,1,63);
  xo4= cs->GetY();
  // view y plane 1 strip 64
  cs->Set(1,1,64);
  xo5= cs->GetY();
  // view y plane 1 strip 95
  cs->Set(1,1,95);
  xo6= cs->GetY();
  
  // view x plane 1 strip 0
  cs->Set(0,1,0);
  ye1= cs->GetX();
  //  z4= cs->GetZ();
  // view x plane 1 strip 31
  cs->Set(0,1,31);
  ye2= cs->GetX();
  // view x plane 1 strip 32
  cs->Set(0,1,32);
  ye3= cs->GetX();
  // view x plane 1 strip 63
  cs->Set(0,1,63);
  ye4= cs->GetX();
  // view x plane 1 strip 64
  cs->Set(0,1,64);
  ye5= cs->GetX();
  // view x plane 1 strip 95
  cs->Set(0,1,95);
  ye6= cs->GetX();

  //printf(" ye1 %f ye2 %f ye3 %f ye4 %f ye5 %f ye6 %f \n",ye1,ye2,ye3,ye4,ye5,ye6);
  
  //
  for (Int_t p = 0; p<22; p ++){
    for (Int_t v = 0; v<2; v++ ){   
      cs->Set(v,p,0);
      trk_z[p][v]= cs->GetZ();                  //  Z coord for each plane
    };
  };
  //
  delete cs;
  //
}

void CaloEnergy::Clear(){
  //
  // clear variables
  //
  fPartsel = false;
  fSel = false;
  fXosel = false;
  fXesel = false;
  fYosel = false;
  fYesel = false;
  fCount = 0.;
  fEnergy = 0.;
  fEnergyxe = 0.;
  fEnergyxo = 0.;
  fEnergyye = 0.;
  fEnergyyo = 0.;
  fMax_plane = 0;
  fMax_planexo = 0;
  fMax_planexe = 0;
  fMax_planeyo = 0;
  fMax_planeye = 0; 
  xomax_en= 0.; 
  xemax_en= 0.; 
  yomax_en= 0.; 
  yemax_en= 0.; 
  //
  memset(enstrip,0,2*22*96*(sizeof(Float_t)));  
  en = 0.;
  view = 0;
  plane = 0;
  strip = 0;
  energyxe = 0.; 
  energyyo = 0.; 
  energyxo = 0.; 
  energyye = 0.; 
  fYoout = 0;
  fYeout = 0;
  fXoout = 0;
  fXeout = 0;
  fXEen_maxplane = 0.;  
  fXOen_maxplane = 0.;  
  fYEen_maxplane = 0.;  
  fYOen_maxplane = 0.;  
  memset(en_xep,0,11*sizeof(Float_t));
  memset(en_yep,0,11*sizeof(Float_t));
  memset(en_xop,0,11*sizeof(Float_t));
  memset(en_yop,0,11*sizeof(Float_t));
  //
}

void CaloEnergy::Print(){
  //
  printf("========================================================================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
  printf(" fEnergy :.............. %f \n",fEnergy);  
  printf(" fMax_plane :........... %f \n",fMax_plane); 
  printf(" fMax_planexo :......... %i \n",fMax_planexo); 
  printf(" fMax_planexe :......... %i \n",fMax_planexe); 
  printf(" fMax_planeyo :......... %i \n",fMax_planeyo); 
  printf(" fMax_planeye :......... %i \n",fMax_planeye); 
  printf(" fCount  :.............. %f \n",fCount); 
  printf(" fSel    :.............. %i \n",fSel); 
  printf(" fPartsel:.............. %i \n",fPartsel); 
  printf(" fXesel  :.............. %i \n",fXesel); 
  printf(" fXosel  :.............. %i \n",fXosel); 
  printf(" fYesel  :.............. %i \n",fYesel); 
  printf(" fYosel  :.............. %i \n",fYosel); 
  printf(" fXemin  :.............. %i \n",fXemin); 
  printf(" fXomin  :.............. %i \n",fXomin); 
  printf(" fYemin  :.............. %i \n",fYemin); 
  printf(" fYomin  :.............. %i \n",fYomin); 
  printf(" fXeout  :.............. %i \n",fXeout); 
  printf(" fXoout  :.............. %i \n",fXoout); 
  printf(" fYeout  :.............. %i \n",fYeout); 
  printf(" fYoout  :.............. %i \n",fYoout); 
  printf(" fSimu   :.............. %i \n",fSimu); 
  printf(" fM      :.............. %f \n",fM); 
  printf(" fM1     :.............. %f \n",fM1); 
  printf(" fRad    :.............. %i \n",fRad); 
  printf(" fPl     :.............. %i \n",fPl); 
  printf(" fConv_rxe ............. %f \n",fConv_rxe); 
  printf(" fConv_rxo ............. %f \n",fConv_rxo); 
  printf(" fConv_ryo ............. %f \n",fConv_ryo); 
  printf(" fConv_rye ............. %f \n",fConv_rye); 
  printf(" fLong   :.............. %i \n",fLong); 
  printf(" energyxe:.............. %f \n",energyxe); 
  printf(" energyxo:.............. %f \n",energyxo); 
  printf(" energyye:.............. %f \n",energyye); 
  printf(" energyyo:.............. %f \n",energyyo); 
  printf(" fXEen_maxplane:........ %f \n",fXEen_maxplane);
  printf(" fXOen_maxplane:........ %f \n",fXOen_maxplane);
  printf(" fYEen_maxplane:........ %f \n",fYEen_maxplane);
  printf(" fYOen_maxplane:........ %f \n",fYOen_maxplane);
  printf(" x0max   :.............. %f \n",x0max); 
  printf(" debug   :.............. %i \n",debug); 

  printf("========================================================================\n");
  //
}

void CaloEnergy::Delete(){
  Clear();
  delete this;
}


Bool_t CaloEnergy::IsInsideAcceptance(TString section){
  //
  // check if the event is inside the acceptance of the given section(s)
  //  
  TString ntr = section;
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloEnergy variables not filled \n");
    return false;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != APKT || L2->GetOrbitalInfo()->OBT != AOBT || L2->GetOrbitalInfo()->absTime != aatime || strcmp(ntr.Data(),asntr.Data()) ){
      newentry = true;
      AOBT = L2->GetOrbitalInfo()->OBT;
      APKT = L2->GetOrbitalInfo()->pkt_num;
      aatime = L2->GetOrbitalInfo()->absTime;
      asntr = ntr;
    };
  } else {
    newentry = true;
  };
  //
  // if we have already called this method for this event and no input changed then return fSel and exit
  //
  if ( !newentry ) return fSel;
  //
  // process the event
  //
  if ( debug ) printf(" ########## IsInsideAcceptance ######### \n");
  //
  // clear variables
  //
  this->Clear();
  //
  section.ToUpper();
  //
  // Count the number of section(s) given as input
  //
  Int_t fNumSec = Int_t(section.Contains("XO"))+Int_t(section.Contains("XE"))+Int_t(section.Contains("YO"))+Int_t(section.Contains("YE"));
  if ( !fNumSec ){
    printf(" ERROR: section must be XO or XE or YO or YE while it is %s \n",section.Data());
    return false;
  };
  //
  // If the presampler object exists then use the presampler output instead of the level2 output
  //
  CaloLevel2 *cl2 = NULL;
  if ( cp ){
    cl2 = cp->GetCaloLevel2();
  } else {
    cl2 = L2->GetCaloLevel2();
  };
  //
  // get the energy for every strip of the calorimeter
  //
  for (Int_t ch=0; ch< L2->GetCaloLevel1()->istrip; ch++){
    en = L2->GetCaloLevel1()->DecodeEstrip(ch,view,plane,strip);
    enstrip[view][plane][strip]=en;
  };
  //
  if ( debug && ((fM1+0.122-0.244*(Float_t)fRad) < 0.) ) printf("Error: (fM1+0.122-0.244*(Float_t)fRad) < 0. fM1 %f fRad %i %f \n",fM1,fRad,(fM1+0.122-0.244*(Float_t)fRad));
  //
  // sum energy plane by plane for each sections 
  //
  for (Int_t i=0;i<11;i++){
    for(strip=0; strip<96; strip++) {
      if ( fRad < 0 ){
        //
        // run over all the strips of the plane
        //
        en_xep[i] += enstrip[1][2*i][strip];
        en_yop[i] += enstrip[0][2*i][strip];
        en_xop[i] += enstrip[1][(2*i)+1][strip];
        en_yep[i] += enstrip[0][(2*i)+1][strip];
      } else {
        //
        // use only the strips inside a cylinder of given radius fRad
        //
        if ( cl2->cibar[2*i][1] >= 1 && cl2->cibar[2*i][1] <= 96 && (strip >= (cl2->cibar[2*i][1]-1-fRad)) &&  (strip <= (cl2->cibar[2*i][1]-1+fRad)) ) en_xep[i] += enstrip[1][2*i][strip]; 
        if ( cl2->cibar[2*i][0] >= 1 && cl2->cibar[2*i][0] <= 96 && (strip >= (cl2->cibar[2*i][0]-1-fRad)) &&  (strip <= (cl2->cibar[2*i][0]-1+fRad)) ) en_yop[i] += enstrip[0][2*i][strip]; 
        if ( cl2->cibar[(2*i)+1][1] >= 1 && cl2->cibar[(2*i)+1][1] <= 96 && (strip >= (cl2->cibar[(2*i)+1][1]-1-fRad)) &&  (strip <= (cl2->cibar[(2*i)+1][1]-1+fRad)) ) en_xop[i] += enstrip[1][(2*i)+1][strip]; 
        if ( cl2->cibar[(2*i)+1][0] >= 1 && cl2->cibar[(2*i)+1][0] <= 96 && (strip >= (cl2->cibar[(2*i)+1][0]-1-fRad)) &&  (strip <= (cl2->cibar[(2*i)+1][0]-1+fRad)) ) en_yep[i] += enstrip[0][(2*i)+1][strip]; 
      };
    };
    energyxe += en_xep[i];
    energyyo += en_yop[i];
    energyxo += en_xop[i];
    energyye += en_yep[i];
  };
  //
  // Find the plane of maximum for each section
  //
  //
  Int_t xen = 0;
  Int_t yon = 0;
  Int_t xon = 0;
  Int_t yen = 0;
  //
  if ( section.Contains("XE") ){
    yon++;
    xon++;
    yen++;
    for (Int_t ipl =0; ipl < 11; ipl ++) {
      if(en_xep[ipl] > xemax_en) {
        xemax_en = en_xep[ipl];
        fMax_planexe = ipl;
      };
    };
  };
  //
  if ( section.Contains("YO") ){
    xon++;
    yen++;
    for (Int_t ipl =0; ipl < 11; ipl ++) {
      if(en_yop[ipl] > yomax_en) {
        yomax_en = en_yop[ipl];
        fMax_planeyo = ipl;
      };
    };
  };
  //
  if ( section.Contains("XO") ){
    yen++;
    for (Int_t ipl =0; ipl < 11; ipl ++) {
      if(en_xop[ipl] > xomax_en) {
        xomax_en = en_xop[ipl];
        fMax_planexo = ipl;
      };
    };
  };
  //
  if ( section.Contains("YE") ){
    for (Int_t ipl =0; ipl < 11; ipl ++) {
      if(en_yep[ipl] > yemax_en) {
        yemax_en = en_yep[ipl];
        fMax_planeye = ipl;
      };
    };
  };
  //
  // the maximum is given externally as X0, convert it to plane and section
  //
  if ( x0max > 0. ){
    if ( debug ) printf(" CALCULATE MAX PLANE GIVEN X0 ASSUMING PERPENDICULAR TRACK \n");
    Int_t wpl = (Int_t)roundf(x0max/0.76);
    Bool_t isY = false;
    if ( ((x0max/0.76)-(Float_t)wpl) > 0. ) isY = true;
    xomax_en = 0.;
    yemax_en = 0.;
    xemax_en = 0.;
    yomax_en = 0.;
    //
    if ( !(wpl%2) ){
      // 0, 2, 4, ...
      if ( isY ){
        if ( section.Contains("YO") ) yomax_en = 1000.;
        if ( section.Contains("XE") ) xemax_en = 500.;
        fMax_planeyo=wpl/2;
        fMax_planexe=wpl/2;
        if ( section.Contains("XO") ) xomax_en = 10.;
        if ( section.Contains("YE") ) yemax_en = 5.;
      } else {
        if ( section.Contains("YO") ) yomax_en = 500.;
        if ( section.Contains("XE") ) xemax_en = 1000.;
        fMax_planeyo=wpl/2;
        fMax_planexe=wpl/2;
        if ( section.Contains("XO") ) xomax_en = 5.;
        if ( section.Contains("YE") ) yemax_en = 10.;
      };
    } else {
      // 1, 3, 5, ...
      if ( isY ){
        if ( section.Contains("YE") ) yemax_en = 1000.;
        if ( section.Contains("XO") ) xomax_en = 500.;
        fMax_planeye=(wpl-1)/2;
        fMax_planexo=(wpl-1)/2;
        if ( section.Contains("XE") ) xemax_en = 10.;
        if ( section.Contains("YO") ) yomax_en = 5.;
      } else {
        if ( section.Contains("YE") ) yemax_en = 500.;
        if ( section.Contains("XO") ) xomax_en = 1000.;
        fMax_planeye=(wpl-1)/2;
        fMax_planexo=(wpl-1)/2;
        if ( section.Contains("XE") ) xemax_en = 5.;
        if ( section.Contains("YO") ) yomax_en = 10.;
      };
    };
    if ( debug ) printf(" x0max %f x0max/0.76 %f wpl %i isY %i yomax_en %f xemax_en %f yemax_en %f xomax_en %f fMaxplane %i %i %i %i\n",x0max,(x0max/0.76),wpl,isY,yomax_en,xemax_en,yemax_en,xomax_en,fMax_planeyo,fMax_planexe,fMax_planeye,fMax_planexo);
  };
  //
  Int_t nPl = fPl;
  //
  // Set the maximum in case of coherent mode was selected
  //
  if ( !indep ){
    nPl = 0;
    if ( debug ) printf(" A: Check maximum, coherent mode: xoen %f yoen %f xeen %f yeen %f xomax %i yomax %i xemax %i yemax %i\n",xomax_en,yomax_en,xemax_en,yemax_en,fMax_planexo,fMax_planeyo,fMax_planexe,fMax_planeye);
    Int_t nummod = 0;
    Int_t numexpl = 0;
    if ( xomax_en > xemax_en && xomax_en > yemax_en && xomax_en > yomax_en ){
      //
      // Section XO contains the maximum energy release per plane of the whole calorimeter
      //
      if ( debug ) printf(" XO is MAX %i %i %i %i\n",xen,yon,xon,yen);
      //
      // fMax_plane is the plane of maximum + number of additional dE/dx measurement counting planes from 0 to 43
      //
      fMax_plane = (fNumSec * fMax_planexo) +(Float_t)xon + fPl;
      //
      // nummod is the integer part of the number of modules in which the maximum is contained
      //
      nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
      //
      // numexpl is the number of additional planes (0,1,2) inside the module
      //
      numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
      //
    };
    if ( xemax_en > xomax_en && xemax_en > yemax_en && xemax_en > yomax_en ){ 
      if ( debug ) printf(" XE is MAX %i %i %i %i\n",xen,yon,xon,yen);
      fMax_plane = (fNumSec * fMax_planexe) +(Float_t)xen + fPl;
      nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
      numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
      //
   };

    if ( yemax_en > xomax_en && yemax_en > xemax_en && yemax_en > yomax_en ){
      if ( debug ) printf(" YE is MAX %i %i %i %i\n",xen,yon,xon,yen);
      fMax_plane = (fNumSec * fMax_planeye) +(Float_t)yen + fPl;
      nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
      numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
      //
    };
    if ( yomax_en > xemax_en && yomax_en > yemax_en && yomax_en > xomax_en ){
      if ( debug ) printf(" YO is MAX %i %i %i %i\n",xen,yon,xon,yen);
      fMax_plane = (fNumSec * fMax_planeyo) +(Float_t)yon + fPl;
      nummod = (Int_t)(((Float_t)fMax_plane)/(Float_t)fNumSec);
      numexpl = (Int_t)((Float_t)fMax_plane-(Float_t)fNumSec*(Float_t)nummod);
      //
    };
    //
    // find the plane up to which is necessary to integrate the energy for each section
    //
    Int_t a = 0;
    Int_t b = 0;
    Int_t c = 0;
    if ( numexpl > xen ) a = 1;
    if ( numexpl > yon ) b = 1;
    if ( numexpl > xon ) c = 1;
    fMax_planexe = nummod;
    fMax_planeyo = nummod - 1 + a;
    fMax_planexo = nummod - 1 + b;
    fMax_planeye = nummod - 1 + c;       
    if ( debug ) printf(" fMax_plane %f nummod %i numexpl %i a %i b %i c %i \n",fMax_plane,nummod,numexpl,a,b,c);
    if ( debug ) printf(" DONE: Check maximum, coherent mode: xoen %f yoen %f xeen %f yeen %f xomax %i yomax %i xemax %i yemax %i\n",xomax_en,yomax_en,xemax_en,yemax_en,fMax_planexo,fMax_planeyo,fMax_planexe,fMax_planeye);
  };
  //
  // for each plane of the calorimeter find the position of the track in the direction along the strip (where we do not have a measurement from the selected plane) by looking at the plane above/below of the other view and extrapolating the trajectory to the given plane
  //
  //
  Float_t tgx_cl2;
  Float_t tgy_cl2;  
  tgx_cl2 = cl2->tanx[0];
  tgy_cl2 = cl2->tany[0];
  //
  for (Int_t p=0; p<22; p++){    
    track_coordy[p][1] = cl2->cbar[p][1]; 
    track_coordx[p][1] = cl2->cbar[p][0] - fabs(trk_z[p][1]-trk_z[p][0])*tgx_cl2;
    //    track_coordx[p][1] = cl2->cbar[p][0] + fabs(trk_z[p][1]-trk_z[p][0])*tgx_cl2;
    track_coordx[p][0] = cl2->cbar[p][0];
    track_coordy[p][0] = cl2->cbar[p][1] - fabs(trk_z[p][1]-trk_z[p][0])*tgy_cl2;
    //    track_coordy[p][0] = cl2->cbar[p][1] + fabs(trk_z[p][1]-trk_z[p][0])*tgy_cl2;
    if ( debug ) printf(" p %i track_coordy[p][1] %f track_coordx[p][1] %f track_coordx[p][0] %f  track_coordy[p][0] %f \n",p,track_coordy[p][1],track_coordx[p][1],track_coordx[p][0],track_coordy[p][0]);
  };
  //
  if ( debug ) printf(" acceptance fNumSec %i tgx %f tgy %f\n",fNumSec,tgx_cl2,tgy_cl2);
  //
  if ( section.Contains("XO") ){
    //
    // check event is inside XO acceptance
    //
    for (Int_t i=0; i<11; i++) {
      if
        (
         (((track_coordx[(2*i)+1][1]>=(-12.054+fM))&&
           (track_coordx[(2*i)+1][1]<=(-4.246-fM)))&&
          (((cl2->cbar[(2*i)+1][1]>=xo1 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo2 - fM1 ))||
           ((cl2->cbar[(2*i)+1][1]>=xo3 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo4 - fM1 ))||
           ((cl2->cbar[(2*i)+1][1]>=xo5 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo6 - fM1 ))))||
         
         (((track_coordx[(2*i)+1][1]>=(-4.004+fM))&&
           (track_coordx[(2*i)+1][1]<=(3.804-fM)))&&
          (((cl2->cbar[(2*i)+1][1]>=xo1 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo2 - fM1 ))||
           ((cl2->cbar[(2*i)+1][1]>=xo3 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo4 - fM1))||
           ((cl2->cbar[(2*i)+1][1]>=xo5 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo6 - fM1 ))))||
         
         (((track_coordx[(2*i)+1][1]>=(4.046+fM))&&
           (track_coordx[(2*i)+1][1]<=(11.854-fM)))&&
          (((cl2->cbar[(2*i)+1][1]>=xo1 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo2 - fM1))||
           ((cl2->cbar[(2*i)+1][1]>=xo3 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo4 - fM1 ))||
           ((cl2->cbar[(2*i)+1][1]>=xo5 + fM1)&&
            (cl2->cbar[(2*i)+1][1]<=xo6 - fM1 ))))
         ){
        fXosel = true;
        fXoout = i;
      } else {
        fXosel = false;
        break;
      };
    };
    //
    // if it goes out of the acceptance BUT the plane up to which we are integrating the energy is contained then the event is "partially" contained
    //
    if ( !fXosel && fXoout >= fXomin && fXoout >= (Int_t)(fMax_planexo+nPl) ){
      if ( debug ) printf(" XO: this event is only partially contained: fXoout %i fXomin %i fMax_planexo + nPl %i \n",fXoout,fXomin,(Int_t)(fMax_planexo+nPl));
      fPartsel = true;
      fXosel = true;
    };
    //
    // event is contained (or partially contained) hence we can integrate energy up to the maximum and calculate the energy as measured by this section
    //
    if ( fXosel ){
      for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planexo+nPl)) ;iplm++){        
        fXOen_maxplane += en_xop[iplm]; 
        if ( debug ) printf(" XO iplm %i fXOen_maxplane %f  en_xop[iplm] %f\n",iplm,fXOen_maxplane,en_xop[iplm]);
      };
      fEnergyxo = fXOen_maxplane/fConv_rxo;
    };
  };
  //
  if ( section.Contains("XE") ){
    for (Int_t i=0; i<11; i++) {
      if
        (
         (((track_coordx[2*i][1]>=(-11.854+fM))&& 
           (track_coordx[2*i][1]<=(-4.046-fM)))&&
          (((cl2->cbar[2*i][1]>=xe1 + fM1)&&
            (cl2->cbar[2*i][1]<=xe2 - fM1 ))||
           ((cl2->cbar[2*i][1]>=xe3 + fM1)&&
            (cl2->cbar[2*i][1]<=xe4 - fM1 ))||
           ((cl2->cbar[2*i][1]>=xe5 + fM1)&&
            (cl2->cbar[2*i][1]<=xe6 - fM1 ))))||
         
         (((track_coordx[2*i][1]>=(-3.804+fM))&&
           (track_coordx[2*i][1]<=(4.004-fM)))&&
          (((cl2->cbar[2*i][1]>=xe1 + fM1)&&
            (cl2->cbar[2*i][1]<=xe2 - fM1 ))||
           ((cl2->cbar[2*i][1]>=xe3 + fM1)&&
            (cl2->cbar[2*i][1]<=xe4 - fM1))||
           ((cl2->cbar[2*i][1]>=xe5 + fM1)&&
            (cl2->cbar[2*i][1]<=xe6 - fM1 ))))||
         
         (((track_coordx[2*i][1]>=(4.246+fM))&&
           (track_coordx[2*i][1]<=(12.054-fM)))&&
          (((cl2->cbar[2*i][1]>=xe1 + fM1)&&
            (cl2->cbar[2*i][1]<=xe2 - fM1))||
           ((cl2->cbar[2*i][1]>=xe3 + fM1)&&
            (cl2->cbar[2*i][1]<=xe4 - fM1 ))||
           ((cl2->cbar[2*i][1]>=xe5 + fM1)&&
            (cl2->cbar[2*i][1]<=xe6 - fM1 ))))
         ){
        fXesel = true;
        fXeout = i;
      } else {
        fXesel = false;
        break;
      };
    };
    if ( !fXesel && fXeout >= fXemin && fXeout >= (Int_t)(fMax_planexe+nPl) ){
      if ( debug ) printf(" XE: this event is only partially contained: fXeout %i fXemin %i fMax_planexe + nPl %i \n",fXeout,fXemin,(Int_t)(fMax_planexe+nPl));
      fPartsel = true;
      fXesel = true;
    };
    if ( fXesel ){
      for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planexe+nPl)) ;iplm++){
        fXEen_maxplane += en_xep[iplm]; 
        if ( debug ) printf(" XE iplm %i fXOen_maxplane %f  en_xop[iplm] %f\n",iplm,fXEen_maxplane,en_xep[iplm]);
      };
      fEnergyxe = fXEen_maxplane/fConv_rxe;
    };
  };  
  //
  if ( section.Contains("YE") ){
    for (Int_t i=0; i<11; i++) {
      if
        (
         (((track_coordy[(2*i)+1][0]>=(-12.154+fM))&& 
           (track_coordy[(2*i)+1][0]<=(-4.346-fM)))&&
          (((cl2->cbar[(2*i)+1][0]>=ye1 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye2 - fM1 ))||
           ((cl2->cbar[(2*i)+1][0]>=ye3 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye4 - fM1 ))||
           ((cl2->cbar[(2*i)+1][0]>=ye5 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye6 - fM1 ))))||
         
         (((track_coordy[(2*i)+1][0]>=(-4.104+fM))&&
           (track_coordy[(2*i)+1][0]<=(3.704-fM)))&&
          (((cl2->cbar[(2*i)+1][0]>=ye1 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye2 - fM1 ))||
           ((cl2->cbar[(2*i)+1][0]>=ye3 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye4 - fM1))||
           ((cl2->cbar[(2*i)+1][0]>=ye5 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye6 - fM1 ))))||
         
         (((track_coordy[(2*i)+1][0]>=(3.946+fM))&&
           (track_coordy[(2*i)+1][0]<=(11.754-fM)))&&
          (((cl2->cbar[(2*i)+1][0]>=ye1 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye2 - fM1))||
           ((cl2->cbar[(2*i)+1][0]>=ye3 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye4 - fM1 ))||
           ((cl2->cbar[(2*i)+1][0]>=ye5 + fM1)&&
            (cl2->cbar[(2*i)+1][0]<=ye6 - fM1 ))))
         ){
        fYesel = true;
        fYeout = i;
      } else {
        fYesel = false;
        break;
      };
    };
    if ( !fYesel && fYeout >= fYemin && fYeout >= (Int_t)(fMax_planeye+nPl) ){
      if ( debug ) printf(" YE: this event is only partially contained: fYeout %i fYemin %i fMax_planeye + nPl %i \n",fYeout,fYemin,(Int_t)(fMax_planeye+nPl));
      fPartsel = true;
      fYesel = true;
    };
    if ( fYesel ){
      for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planeye+nPl)) ;iplm++) fYEen_maxplane += en_yep[iplm]; 
      fEnergyye = fYEen_maxplane/fConv_rye;
    };
  };  
  //
  if ( section.Contains("YO") ){
    for (Int_t i=0; i<11; i++) {
      if ( debug ) printf(" i %i track_coordy[2*i][0] %f cl2->cbar[2*i][0] %f \n",i,track_coordy[2*i][0],cl2->cbar[2*i][0]);
      if ( debug ) printf(" i %i fm %f fm1 %f yo1 %g yo2 %f yo3 %f yo4 %f yo5 %f yo6 %f \n",i,fM,fM1,yo1,yo2,yo3,yo4,yo5,yo6);
      if
        (
         (((track_coordy[2*i][0]>=(-11.954+fM))&& 
           (track_coordy[2*i][0]<=(-4.146-fM)))&&
          (((cl2->cbar[2*i][0]>=yo1 + fM1)&&
            (cl2->cbar[2*i][0]<=yo2 - fM1 ))||
           ((cl2->cbar[2*i][0]>=yo3 + fM1)&&
            (cl2->cbar[2*i][0]<=yo4 - fM1 ))||
           ((cl2->cbar[2*i][0]>=yo5 + fM1)&&
            (cl2->cbar[2*i][0]<=yo6 - fM1 ))))||
         
         (((track_coordy[2*i][0]>=(-3.904+fM))&&
           (track_coordy[2*i][0]<=(+3.904-fM)))&&
          (((cl2->cbar[2*i][0]>=yo1 + fM1)&&
            (cl2->cbar[2*i][0]<=yo2 - fM1 ))||
           ((cl2->cbar[2*i][0]>=yo3 + fM1)&&
            (cl2->cbar[2*i][0]<=yo4 - fM1))||
           ((cl2->cbar[2*i][0]>=yo5 + fM1)&&
            (cl2->cbar[2*i][0]<=yo6 - fM1 ))))||
         
         (((track_coordy[2*i][0]>=(4.146+fM))&&
           (track_coordy[2*i][0]<=(11.954-fM)))&&
          (((cl2->cbar[2*i][0]>=yo1 + fM1)&&
            (cl2->cbar[2*i][0]<=yo2 - fM1))||
           ((cl2->cbar[2*i][0]>=yo3 + fM1)&&
            (cl2->cbar[2*i][0]<=yo4 - fM1 ))||
           ((cl2->cbar[2*i][0]>=yo5 + fM1)&&
            (cl2->cbar[2*i][0]<=yo6 - fM1 ))))
         ){
        fYosel = true;
        fYoout = i;
      } else {
        fYosel = false;
        break;
      };
    };
    if ( !fYosel && fYoout >= fYomin && fYoout >= (Int_t)(fMax_planeyo+nPl) ){
      if ( debug ) printf(" YO: this event is only partially contained: fYoout %i fYomin %i fMax_planeyo + nPl %i \n",fYoout,fYomin,(Int_t)(fMax_planeyo+nPl));
      fPartsel = true;
      fYosel = true;
    };
    if ( fYosel ){
      for (Int_t iplm=0;iplm<=TMath::Min(10,(Int_t)(fMax_planeyo+nPl)) ;iplm++) fYOen_maxplane += en_yop[iplm]; 
      fEnergyyo = fYOen_maxplane/fConv_ryo;
    };
  };  
  //
  // Count the number of sections in which the event is contained
  //
  fCount = (Float_t)((Int_t)fXesel+(Int_t)fXosel+(Int_t)fYesel+(Int_t)fYosel);
  //
  if ( indep ){
    //
    // independent mode, average the energy measurement and max plane of the contained sections
    //
    fSel = ( fXesel || fYesel || fXosel || fYosel );
    fMax_plane = (Float_t)(fMax_planeyo+fMax_planeye+fMax_planexo+fMax_planexe)/fCount;
    fEnergy = (fEnergyxe+fEnergyyo+fEnergyye+fEnergyxo)/fCount;
    //
  } else {
    //
    // coherent mode, sum the energy [MIP] of the given sections and convert using fConv_rxo. **** NB: it is assumed that the conversion factor is unique and the method SetConvertionFactor(Float_t) has been used**** The event is selected only if it is contained in all the given sections
    //
    if ( fCount != fNumSec ){
      fSel = false;
    } else {
      fSel = true;
    };
    fEnergy = (fXEen_maxplane+fYOen_maxplane+fYEen_maxplane+fXOen_maxplane)/fConv_rxo;
  };
  //
  if ( debug ) printf("sel %i indep %i fMax_plane %f conv_r %f en_maxplane %f encalo %f  \n",fSel,indep,fMax_plane,fConv_rxo,fXOen_maxplane,fEnergy);
  if ( debug ) printf(" IsInside XE %i XO %i YE %i YO %i => SEL %i \n",fXesel,fXosel,fYesel,fYosel,fSel);
  //
  // finish exit
  //
  return fSel;
  // 
}

void CaloEnergy::Process(){
  TString xo = "XO";
  this->Process(xo);
}


void CaloEnergy::Process(TString section){
  //  
  // process the event
  //
  TString ntr = section;
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloEnergy variables not filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || strcmp(ntr.Data(),sntr.Data()) ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
      sntr = ntr;
    };
  } else {
    newentry = true;
  };
  //
  // if we have already called this method for this event and no input changed then return fSel and exit
  //
  if ( !newentry ) return;
  //
  // process the event
  //
  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u section %s\n",OBT,PKT,atime,section.Data());
  //
  // check if the cylinder of integration can go out of the sensor given the frame which has been set (if we use all the calorimeter fRad is < 0 and the printout is suppressed)
  //
  if ( (fM1+0.122-0.244*(Float_t)fRad) < 0. ) printf("Error: (fM1+0.122-0.244*(Float_t)fRad) < 0. fM1 %f fRad %i %f \n",fM1,fRad,(fM1+0.122-0.244*(Float_t)fRad));
  //
  if ( fLong ){
    if ( debug ) printf(" ==================================================================> LONGITUDINAL FIT! \n"); 
    //
    // use long fit to measure energy 
    //
    if ( this->IsInsideAcceptance(section) ){
      //
      if ( debug ) printf(" ==================================================================> LONG INSIDE! \n"); 
      //
      Float_t myene[2][22];
      memset(myene,0,(sizeof(Float_t))*2*22);
      for (Int_t j=0; j<11; j++){
        if ( section.Contains("XE") ) myene[1][2*j] = en_xep[j];
        if ( section.Contains("YO") ) myene[0][2*j] = en_yop[j];
        if ( section.Contains("XO") ) myene[1][(2*j)+1] = en_xop[j];
        if ( section.Contains("YE") ) myene[0][(2*j)+1] = en_yep[j];    
      };
      clong->UnMaskSections();
      if ( !(section.Contains("YE")) ) clong->MaskSection("YE");
      if ( !(section.Contains("YO")) ) clong->MaskSection("YO");
      if ( !(section.Contains("XO")) ) clong->MaskSection("XO");
      if ( !(section.Contains("XE")) ) clong->MaskSection("XE");
      clong->ForceNextFit();
      clong->SetEnergies(myene);
      if ( debug ){
        clong->Fit(true);
      } else {
        clong->Fit();
      };      
      if ( clong->GetLowerLimit() != 0. || clong->GetUpperLimit() != 0. ){
        fXOen_maxplane = clong->Get_defE0();
      } else {
        fXOen_maxplane = clong->Get_E0();
      };
      fMax_plane = clong->Get_tmax();
      fYOen_maxplane  = 0.;
      fYEen_maxplane = 0.; 
      fXEen_maxplane = 0.;
      fEnergy=fXOen_maxplane/fConv_rxo;
      if ( fEnergy != fEnergy || clong->Get_fitresult() != 0 ) fEnergy = -1.;
      //      if ( fEnergy != fEnergy ) fEnergy = 1.;
      //
    } else {
      //
      // if the event is not in the acceptance, return a negative energy.
      //
      if ( debug ) printf(" Outside acceptance \n");
      fEnergy *= -1.;
      //
    };
    //
  } else {
    //
    // use the energy measurement
    //
    if ( this->IsInsideAcceptance(section) ){
      //
      // the event is good
      //
      if ( debug ) printf(" XE %i XO %i YE %i YO %i \n",fXesel,fXosel,fYesel,fYosel);
      //
    } else {
      //
      // if the event is not in the acceptance, return a negative energy.
      //
      if ( debug ) printf(" Outside acceptance \n");
      fEnergy *= -1.;
      //
    };
  };
  //
}