EndcapPiZeroDiscriminatorAlgo.cc

Go to the documentation of this file.

#include "RecoEcal/EgammaClusterAlgos/interface/EndcapPiZeroDiscriminatorAlgo.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include <fstream>
#include <iostream>

using namespace std;

namespace {
  constexpr int Nfiles_EE = 5;
}

EndcapPiZeroDiscriminatorAlgo::EndcapPiZeroDiscriminatorAlgo(double stripEnergyCut, int nStripCut, const string& path) :
   preshStripEnergyCut_(stripEnergyCut),  preshSeededNstr_(nStripCut)
{   
   
     // Read all Weight files
     constexpr std::array<char const*,Nfiles_EE> file_pt{ {"20","30","40","50","60"} };
     constexpr std::array<char const*,5> file_barrel_pt { {"20","30","40","50","60"} };   

     for(auto ptName: file_pt) {
       
       string nn_paterns_file ="endcapPiZeroDiscriminatorWeights_et";
       nn_paterns_file +=ptName;
       nn_paterns_file +=".wts";
       edm::FileInPath WFile(path+nn_paterns_file);
       readWeightFile(WFile.fullPath().c_str(), EE_Layers, EE_Indim, EE_Hidden,EE_Outdim); // read the weights' file
     }

     for(auto ptName: file_barrel_pt) {
       string nn_paterns_file = "barrelPiZeroDiscriminatorWeights_et";
       nn_paterns_file +=ptName;
       nn_paterns_file +=".wts";
       edm::FileInPath WFile(path+nn_paterns_file);
       readWeightFile(WFile.fullPath().c_str(), EB_Layers, EB_Indim, EB_Hidden, EB_Outdim); // read the weights' file
     }
}


vector<float> EndcapPiZeroDiscriminatorAlgo::findPreshVector(ESDetId strip,  RecHitsMap *rechits_map,
                                                     CaloSubdetectorTopology *topology_p)
{
  vector<float> vout_stripE;

  // skip if rechits_map contains no hits
  if ( rechits_map->empty() ) {
          edm::LogWarning("EndcapPiZeroDiscriminatorAlgo") << "RecHitsMap has size 0.";
          return vout_stripE;
  }

  vout_stripE.clear();

  vector<ESDetId> road_2d;
  road_2d.clear();

  int plane = strip.plane();

  LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: Preshower Seeded Algorithm - looking for clusters" << "n"<< "findPreshVectors: Preshower is intersected at strip " << strip.strip() << ", at plane " << plane ;


  if ( strip == ESDetId(0) ) { //works in case of no intersected strip found
    for(int i=0;i<11;i++) {
       vout_stripE.push_back(-100.);
    }
  }

  // Add to the road the central strip
  road_2d.push_back(strip);

  //Make a navigator, and set it to the strip cell.
  EcalPreshowerNavigator navigator(strip, topology_p);
  navigator.setHome(strip);
 //search for neighbours in the central road
  findPi0Road(strip, navigator, plane, road_2d);
  
  LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo:findPreshVectors: Total number of strips in the central road: " << road_2d.size() ;
  
  // Find the energy of each strip
  RecHitsMap::iterator final_strip =  rechits_map->end();
  // very dangerous, added a protection on the rechits_map->size()
  // at the beginning of the method
  final_strip--;
  ESDetId last_stripID = final_strip->first;

  vector<ESDetId>::iterator itID;
  for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: ID = " << *itID ;
  
    float E = 0.;
    RecHitsMap::iterator strip_it = rechits_map->find(*itID);
    if(goodPi0Strip(strip_it,last_stripID)) { // continue if strip not found in rechit_map
      E = strip_it->second.energy();
    }
    vout_stripE.push_back(E);
    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: E = " << E ;
    
  }
 
  // ***ML beg***
  // vector of size=11, content of vout_stripE is copied into vout_ElevenStrips_Energy
  // to avoid problem in case number of strips is less than 11
  vector<float> vout_ElevenStrips_Energy;
  for(int i=0;i<11;i++)
  {
    vout_ElevenStrips_Energy.push_back(0.);
  }
  
  for(unsigned int i=0;i<vout_stripE.size();i++)
  {
    vout_ElevenStrips_Energy[i] = vout_stripE.at(i);
  }

  //return vout_stripE;
  return vout_ElevenStrips_Energy;
  // ***ML end***

}

