EgammaClusterAlgos/src/EndcapPiZeroDiscriminatorAlgo.cc

0001
0002 #include "RecoEcal/EgammaClusterAlgos/interface/EndcapPiZeroDiscriminatorAlgo.h"
0003 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/FileInPath.h"
0008 #include <fstream>
0009 #include <iostream>
0010
0011 using namespace std;
0012
0013 namespace {
0014   constexpr int Nfiles_EE = 5;
0015 }
0016
0017 EndcapPiZeroDiscriminatorAlgo::EndcapPiZeroDiscriminatorAlgo(double stripEnergyCut, int nStripCut, const string& path)
0018     : preshStripEnergyCut_(stripEnergyCut), preshSeededNstr_(nStripCut) {
0019   // Read all Weight files
0020   constexpr std::array<char const*, Nfiles_EE> file_pt{{"20", "30", "40", "50", "60"}};
0021   constexpr std::array<char const*, 5> file_barrel_pt{{"20", "30", "40", "50", "60"}};
0022
0023   for (auto ptName : file_pt) {
0024     string nn_paterns_file = "endcapPiZeroDiscriminatorWeights_et";
0025     nn_paterns_file += ptName;
0026     nn_paterns_file += ".wts";
0027     edm::FileInPath WFile(path + nn_paterns_file);
0028     readWeightFile(WFile.fullPath().c_str(), EE_Layers, EE_Indim, EE_Hidden, EE_Outdim);  // read the weights' file
0029   }
0030
0031   for (auto ptName : file_barrel_pt) {
0032     string nn_paterns_file = "barrelPiZeroDiscriminatorWeights_et";
0033     nn_paterns_file += ptName;
0034     nn_paterns_file += ".wts";
0035     edm::FileInPath WFile(path + nn_paterns_file);
0036     readWeightFile(WFile.fullPath().c_str(), EB_Layers, EB_Indim, EB_Hidden, EB_Outdim);  // read the weights' file
0037   }
0038 }
0039
0040 vector<float> EndcapPiZeroDiscriminatorAlgo::findPreshVector(ESDetId strip,
0041                                                              RecHitsMap* rechits_map,
0042                                                              CaloSubdetectorTopology* topology_p) {
0043   vector<float> vout_stripE;
0044
0045   // skip if rechits_map contains no hits
0046   if (rechits_map->empty()) {
0047     edm::LogWarning("EndcapPiZeroDiscriminatorAlgo") << "RecHitsMap has size 0.";
0048     return vout_stripE;
0049   }
0050
0051   vout_stripE.clear();
0052
0053   vector<ESDetId> road_2d;
0054   road_2d.clear();
0055
0056   int plane = strip.plane();
0057
0058   LogTrace("EcalClusters")
0059       << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: Preshower Seeded Algorithm - looking for clusters"
0060       << "n"
0061       << "findPreshVectors: Preshower is intersected at strip " << strip.strip() << ", at plane " << plane;
0062
0063   if (strip == ESDetId(0)) {  //works in case of no intersected strip found
0064     for (int i = 0; i < 11; i++) {
0065       vout_stripE.push_back(-100.);
0066     }
0067   }
0068
0069   // Add to the road the central strip
0070   road_2d.push_back(strip);
0071
0072   //Make a navigator, and set it to the strip cell.
0073   EcalPreshowerNavigator navigator(strip, topology_p);
0074   navigator.setHome(strip);
0075   //search for neighbours in the central road
0076   findPi0Road(strip, navigator, plane, road_2d);
0077
0078   LogTrace("EcalClusters")
0079       << "EndcapPiZeroDiscriminatorAlgo:findPreshVectors: Total number of strips in the central road: "
0080       << road_2d.size();
0081
0082   // Find the energy of each strip
0083   RecHitsMap::iterator final_strip = rechits_map->end();
0084   // very dangerous, added a protection on the rechits_map->size()
0085   // at the beginning of the method
0086   final_strip--;
0087   ESDetId last_stripID = final_strip->first;
0088
0089   vector<ESDetId>::iterator itID;
0090   for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
0091     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: ID = " << *itID;
0092
0093     float E = 0.;
0094     RecHitsMap::iterator strip_it = rechits_map->find(*itID);
0095     if (goodPi0Strip(strip_it, last_stripID)) {  // continue if strip not found in rechit_map
0096       E = strip_it->second.energy();
0097     }
0098     vout_stripE.push_back(E);
0099     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: E = " << E;
0100   }
0101
0102   // ***ML beg***
0103   // vector of size=11, content of vout_stripE is copied into vout_ElevenStrips_Energy
0104   // to avoid problem in case number of strips is less than 11
0105   vector<float> vout_ElevenStrips_Energy;
0106   vout_ElevenStrips_Energy.reserve(11);
0107   for (int i = 0; i < 11; i++) {
0108     vout_ElevenStrips_Energy.push_back(0.);
0109   }
0110
0111   for (unsigned int i = 0; i < vout_stripE.size(); i++) {
0112     vout_ElevenStrips_Energy[i] = vout_stripE.at(i);
0113   }
0114
0115   //return vout_stripE;
0116   return vout_ElevenStrips_Energy;
0117   // ***ML end***
0118 }
0119
0120 // returns true if the candidate strip fulfills the requirements to be added to the cluster:
0121 bool EndcapPiZeroDiscriminatorAlgo::goodPi0Strip(RecHitsMap::iterator candidate_it, ESDetId lastID) {
0122   RecHitsMap::iterator candidate_tmp = candidate_it;
0123   candidate_tmp--;
0124
0125   // crystal should not be included...
0126   if (candidate_tmp->first == lastID)  // ...if it corresponds to a hit
0127   {
0128     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: goodPi0Strip No such a strip in rechits_map ";
0129     return false;
0130   } else if (candidate_it->second.energy() <= preshStripEnergyCut_)  // ...if it has a negative or zero energy
0131   {
0132     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: goodPi0Strip Strip energy "
0133                              << candidate_it->second.energy() << " is below threshold ";
0134     return false;
0135   }
0136
0137   return true;
0138 }
0139
0140 // find strips in the road of size +/- preshSeededNstr_ from the central strip
0141 void EndcapPiZeroDiscriminatorAlgo::findPi0Road(ESDetId strip,
0142                                                 EcalPreshowerNavigator& theESNav,
0143                                                 int plane,
0144                                                 vector<ESDetId>& vout) {
0145   if (strip == ESDetId(0))
0146     return;
0147   ESDetId next;
0148   theESNav.setHome(strip);
0149   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: starts from strip " << strip;
0150
0151   if (plane == 1) {
0152     // east road
0153     int n_east = 0;
0154     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the East ";
0155
0156     while (((next = theESNav.east()) != ESDetId(0) && next != strip)) {
0157       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: East: " << n_east << " current strip is "
0158                                << next;
0159
0160       vout.push_back(next);
0161       ++n_east;
0162       if (n_east == preshSeededNstr_)
0163         break;
0164     }
0165     // west road
0166     int n_west = 0;
0167     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the West ";
0168
0169     theESNav.home();
0170     while (((next = theESNav.west()) != ESDetId(0) && next != strip)) {
0171       LogTrace("EcalClusters") << "findPi0Road: West: " << n_west << " current strip is " << next;
0172
0173       vout.push_back(next);
0174       ++n_west;
0175       if (n_west == preshSeededNstr_)
0176         break;
0177     }
0178     LogTrace("EcalClusters")
0179         << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 1-st plane is "
0180         << n_east + n_west;
0181
0182   } else if (plane == 2) {
0183     // north road
0184     int n_north = 0;
0185     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the North ";
0186
0187     while (((next = theESNav.north()) != ESDetId(0) && next != strip)) {
0188       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: North: " << n_north
0189                                << " current strip is " << next;
0190
0191       vout.push_back(next);
0192       ++n_north;
0193       if (n_north == preshSeededNstr_)
0194         break;
0195     }
0196     // south road
0197     int n_south = 0;
0198     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the South ";
0199
0200     theESNav.home();
0201     while (((next = theESNav.south()) != ESDetId(0) && next != strip)) {
0202       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: South: " << n_south
0203                                << " current strip is " << next;
0204
0205       vout.push_back(next);
0206       ++n_south;
0207       if (n_south == preshSeededNstr_)
0208         break;
0209     }
0210     LogTrace("EcalClusters")
0211         << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 2-nd plane is "
0212         << n_south + n_north;
0213
0214   } else {
0215     LogTrace("EcalClusters")
0216         << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Wrong plane number, null cluster will be returned! ";
0217
0218   }  // end of if
0219
0220   theESNav.home();
0221 }
0222
0223 //===================================================================
0224 // EndcapPiZeroDiscriminatorAlgo::readWeightFile(...), a method that reads the weigths of the NN
0225 // INPUT: Weights_file
0226 // OUTPUT: I_H_Weight, H_Thresh, H_O_Weight, O_Thresh arrays
0227 //===================================================================
0228 void EndcapPiZeroDiscriminatorAlgo::readWeightFile(
0229     const char* Weights_file, int& Layers, int& Indim, int& Hidden, int& Outdim) {
0230   FILE* weights = nullptr;
0231
0232   char line[80];
0233
0234   bool checkinit = false;
0235   // Open the weights file, generated by jetnet, and read
0236   // in the nodes and weights
0237   //*******************************************************
0238   weights = fopen(Weights_file, "r");
0239   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: I opeded the Weights file  = " << Weights_file;
0240   if (weights == nullptr) {
0241     throw cms::Exception("MissingWeightFile") << "Could not open the weights file: " << Weights_file;
0242   }
0243
0244   const auto I_H_W_offset = I_H_Weight_all.size();
0245   const auto H_O_W_offset = H_O_Weight_all.size();
0246   const auto H_T_offset = H_Thresh_all.size();
0247   const auto O_T_offset = O_Thresh_all.size();
0248
0249   while (!feof(weights)) {
0250     fscanf(weights, "%s", line);
0251     if (line[0] == 'A') {              //Read in ANN nodes: Layers, input , Hidden, Output
0252       fscanf(weights, "%d", &Layers);  // # of NN Layers  used
0253       fscanf(weights, "%d", &Indim);   // # of Inputs actually used
0254       fscanf(weights, "%d", &Hidden);  // # of hidden nodes
0255       fscanf(weights, "%d", &Outdim);  // # of output nodes
0256
0257       I_H_Weight_all.resize(I_H_W_offset + Indim * Hidden);
0258       H_Thresh_all.resize(H_T_offset + Hidden);
0259       H_O_Weight_all.resize(H_O_W_offset + Hidden * Outdim);
0260       O_Thresh_all.resize(O_T_offset + Outdim);
0261       checkinit = true;
0262     } else if (line[0] == 'B') {  // read in weights between hidden and intput nodes
0263       assert(checkinit);
0264       for (int i = 0; i < Indim; i++) {
0265         for (int j = 0; j < Hidden; j++) {
0266           fscanf(weights, "%f", &I_H_Weight_all[I_H_W_offset + i * Hidden + j]);
0267         }
0268       }
0269     } else if (line[0] == 'C') {  // Read in the thresholds for hidden nodes
0270       assert(checkinit);
0271       for (int i = 0; i < Hidden; i++) {
0272         fscanf(weights, "%f", &H_Thresh_all[H_T_offset + i]);
0273       }
0274     } else if (line[0] == 'D') {  // read in weights between hidden and output nodes
0275       assert(checkinit);
0276       for (int i = 0; i < Hidden * Outdim; i++) {
0277         fscanf(weights, "%f", &H_O_Weight_all[H_O_W_offset + i]);
0278       }
0279     } else if (line[0] == 'E') {  // read in the threshold for the output nodes
0280       assert(checkinit);
0281       for (int i = 0; i < Outdim; i++) {
0282         fscanf(weights, "%f", &O_Thresh_all[O_T_offset + i]);
0283       }
0284     } else {
0285       edm::LogError("EEPi0Discrim") << "EndcapPiZeroDiscriminatorAlgo: Not a Net file of Corrupted Net file " << endl;
0286     }
0287   }
0288   fclose(weights);
0289 }
0290
0291 //=====================================================================================
0292 // EndcapPiZeroDiscriminatorAlgo::getNNoutput(int sel_wfile), a method that calculated the NN output
0293 // INPUT: sel_wfile -> Weight file selection
0294 // OUTPUT : nnout -> the NN output
0295 //=====================================================================================
0296
0297 float EndcapPiZeroDiscriminatorAlgo::getNNoutput(
0298     int sel_wfile, int Layers, int Indim, int Hidden, int Outdim, int barrelstart) const {
0299   float nnout = 0.0;
0300   int mij;
0301
0302   std::vector<float> I_SUM(size_t(Hidden), 0.0);
0303   std::vector<float> OUT(size_t(Outdim), 0.0);
0304
0305   for (int h = 0; h < Hidden; h++) {
0306     mij = h - Hidden;
0307     for (int i = 0; i < Indim; i++) {
0308       mij = mij + Hidden;
0309       I_SUM[h] += I_H_Weight_all[mij + sel_wfile * Indim * Hidden + barrelstart * Nfiles_EE * EE_Indim * EE_Hidden] *
0310                   input_var[i];
0311     }
0312     I_SUM[h] += H_Thresh_all[h + sel_wfile * Hidden + barrelstart * Nfiles_EE * EE_Hidden];
0313     for (int o1 = 0; o1 < Outdim; o1++) {
0314       OUT[o1] += H_O_Weight_all[barrelstart * Nfiles_EE * EE_Outdim * EE_Hidden + h * Outdim + o1 +
0315                                 sel_wfile * Outdim * Hidden] *
0316                  Activation_fun(I_SUM[h]);
0317     }
0318   }
0319   for (int o2 = 0; o2 < Outdim; o2++) {
0320     OUT[o2] += O_Thresh_all[barrelstart * Nfiles_EE * EE_Outdim + o2 + sel_wfile * Outdim];
0321   }
0322   nnout = Activation_fun(OUT[0]);
0323   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: getNNoutput :: -> NNout = " << nnout;
0324
0325   return (nnout);
0326 }
0327
0328 float EndcapPiZeroDiscriminatorAlgo::Activation_fun(float SUM) const { return (1.0 / (1.0 + exp(-2.0 * SUM))); }
0329 //=====================================================================================
0330 // EndcapPiZeroDiscriminatorAlgo::calculateNNInputVariables(...), a method that calculates the 25 input variables
0331 // INPUTS:
0332 // vph1 -> vector of the stip energies in 1st Preshower plane
0333 // vph2 -> vector of the stip energies in 2nd Preshower plane
0334 // pS1_max -> E1
0335 // pS9_max -> E9
0336 // pS25_max -> E25
0337 // OUTPUT:
0338 // input_var[25] -> the 25 input to the NN variables array
0339 //=====================================================================================
0340 bool EndcapPiZeroDiscriminatorAlgo::calculateNNInputVariables(
0341     vector<float>& vph1, vector<float>& vph2, float pS1_max, float pS9_max, float pS25_max, int EScorr) {
0342   input_var.resize(EE_Indim);
0343   bool valid_NNinput = true;
0344
0345   /*
0346    for(int i = 0; i<11;i++) {
0347    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Energies of the Preshower Strips in X plane = " << vph1[i] ;
0348    }
0349
0350    for(int i = 0; i<11;i++) {
0351    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Energies of the Preshower Strips in Y plane = " << vph2[i] ;
0352    }
0353    */
0354
0355   // check if all Preshower info is availabla - If NOT use remaning info
0356   for (int k = 0; k < 11; k++) {
0357     if (vph1[k] < 0) {
0358       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " << k
0359                                << " of X plane Do not exists";
0360
0361       vph1[k] = 0.0;
0362     }
0363     if (vph2[k] < 0) {
0364       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " << k
0365                                << " of Y plane Do not exists";
0366
0367       vph2[k] = 0.0;
0368     }
0369   }
0370
0371   /*
0372    for(int i = 0; i<11;i++) {
0373      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: After: Energies of the Preshower Strips in X plane = " << vph1[i] ;
0374    }
0375
0376    for(int i = 0; i<11;i++) {
0377
0378      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: After: Energies of the Preshower Strips in Y plane = " << vph2[i] ;
0379    }
0380    */
0381
0382   // FIRST : Produce the 22 NN variables related with the Preshower
0383   // --------------------------------------------------------------
0384   // New normalization of the preshower strip energies Aris 8/11/2004
0385   for (int kk = 0; kk < 11; kk++) {
0386     input_var[kk] = fabs(vph1[kk] / 0.01);
0387     input_var[kk + 11] = fabs(vph2[kk] / 0.02);
0388     if (input_var[kk] < 0.0001)
0389       input_var[kk] = 0.;
0390     if (input_var[kk + 11] < 0.0001)
0391       input_var[kk + 11] = 0.;
0392   }
0393   input_var[0] = fabs(input_var[0] / 2.);
0394   input_var[1] = fabs(input_var[1] / 2.);
0395   input_var[6] = fabs(input_var[6] / 2.);
0396   input_var[11] = fabs(input_var[11] / 2.);
0397   input_var[12] = fabs(input_var[12] / 2.);
0398   input_var[17] = fabs(input_var[17] / 2.);
0399
0400   // correction for version > CMSSW_3_1_0_pre5 where extra enegry is given to the ES strips
0401   // Aris 18/5/2009
0402   if (EScorr == 1) {
0403     input_var[0] -= 0.05;
0404     input_var[1] -= 0.035;
0405     input_var[2] -= 0.035;
0406     input_var[3] -= 0.02;
0407     input_var[4] -= 0.015;
0408     input_var[5] -= 0.0075;
0409     input_var[6] -= 0.035;
0410     input_var[7] -= 0.035;
0411     input_var[8] -= 0.02;
0412     input_var[9] -= 0.015;
0413     input_var[10] -= 0.0075;
0414
0415     input_var[11] -= 0.05;
0416     input_var[12] -= 0.035;
0417     input_var[13] -= 0.035;
0418     input_var[14] -= 0.02;
0419     input_var[15] -= 0.015;
0420     input_var[16] -= 0.0075;
0421     input_var[17] -= 0.035;
0422     input_var[18] -= 0.035;
0423     input_var[19] -= 0.02;
0424     input_var[20] -= 0.015;
0425     input_var[21] -= 0.0075;
0426
0427     for (int kk1 = 0; kk1 < 22; kk1++) {
0428       if (input_var[kk1] < 0)
0429         input_var[kk1] = 0.0;
0430     }
0431   }
0432   // SECOND: Take the final NN variable related to the ECAL
0433   // -----------------------------------------------
0434   float ECAL_norm_factor = 500.;
0435   if (pS25_max > 500 && pS25_max <= 1000)
0436     ECAL_norm_factor = 1000;
0437   if (pS25_max > 1000)
0438     ECAL_norm_factor = 7000;
0439
0440   input_var[22] = pS1_max / ECAL_norm_factor;
0441   input_var[23] = pS9_max / ECAL_norm_factor;
0442   input_var[24] = pS25_max / ECAL_norm_factor;
0443
0444   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S1/ECAL_norm_factor = " << input_var[22];
0445   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S9/ECAL_norm_factor = " << input_var[23];
0446   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S25/ECAL_norm_factor = " << input_var[24];
0447
0448   for (int i = 0; i < EE_Indim; i++) {
0449     if (input_var[i] > 1.0e+00) {
0450       valid_NNinput = false;
0451       break;
0452     }
0453   }
0454
0455   LogTrace("EcalClusters") << " valid_NNinput = " << valid_NNinput;
0456
0457   return valid_NNinput;
0458 }
0459
0460 //=====================================================================================
0461 // EndcapPiZeroDiscriminatorAlgo::calculateBarrelNNInputVariables(...), a method that calculates
0462 // the 12 barrel NN input
0463 // OUTPUT:
0464 // input_var[12] -> the 12 input to the barrel NN variables array
0465 //=====================================================================================
0466
0467 void EndcapPiZeroDiscriminatorAlgo::calculateBarrelNNInputVariables(float et,
0468                                                                     double s1,
0469                                                                     double s9,
0470                                                                     double s25,
0471                                                                     double m2,
0472                                                                     double cee,
0473                                                                     double cep,
0474                                                                     double cpp,
0475                                                                     double s4,
0476                                                                     double s6,
0477                                                                     double ratio,
0478                                                                     double xcog,
0479                                                                     double ycog) {
0480   input_var.resize(EB_Indim);
0481
0482   double lam, lam1, lam2;
0483
0484   if (xcog < 0.) {
0485     input_var[0] = -xcog / s25;
0486   } else {
0487     input_var[0] = xcog / s25;
0488   }
0489
0490   input_var[1] = cee / 0.0004;
0491
0492   if (cpp < .001) {
0493     input_var[2] = cpp / .001;
0494   } else {
0495     input_var[2] = 0.;
0496   }
0497
0498   if (s9 != 0.) {
0499     input_var[3] = s1 / s9;
0500     input_var[8] = s6 / s9;
0501     input_var[10] = (m2 + s1) / s9;
0502   } else {
0503     input_var[3] = 0.;
0504     input_var[8] = 0.;
0505     input_var[10] = 0.;
0506   }
0507
0508   if (s25 - s1 > 0.) {
0509     input_var[4] = (s9 - s1) / (s25 - s1);
0510   } else {
0511     input_var[4] = 0.;
0512   }
0513
0514   if (s25 > 0.) {
0515     input_var[5] = s4 / s25;
0516   } else {
0517     input_var[5] = 0.;
0518   }
0519
0520   if (ycog < 0.) {
0521     input_var[6] = -ycog / s25;
0522   } else {
0523     input_var[6] = ycog / s25;
0524   }
0525
0526   input_var[7] = ratio;
0527
0528   lam = sqrt((cee - cpp) * (cee - cpp) + 4 * cep * cep);
0529   lam1 = (cee + cpp + lam) / 2;
0530   lam2 = (cee + cpp - lam) / 2;
0531
0532   if (lam1 == 0) {
0533     input_var[9] = .0;
0534   } else {
0535     input_var[9] = lam2 / lam1;
0536   }
0537   if (s4 != 0.) {
0538     input_var[11] = (m2 + s1) / s4;
0539   } else {
0540     input_var[11] = 0.;
0541   }
0542 }
0543
0544 //=====================================================================================
0545 // EndcapPiZeroDiscriminatorAlgo::GetNNOutput(...), a method that calculates the NNoutput
0546 // INPUTS: Super Cluster Energy
0547 // OUTPUTS : NNoutput
0548 //=====================================================================================
0549 float EndcapPiZeroDiscriminatorAlgo::GetNNOutput(float EE_Et) {
0550   float nnout = -1;
0551   // Print the  NN input variables that are related to the Preshower + ECAL
0552   // ------------------------------------------------------------------------
0553   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo::GetNNoutput :nn_invar_presh = ";
0554
0555   LogTrace("EcalClusters").log([&](auto& lt) {
0556     for (auto const v : input_var) {
0557       lt << v << " ";
0558     }
0559   });
0560   LogTrace("EcalClusters") << " ";
0561
0562   // select the appropriate Weigth file
0563   int sel_wfile;
0564   if (EE_Et < 25.0) {
0565     sel_wfile = 0;
0566   } else if (EE_Et >= 25.0 && EE_Et < 35.0) {
0567     sel_wfile = 1;
0568   } else if (EE_Et >= 35.0 && EE_Et < 45.0) {
0569     sel_wfile = 2;
0570   } else if (EE_Et >= 45.0 && EE_Et < 55.0) {
0571     sel_wfile = 3;
0572   } else {
0573     sel_wfile = 4;
0574   }
0575
0576   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Et_SC = " << EE_Et
0577                            << " and I select Weight file Number = " << sel_wfile;
0578
0579   nnout = getNNoutput(
0580       sel_wfile, EE_Layers, EE_Indim, EE_Hidden, EE_Outdim, 0);  // calculate the nnoutput for the given ECAL object
0581
0582   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout;
0583
0584   return nnout;
0585 }
0586
0587 //=====================================================================================
0588 // EndcapPiZeroDiscriminatorAlgo::GetBarrelNNOutput(...), a method that calculates the barrel NNoutput
0589 // INPUTS: Super Cluster Energy
0590 // OUTPUTS : NNoutput
0591 //=====================================================================================
0592 float EndcapPiZeroDiscriminatorAlgo::GetBarrelNNOutput(float EB_Et) {
0593   float nnout = -1;
0594   // Print the  NN input variables that are related to the ECAL Barrel
0595   // ------------------------------------------------------------------------
0596   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo::GetBarrelNNoutput :nn_invar_presh = ";
0597
0598   LogTrace("EcalCluster").log([&](auto& lt) {
0599     for (auto const v : input_var) {
0600       lt << v << " ";
0601     }
0602   });
0603   LogTrace("EcalClusters") << " ";
0604
0605   // select the appropriate Weigth file
0606   int sel_wfile;
0607   if (EB_Et < 25.0) {
0608     sel_wfile = 0;
0609   } else if (EB_Et >= 25.0 && EB_Et < 35.0) {
0610     sel_wfile = 1;
0611   } else if (EB_Et >= 35.0 && EB_Et < 45.0) {
0612     sel_wfile = 2;
0613   } else if (EB_Et >= 45.0 && EB_Et < 55.0) {
0614     sel_wfile = 3;
0615   } else {
0616     sel_wfile = 4;
0617   }
0618   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: E_SC = " << EB_Et
0619                            << " and I select Weight file Number = " << sel_wfile;
0620
0621   nnout = getNNoutput(
0622       sel_wfile, EB_Layers, EB_Indim, EB_Hidden, EB_Outdim, 1);  // calculate the nnoutput for the given ECAL object
0623
0624   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout;
0625
0626   return nnout;
0627 }