// returns true if the candidate strip fulfills the requirements to be added to the cluster:
bool EndcapPiZeroDiscriminatorAlgo::goodPi0Strip(RecHitsMap::iterator candidate_it, ESDetId lastID)
{
  RecHitsMap::iterator candidate_tmp = candidate_it;
  candidate_tmp--;

  // crystal should not be included...
  if ( candidate_tmp->first == lastID )                              // ...if it corresponds to a hit
    {
      LogTrace("EcalClusters") <<"EndcapPiZeroDiscriminatorAlgo: goodPi0Strip No such a strip in rechits_map " ;
      return false;
    }
  else if (candidate_it->second.energy() <= preshStripEnergyCut_ )   // ...if it has a negative or zero energy
    {
      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: goodPi0Strip Strip energy " << candidate_it->second.energy() <<" is below threshold " ;
      return false;
    }

  return true;
}

// find strips in the road of size +/- preshSeededNstr_ from the central strip
void EndcapPiZeroDiscriminatorAlgo::findPi0Road(ESDetId strip, EcalPreshowerNavigator& theESNav,
                                                                int plane, vector<ESDetId>& vout) {
  if ( strip == ESDetId(0) ) return;
   ESDetId next;
   theESNav.setHome(strip);
   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: starts from strip " << strip ;
   
   if (plane == 1) {
     // east road
     int n_east= 0;
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the East " ;
   
     while ( ((next=theESNav.east()) != ESDetId(0) && next != strip) ) {
        LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: East: " << n_east << " current strip is " << next ;
        
        vout.push_back(next);
        ++n_east;
        if (n_east == preshSeededNstr_) break;
     }
     // west road
     int n_west= 0;
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the West " ;
     
     theESNav.home();
     while ( ((next=theESNav.west()) != ESDetId(0) && next != strip )) {
        LogTrace("EcalClusters") << "findPi0Road: West: " << n_west << " current strip is " << next ;
        
        vout.push_back(next);
        ++n_west;
        if (n_west == preshSeededNstr_) break;
     }
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 1-st plane is " << n_east+n_west ;
     
  }
  else if (plane == 2) {
    // north road
    int n_north= 0;
    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the North " ;
    
    while ( ((next=theESNav.north()) != ESDetId(0) && next != strip) ) {
       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: North: " << n_north << " current strip is " << next;
       
       vout.push_back(next);
       ++n_north;
       if (n_north == preshSeededNstr_) break;
    }
    // south road
    int n_south= 0;
    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the South " ;
    
    theESNav.home();
    while ( ((next=theESNav.south()) != ESDetId(0) && next != strip) ) {
       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: South: " << n_south << " current strip is " << next ;
       
       vout.push_back(next);
       ++n_south;
       if (n_south == preshSeededNstr_) break;
    }
    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 2-nd plane is " << n_south+n_north ;
    
  }
  else {
    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Wrong plane number, null cluster will be returned! " ;
    
  } // end of if

  theESNav.home();
}


//===================================================================
// EndcapPiZeroDiscriminatorAlgo::readWeightFile(...), a method that reads the weigths of the NN
// INPUT: Weights_file
// OUTPUT: I_H_Weight, H_Thresh, H_O_Weight, O_Thresh arrays
//===================================================================
void EndcapPiZeroDiscriminatorAlgo::readWeightFile(const char *Weights_file, int&Layers, int& Indim, int& Hidden, int& Outdim){
   FILE *weights =nullptr;

   char line[80];

   bool checkinit=false;
// Open the weights file, generated by jetnet, and read
// in the nodes and weights
//*******************************************************
  weights = fopen(Weights_file, "r");
  LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: I opeded the Weights file  = " << Weights_file ;
  if(weights == nullptr) {
    throw cms::Exception("MissingWeightFile")<<"Could not open the weights file: "<<Weights_file;
  }

  const auto I_H_W_offset = I_H_Weight_all.size();
  const auto H_O_W_offset = H_O_Weight_all.size();
  const auto H_T_offset = H_Thresh_all.size();
  const auto O_T_offset = O_Thresh_all.size();

  while( !feof(weights) ){
    fscanf(weights, "%s", line);
    if (line[0] == 'A') { //Read in ANN nodes: Layers, input , Hidden, Output
       fscanf(weights, "%d", &Layers);   // # of NN Layers  used
       fscanf(weights, "%d", &Indim);    // # of Inputs actually used
       fscanf(weights, "%d", &Hidden);   // # of hidden nodes
       fscanf(weights, "%d", &Outdim);   // # of output nodes

           I_H_Weight_all.resize(I_H_W_offset+Indim*Hidden);
           H_Thresh_all.resize(H_T_offset+Hidden);
           H_O_Weight_all.resize(H_O_W_offset+Hidden*Outdim);
           O_Thresh_all.resize(O_T_offset+Outdim);
       checkinit=true;
    }else if (line[0] == 'B') { // read in weights between hidden and intput nodes
        assert(checkinit);
        for (int i = 0; i<Indim; i++){
                for (int j = 0; j<Hidden; j++){
                    fscanf(weights, "%f", &I_H_Weight_all[I_H_W_offset+i*Hidden+j]);
                }
        }
    }else if (line[0] == 'C'){   // Read in the thresholds for hidden nodes
        assert(checkinit);
        for (int i = 0; i<Hidden; i++){
        fscanf(weights, "%f", &H_Thresh_all[H_T_offset+i]);
        }
    }else if (line[0] == 'D'){ // read in weights between hidden and output nodes
        assert(checkinit);
        for (int i = 0; i<Hidden*Outdim; i++){
        fscanf(weights, "%f", &H_O_Weight_all[H_O_W_offset+i]);
        }
    }else if (line[0] == 'E'){ // read in the threshold for the output nodes
        assert(checkinit);
        for (int i = 0; i<Outdim; i++){
        fscanf(weights, "%f", &O_Thresh_all[O_T_offset+i]);

        }
        }
    else{edm::LogError("EEPi0Discrim")<< "EndcapPiZeroDiscriminatorAlgo: Not a Net file of Corrupted Net file " << endl;
        }
   }
   fclose(weights);
}

//=====================================================================================
// EndcapPiZeroDiscriminatorAlgo::getNNoutput(int sel_wfile), a method that calculated the NN output
// INPUT: sel_wfile -> Weight file selection
// OUTPUT : nnout -> the NN output
//=====================================================================================

float EndcapPiZeroDiscriminatorAlgo::getNNoutput(int sel_wfile, int Layers, int Indim, int Hidden, int Outdim, int barrelstart) const
{
 float nnout=0.0;
 int mij;

 std::vector<float> I_SUM(size_t(Hidden), 0.0);
 std::vector<float> OUT(size_t(Outdim), 0.0);

 for (int h = 0; h<Hidden; h++){
     mij = h - Hidden;
     for (int i = 0; i<Indim; i++){
         mij = mij + Hidden;
         I_SUM[h] += I_H_Weight_all[mij+sel_wfile*Indim*Hidden + barrelstart*Nfiles_EE*EE_Indim*EE_Hidden] * input_var[i];
     }
     I_SUM[h] += H_Thresh_all[h+sel_wfile*Hidden + barrelstart*Nfiles_EE*EE_Hidden];
     for (int o1 = 0; o1<Outdim; o1++) {
        OUT[o1] += H_O_Weight_all[barrelstart*Nfiles_EE*EE_Outdim*EE_Hidden + h*Outdim+o1 + sel_wfile*Outdim*Hidden]*Activation_fun(I_SUM[h]); 

     }
 }
 for (int o2 = 0; o2<Outdim; o2++){      
            OUT[o2] += O_Thresh_all[barrelstart*Nfiles_EE*EE_Outdim + o2 + sel_wfile*Outdim]; 
  }
  nnout = Activation_fun(OUT[0]);
  LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: getNNoutput :: -> NNout = " <<  nnout ;
  

   
  return (nnout);
}


float EndcapPiZeroDiscriminatorAlgo::Activation_fun(float SUM) const{
        return( 1.0 / ( 1.0 + exp(-2.0*SUM) ) );
}
//=====================================================================================
// EndcapPiZeroDiscriminatorAlgo::calculateNNInputVariables(...), a method that calculates the 25 input variables
// INPUTS:
// vph1 -> vector of the stip energies in 1st Preshower plane
// vph2 -> vector of the stip energies in 2nd Preshower plane
// pS1_max -> E1
// pS9_max -> E9
// pS25_max -> E25
// OUTPUT: 
// input_var[25] -> the 25 input to the NN variables array
//=====================================================================================
bool EndcapPiZeroDiscriminatorAlgo::calculateNNInputVariables(vector<float>& vph1, vector<float>& vph2, 
                                          float pS1_max, float pS9_max, float pS25_max, int EScorr)
{
   input_var.resize(EE_Indim);
   bool valid_NNinput = true;

   /*   
   for(int i = 0; i<11;i++) {
   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Energies of the Preshower Strips in X plane = " << vph1[i] ;
   }
   
   for(int i = 0; i<11;i++) {
   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Energies of the Preshower Strips in Y plane = " << vph2[i] ;
   } 
   */

   // check if all Preshower info is availabla - If NOT use remaning info
   for(int k = 0; k<11; k++) {
      if(vph1[k] < 0 ) {
      
      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " << k << " of X plane Do not exists" ;
     
         vph1[k] = 0.0;
      }      
      if(vph2[k] < 0 ) { 
    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " << k << " of Y plane Do not exists" ;   
        
        vph2[k] = 0.0;
      }
   }   

   /*
   for(int i = 0; i<11;i++) {  
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: After: Energies of the Preshower Strips in X plane = " << vph1[i] ;
   }

   for(int i = 0; i<11;i++) {
 
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: After: Energies of the Preshower Strips in Y plane = " << vph2[i] ;
   }
   */


// FIRST : Produce the 22 NN variables related with the Preshower 
// --------------------------------------------------------------
// New normalization of the preshower strip energies Aris 8/11/2004
   for(int kk=0;kk<11;kk++){
     input_var[kk] = fabs(vph1[kk]/0.01);
     input_var[kk + 11] = fabs(vph2[kk]/0.02);       
     if(input_var[kk] < 0.0001) input_var[kk] = 0.;
     if(input_var[kk + 11] < 0.0001) input_var[kk + 11] = 0.;
   }
   input_var[0] = fabs(input_var[0]/2.); 
   input_var[1] = fabs(input_var[1]/2.); 
   input_var[6] = fabs(input_var[6]/2.); 
   input_var[11] = fabs(input_var[11]/2.); 
   input_var[12] = fabs(input_var[12]/2.); 
   input_var[17] = fabs(input_var[17]/2.); 
   
// correction for version > CMSSW_3_1_0_pre5 where extra enegry is given to the ES strips 
// Aris 18/5/2009   
   if( EScorr == 1) { 
     input_var[0] -= 0.05;
     input_var[1] -= 0.035;
     input_var[2] -= 0.035;
     input_var[3] -= 0.02;
     input_var[4] -= 0.015;
     input_var[5] -= 0.0075;
     input_var[6] -= 0.035;
     input_var[7] -= 0.035;
     input_var[8] -= 0.02;
     input_var[9] -= 0.015;
     input_var[10] -= 0.0075;
    
     input_var[11] -= 0.05;
     input_var[12] -= 0.035;
     input_var[13] -= 0.035;
     input_var[14] -= 0.02;
     input_var[15] -= 0.015;
     input_var[16] -= 0.0075;
     input_var[17] -= 0.035;
     input_var[18] -= 0.035;
     input_var[19] -= 0.02;
     input_var[20] -= 0.015;
     input_var[21] -= 0.0075;   
     
     for(int kk1=0;kk1<22;kk1++){
       if(input_var[kk1] < 0 ) input_var[kk1] = 0.0;   
     }
   }
// SECOND: Take the final NN variable related to the ECAL
// -----------------------------------------------
   float ECAL_norm_factor = 500.;
   if(pS25_max>500&&pS25_max<=1000) ECAL_norm_factor = 1000;
   if(pS25_max>1000) ECAL_norm_factor = 7000;

   input_var[22] = pS1_max/ECAL_norm_factor;
   input_var[23] = pS9_max/ECAL_norm_factor;
   input_var[24] = pS25_max/ECAL_norm_factor;

   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S1/ECAL_norm_factor = " << input_var[22];
   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S9/ECAL_norm_factor = " << input_var[23];
   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S25/ECAL_norm_factor = " << input_var[24] ;

   for(int i=0;i<EE_Indim;i++){
     if(input_var[i] > 1.0e+00) {
       valid_NNinput = false;
       break;
     }
   }
   
   LogTrace("EcalClusters") << " valid_NNinput = " << valid_NNinput ;
   

   return valid_NNinput;
}


//=====================================================================================
// EndcapPiZeroDiscriminatorAlgo::calculateBarrelNNInputVariables(...), a method that calculates
// the 12 barrel NN input
// OUTPUT:
// input_var[12] -> the 12 input to the barrel NN variables array
//=====================================================================================

void EndcapPiZeroDiscriminatorAlgo::calculateBarrelNNInputVariables(float et, double s1, double s9, double s25,
                                    double m2, double cee, double cep,double cpp,
                                    double s4, double s6, double ratio,
                                    double xcog, double ycog)
{
  input_var.resize(EB_Indim);

  double lam, lam1, lam2;

  if(xcog < 0.) { 
    input_var[0] = -xcog/s25;
  } else { 
    input_var[0] = xcog/s25;
  }    

  input_var[1] = cee/0.0004;

  if(cpp<.001) {
    input_var[2] = cpp/.001;
  } else  { 
    input_var[2] = 0.;
  }

  if(s9!=0.) {
    input_var[3] = s1/s9; 
    input_var[8] = s6/s9; 
    input_var[10] = (m2+s1)/s9;
  }
  else {
    input_var[3] = 0.;
    input_var[8] = 0.;
    input_var[10] = 0.;
  }

  if(s25-s1>0.) {
    input_var[4] = (s9-s1)/(s25-s1); 
  } else {
    input_var[4] = 0.;
  }  

  if(s25>0.) {
    input_var[5] = s4/s25; 
  }  else {
    input_var[5] = 0.;
  } 

  if(ycog < 0.) {
    input_var[6] = -ycog/s25; 
  } else {
    input_var[6] = ycog/s25;
  } 

  input_var[7] = ratio;

  lam=sqrt((cee -cpp)*(cee -cpp)+4*cep*cep);
  lam1=(cee + cpp + lam)/2;
  lam2=(cee + cpp - lam)/2;

  if(lam1 == 0) {
    input_var[9] = .0;
  } else {
    input_var[9] = lam2/lam1;
  }
  if(s4!=0.) {
    input_var[11] = (m2+s1)/s4;
  } else {
    input_var[11] = 0.;
  }

}


//=====================================================================================
// EndcapPiZeroDiscriminatorAlgo::GetNNOutput(...), a method that calculates the NNoutput
// INPUTS: Super Cluster Energy
// OUTPUTS : NNoutput
//=====================================================================================
float EndcapPiZeroDiscriminatorAlgo::GetNNOutput(float EE_Et)
{
     float nnout = -1;
// Print the  NN input variables that are related to the Preshower + ECAL
// ------------------------------------------------------------------------
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo::GetNNoutput :nn_invar_presh = " ;

     LogTrace("EcalClusters").log([&](auto& lt) {
       for(auto const v: input_var) {
         lt << v << " " ;
       }
     });
     LogTrace("EcalClusters") << " " ;
    
     // select the appropriate Weigth file
     int sel_wfile;
     if(EE_Et<25.0)                     {sel_wfile = 0;}
     else if(EE_Et>=25.0 && EE_Et<35.0) {sel_wfile = 1;}
     else if(EE_Et>=35.0 && EE_Et<45.0) {sel_wfile = 2;}
     else if(EE_Et>=45.0 && EE_Et<55.0) {sel_wfile = 3;}
     else                               {sel_wfile = 4;} 

     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Et_SC = " << EE_Et << " and I select Weight file Number = " << sel_wfile ;

     
     nnout = getNNoutput(sel_wfile, EE_Layers, EE_Indim, EE_Hidden, EE_Outdim, 0); // calculate the nnoutput for the given ECAL object
     
     LogTrace("EcalClusters") <<  "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
     
    
   return nnout;
}


//=====================================================================================
// EndcapPiZeroDiscriminatorAlgo::GetBarrelNNOutput(...), a method that calculates the barrel NNoutput
// INPUTS: Super Cluster Energy
// OUTPUTS : NNoutput
//=====================================================================================
float EndcapPiZeroDiscriminatorAlgo::GetBarrelNNOutput(float EB_Et)
{

     float nnout = -1;
// Print the  NN input variables that are related to the ECAL Barrel
// ------------------------------------------------------------------------
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo::GetBarrelNNoutput :nn_invar_presh = " ;

     LogTrace("EcalCluster").log([&](auto& lt) {
       for(auto const v: input_var) {
         lt << v << " " ;  
       }
      });
     LogTrace("EcalClusters") << " " ;
     
     // select the appropriate Weigth file
     int sel_wfile;
     if(EB_Et<25.0)                     {sel_wfile = 0;}
     else if(EB_Et>=25.0 && EB_Et<35.0) {sel_wfile = 1;}
     else if(EB_Et>=35.0 && EB_Et<45.0) {sel_wfile = 2;}
     else if(EB_Et>=45.0 && EB_Et<55.0) {sel_wfile = 3;}
     else                               {sel_wfile = 4;}
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: E_SC = " << EB_Et << " and I select Weight file Number = " << sel_wfile ;
    
     nnout = getNNoutput(sel_wfile, EB_Layers, EB_Indim, EB_Hidden, EB_Outdim, 1); // calculate the nnoutput for the given ECAL object
     
     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
    
   return nnout;
}