CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_0/src/RecoEcal/EgammaClusterAlgos/src/EndcapPiZeroDiscriminatorAlgo.cc

Go to the documentation of this file.
00001 
00002 #include "RecoEcal/EgammaClusterAlgos/interface/EndcapPiZeroDiscriminatorAlgo.h"
00003 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00005 
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/ParameterSet/interface/FileInPath.h"
00008 #include <fstream>
00009 #include <iostream>
00010 
00011 using namespace std;
00012 
00013 EndcapPiZeroDiscriminatorAlgo::EndcapPiZeroDiscriminatorAlgo(double stripEnergyCut, int nStripCut, const string& path) :
00014    preshStripEnergyCut_(stripEnergyCut),  preshSeededNstr_(nStripCut), pathToFiles_(path)
00015 {   
00016    
00017      // Read all Weight files
00018      Nfiles_EB = 5;
00019      Nfiles_EE = 5;
00020      string file_pt[5] = {"20","30","40","50","60"};
00021      string file_barrel_pt[5] = {"20","30","40","50","60"};   
00022 
00023      string nn_paterns_file  = "";
00024      for(int j=0;j<Nfiles_EE;j++) {
00025        nn_paterns_file = "endcapPiZeroDiscriminatorWeights_et"+file_pt[j]+".wts";
00026        edm::FileInPath WFile(pathToFiles_+nn_paterns_file);
00027        readWeightFile(WFile.fullPath().c_str()); // read the weights' file
00028 
00029        EE_Layers = Layers;
00030        EE_Indim = Indim;
00031        EE_Hidden = Hidden;
00032        EE_Outdim = Outdim;
00033           
00034        for(int i=0;i<Indim*Hidden;i++)  I_H_Weight_all.push_back(I_H_Weight[i]);
00035        for(int i=0;i<Hidden;i++)  H_Thresh_all.push_back(H_Thresh[i]);
00036        for(int i=0;i<Outdim*Hidden;i++)  H_O_Weight_all.push_back(H_O_Weight[i]);
00037        for(int i=0;i<Outdim;i++)  O_Thresh_all.push_back(O_Thresh[i]);
00038      }
00039 
00040      for(int k=0;k<Nfiles_EB;k++) {
00041        nn_paterns_file = "barrelPiZeroDiscriminatorWeights_et"+file_barrel_pt[k]+".wts";
00042        edm::FileInPath WFile(pathToFiles_+nn_paterns_file);
00043        readWeightFile(WFile.fullPath().c_str()); // read the weights' file
00044 
00045        EB_Layers = Layers;
00046        EB_Indim = Indim;
00047        EB_Hidden = Hidden;
00048        EB_Outdim = Outdim;
00049 
00050        for(int i=0;i<Indim*Hidden;i++)  I_H_Weight_all.push_back(I_H_Weight[i]);
00051        for(int i=0;i<Hidden;i++)  H_Thresh_all.push_back(H_Thresh[i]);
00052        for(int i=0;i<Outdim*Hidden;i++)  H_O_Weight_all.push_back(H_O_Weight[i]);
00053        for(int i=0;i<Outdim;i++)  O_Thresh_all.push_back(O_Thresh[i]);
00054      }
00055    delete [] I_H_Weight;
00056    delete [] H_Thresh;
00057    delete [] H_O_Weight;
00058    delete [] O_Thresh;
00059 }
00060 
00061 
00062 vector<float> EndcapPiZeroDiscriminatorAlgo::findPreshVector(ESDetId strip,  RecHitsMap *rechits_map,
00063                                                      CaloSubdetectorTopology *topology_p)
00064 {
00065   vector<float> vout_stripE;
00066 
00067   // skip if rechits_map contains no hits
00068   if ( rechits_map->size() == 0 ) {
00069           edm::LogWarning("EndcapPiZeroDiscriminatorAlgo") << "RecHitsMap has size 0.";
00070           return vout_stripE;
00071   }
00072 
00073   vout_stripE.clear();
00074 
00075   vector<ESDetId> road_2d;
00076   road_2d.clear();
00077 
00078   int plane = strip.plane();
00079 
00080   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: Preshower Seeded Algorithm - looking for clusters" << "n"<< "findPreshVectors: Preshower is intersected at strip " << strip.strip() << ", at plane " << plane ;
00081 
00082 
00083   if ( strip == ESDetId(0) ) { //works in case of no intersected strip found
00084     for(int i=0;i<11;i++) {
00085        vout_stripE.push_back(-100.);
00086     }
00087   }
00088 
00089   // Add to the road the central strip
00090   road_2d.push_back(strip);
00091 
00092   //Make a navigator, and set it to the strip cell.
00093   EcalPreshowerNavigator navigator(strip, topology_p);
00094   navigator.setHome(strip);
00095  //search for neighbours in the central road
00096   findPi0Road(strip, navigator, plane, road_2d);
00097   
00098   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo:findPreshVectors: Total number of strips in the central road: " << road_2d.size() ;
00099   
00100   // Find the energy of each strip
00101   RecHitsMap::iterator final_strip =  rechits_map->end();
00102   // very dangerous, added a protection on the rechits_map->size()
00103   // at the beginning of the method
00104   final_strip--;
00105   ESDetId last_stripID = final_strip->first;
00106 
00107   float E = 0;
00108   vector<ESDetId>::iterator itID;
00109   for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
00110     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: ID = " << *itID ;
00111   
00112     RecHitsMap::iterator strip_it = rechits_map->find(*itID);
00113     if(goodPi0Strip(strip_it,last_stripID)) { // continue if strip not found in rechit_map
00114       E = strip_it->second.energy();
00115     } else  E = 0;
00116     vout_stripE.push_back(E);
00117     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPreshVectors: E = " << E ;
00118     
00119   }
00120  
00121   // ***ML beg***
00122   // vector of size=11, content of vout_stripE is copied into vout_ElevenStrips_Energy
00123   // to avoid problem in case number of strips is less than 11
00124   vector<float> vout_ElevenStrips_Energy;
00125   for(int i=0;i<11;i++)
00126   {
00127     vout_ElevenStrips_Energy.push_back(0.);
00128   }
00129   
00130   for(unsigned int i=0;i<vout_stripE.size();i++)
00131   {
00132     vout_ElevenStrips_Energy[i] = vout_stripE.at(i);
00133   }
00134 
00135   //return vout_stripE;
00136   return vout_ElevenStrips_Energy;
00137   // ***ML end***
00138 
00139 }
00140 
00141 // returns true if the candidate strip fulfills the requirements to be added to the cluster:
00142 bool EndcapPiZeroDiscriminatorAlgo::goodPi0Strip(RecHitsMap::iterator candidate_it, ESDetId lastID)
00143 {
00144   RecHitsMap::iterator candidate_tmp = candidate_it;
00145   candidate_tmp--;
00146 
00147   if (candidate_tmp->first == lastID )
00148   LogTrace("EcalClusters") <<"EndcapPiZeroDiscriminatorAlgo: goodPi0Strip No such a strip in rechits_map " ;
00149   if (candidate_it->second.energy() <= preshStripEnergyCut_)
00150   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: goodPi0Strip Strip energy " << candidate_it->second.energy() <<" is below threshold " ;
00151   
00152   // crystal should not be included...
00153   if ( (candidate_tmp->first == lastID )                    ||       // ...if it corresponds to a hit
00154        (candidate_it->second.energy() <= preshStripEnergyCut_ ) )   // ...if it has a negative or zero energy
00155     {
00156       return false;
00157     }
00158 
00159   return true;
00160 }
00161 
00162 // find strips in the road of size +/- preshSeededNstr_ from the central strip
00163 void EndcapPiZeroDiscriminatorAlgo::findPi0Road(ESDetId strip, EcalPreshowerNavigator& theESNav,
00164                                                                 int plane, vector<ESDetId>& vout) {
00165   if ( strip == ESDetId(0) ) return;
00166    ESDetId next;
00167    theESNav.setHome(strip);
00168    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: starts from strip " << strip ;
00169    
00170    if (plane == 1) {
00171      // east road
00172      int n_east= 0;
00173      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the East " ;
00174    
00175      while ( ((next=theESNav.east()) != ESDetId(0) && next != strip) ) {
00176         LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: East: " << n_east << " current strip is " << next ;
00177         
00178         vout.push_back(next);
00179         ++n_east;
00180         if (n_east == preshSeededNstr_) break;
00181      }
00182      // west road
00183      int n_west= 0;
00184      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the West " ;
00185      
00186      theESNav.home();
00187      while ( ((next=theESNav.west()) != ESDetId(0) && next != strip )) {
00188         LogTrace("EcalClusters") << "findPi0Road: West: " << n_west << " current strip is " << next ;
00189         
00190         vout.push_back(next);
00191         ++n_west;
00192         if (n_west == preshSeededNstr_) break;
00193      }
00194      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 1-st plane is " << n_east+n_west ;
00195      
00196   }
00197   else if (plane == 2) {
00198     // north road
00199     int n_north= 0;
00200     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the North " ;
00201     
00202     while ( ((next=theESNav.north()) != ESDetId(0) && next != strip) ) {
00203        LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: North: " << n_north << " current strip is " << next;
00204        
00205        vout.push_back(next);
00206        ++n_north;
00207        if (n_north == preshSeededNstr_) break;
00208     }
00209     // south road
00210     int n_south= 0;
00211     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the South " ;
00212     
00213     theESNav.home();
00214     while ( ((next=theESNav.south()) != ESDetId(0) && next != strip) ) {
00215        LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: South: " << n_south << " current strip is " << next ;
00216        
00217        vout.push_back(next);
00218        ++n_south;
00219        if (n_south == preshSeededNstr_) break;
00220     }
00221     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 2-nd plane is " << n_south+n_north ;
00222     
00223   }
00224   else {
00225     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: findPi0Road: Wrong plane number, null cluster will be returned! " ;
00226     
00227   } // end of if
00228 
00229   theESNav.home();
00230 }
00231 
00232 
00233 //===================================================================
00234 // EndcapPiZeroDiscriminatorAlgo::readWeightFile(...), a method that reads the weigths of the NN
00235 // INPUT: Weights_file
00236 // OUTPUT: I_H_Weight, H_Thresh, H_O_Weight, O_Thresh arrays
00237 //===================================================================
00238 void EndcapPiZeroDiscriminatorAlgo::readWeightFile(const char *Weights_file){
00239    FILE *weights;
00240 
00241    char *line;
00242    line = new char[80];
00243 
00244    bool checkinit=false;
00245 // Open the weights file, generated by jetnet, and read
00246 // in the nodes and weights
00247 //*******************************************************
00248   weights = fopen(Weights_file, "r");
00249   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: I opeded the Weights file  = " << Weights_file ;
00250   
00251 
00252   while( !feof(weights) ){
00253         fscanf(weights, "%s", line);
00254         if (line[0] == 'A') { //Read in ANN nodes: Layers, input , Hidden, Output
00255            fscanf(weights, "%d", &Layers);       // # of NN Layers  used
00256            fscanf(weights, "%d", &Indim);        // # of Inputs actually used
00257            fscanf(weights, "%d", &Hidden);       // # of hidden nodes
00258            fscanf(weights, "%d", &Outdim);   // # of output nodes
00259 
00260            inp_var = Indim + 1;
00261 
00262            I_H_Weight = new float[Indim*Hidden];
00263            H_Thresh = new float[Hidden];
00264            H_O_Weight = new float[Hidden*Outdim];
00265            O_Thresh = new float[Outdim];
00266            checkinit=true;
00267         }else if (line[0] == 'B') { // read in weights between hidden and intput nodes
00268             assert(checkinit);
00269             for (int i = 0; i<Indim; i++){
00270                 for (int j = 0; j<Hidden; j++){
00271                     fscanf(weights, "%f", &I_H_Weight[i*Hidden+j]);
00272                 }
00273             }
00274         }else if (line[0] == 'C'){       // Read in the thresholds for hidden nodes
00275             assert(checkinit);
00276             for (int i = 0; i<Hidden; i++){
00277                 fscanf(weights, "%f", &H_Thresh[i]);
00278             }
00279         }else if (line[0] == 'D'){ // read in weights between hidden and output nodes
00280             assert(checkinit);
00281             for (int i = 0; i<Hidden*Outdim; i++){
00282                 fscanf(weights, "%f", &H_O_Weight[i]);
00283             }
00284         }else if (line[0] == 'E'){ // read in the threshold for the output nodes
00285             assert(checkinit);
00286             for (int i = 0; i<Outdim; i++){
00287                 fscanf(weights, "%f", &O_Thresh[i]);
00288 
00289             }
00290         }
00291         else{edm::LogError("EEPi0Discrim")<< "EndcapPiZeroDiscriminatorAlgo: Not a Net file of Corrupted Net file " << endl;
00292         }
00293    }
00294    fclose(weights);
00295 }
00296 
00297 //=====================================================================================
00298 // EndcapPiZeroDiscriminatorAlgo::getNNoutput(int sel_wfile), a method that calculated the NN output
00299 // INPUT: sel_wfile -> Weight file selection
00300 // OUTPUT : nnout -> the NN output
00301 //=====================================================================================
00302 
00303 float EndcapPiZeroDiscriminatorAlgo::getNNoutput(int sel_wfile)
00304 {
00305  float* I_SUM;
00306  float* OUT;
00307  float nnout=0.0;
00308  int mij;
00309 
00310  I_SUM = new float[Hidden];
00311  OUT = new  float[Outdim];
00312 
00313  for(int k=0;k<Hidden;k++) I_SUM[k]=0.0;
00314  for(int k1=0;k1<Outdim;k1++) OUT[k1]=0.0;
00315 
00316  for (int h = 0; h<Hidden; h++){
00317      mij = h - Hidden;
00318      for (int i = 0; i<Indim; i++){
00319          mij = mij + Hidden;
00320          I_SUM[h] += I_H_Weight_all[mij+sel_wfile*Indim*Hidden + barrelstart*Nfiles_EE*EE_Indim*EE_Hidden] * input_var[i];
00321      }
00322      I_SUM[h] += H_Thresh_all[h+sel_wfile*Hidden + barrelstart*Nfiles_EE*EE_Hidden];
00323      for (int o1 = 0; o1<Outdim; o1++) {
00324         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]); 
00325 
00326      }
00327  }
00328  for (int o2 = 0; o2<Outdim; o2++){      
00329                 OUT[o2] += O_Thresh_all[barrelstart*Nfiles_EE*EE_Outdim + o2 + sel_wfile*Outdim]; 
00330   }
00331   nnout = Activation_fun(OUT[0]);
00332   LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: getNNoutput :: -> NNout = " <<  nnout ;
00333   
00334 
00335   delete I_SUM;
00336   delete OUT;
00337   delete input_var;
00338    
00339   return (nnout);
00340 }
00341 
00342 
00343 float EndcapPiZeroDiscriminatorAlgo::Activation_fun(float SUM){
00344         return( 1.0 / ( 1.0 + exp(-2.0*SUM) ) );
00345 }
00346 //=====================================================================================
00347 // EndcapPiZeroDiscriminatorAlgo::calculateNNInputVariables(...), a method that calculates the 25 input variables
00348 // INPUTS:
00349 // vph1 -> vector of the stip energies in 1st Preshower plane
00350 // vph2 -> vector of the stip energies in 2nd Preshower plane
00351 // pS1_max -> E1
00352 // pS9_max -> E9
00353 // pS25_max -> E25
00354 // OUTPUT: 
00355 // input_var[25] -> the 25 input to the NN variables array
00356 //=====================================================================================
00357 bool EndcapPiZeroDiscriminatorAlgo::calculateNNInputVariables(vector<float>& vph1, vector<float>& vph2, 
00358                                           float pS1_max, float pS9_max, float pS25_max, int EScorr)
00359 {
00360    input_var = new float[EE_Indim];
00361 
00362    bool valid_NNinput = true;
00363 
00364    /*   
00365    for(int i = 0; i<11;i++) {
00366    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Energies of the Preshower Strips in X plane = " << vph1[i] ;
00367    }
00368    
00369    for(int i = 0; i<11;i++) {
00370    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Energies of the Preshower Strips in Y plane = " << vph2[i] ;
00371    } 
00372    */
00373 
00374    // check if all Preshower info is availabla - If NOT use remaning info
00375    for(int k = 0; k<11; k++) {
00376       if(vph1[k] < 0 ) {
00377       
00378       LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " << k << " of X plane Do not exists" ;
00379      
00380          vph1[k] = 0.0;
00381       }          
00382       if(vph2[k] < 0 ) { 
00383         LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " << k << " of Y plane Do not exists" ;   
00384         
00385         vph2[k] = 0.0;
00386       }
00387    }   
00388 
00389    /*
00390    for(int i = 0; i<11;i++) {  
00391      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: After: Energies of the Preshower Strips in X plane = " << vph1[i] ;
00392    }
00393 
00394    for(int i = 0; i<11;i++) {
00395  
00396      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: After: Energies of the Preshower Strips in Y plane = " << vph2[i] ;
00397    }
00398    */
00399 
00400 
00401 // FIRST : Produce the 22 NN variables related with the Preshower 
00402 // --------------------------------------------------------------
00403 // New normalization of the preshower strip energies Aris 8/11/2004
00404    for(int kk=0;kk<11;kk++){
00405      input_var[kk] = fabs(vph1[kk]/0.01);
00406      input_var[kk + 11] = fabs(vph2[kk]/0.02);       
00407      if(input_var[kk] < 0.0001) input_var[kk] = 0.;
00408      if(input_var[kk + 11] < 0.0001) input_var[kk + 11] = 0.;
00409    }
00410    input_var[0] = fabs(input_var[0]/2.); 
00411    input_var[1] = fabs(input_var[1]/2.); 
00412    input_var[6] = fabs(input_var[6]/2.); 
00413    input_var[11] = fabs(input_var[11]/2.); 
00414    input_var[12] = fabs(input_var[12]/2.); 
00415    input_var[17] = fabs(input_var[17]/2.); 
00416    
00417 // correction for version > CMSSW_3_1_0_pre5 where extra enegry is given to the ES strips 
00418 // Aris 18/5/2009   
00419    if( EScorr == 1) { 
00420      input_var[0] -= 0.05;
00421      input_var[1] -= 0.035;
00422      input_var[2] -= 0.035;
00423      input_var[3] -= 0.02;
00424      input_var[4] -= 0.015;
00425      input_var[5] -= 0.0075;
00426      input_var[6] -= 0.035;
00427      input_var[7] -= 0.035;
00428      input_var[8] -= 0.02;
00429      input_var[9] -= 0.015;
00430      input_var[10] -= 0.0075;
00431     
00432      input_var[11] -= 0.05;
00433      input_var[12] -= 0.035;
00434      input_var[13] -= 0.035;
00435      input_var[14] -= 0.02;
00436      input_var[15] -= 0.015;
00437      input_var[16] -= 0.0075;
00438      input_var[17] -= 0.035;
00439      input_var[18] -= 0.035;
00440      input_var[19] -= 0.02;
00441      input_var[20] -= 0.015;
00442      input_var[21] -= 0.0075;   
00443      
00444      for(int kk1=0;kk1<22;kk1++){
00445        if(input_var[kk1] < 0 ) input_var[kk1] = 0.0;   
00446      }
00447    }
00448 // SECOND: Take the final NN variable related to the ECAL
00449 // -----------------------------------------------
00450    float ECAL_norm_factor = 500.;
00451    if(pS25_max>500&&pS25_max<=1000) ECAL_norm_factor = 1000;
00452    if(pS25_max>1000) ECAL_norm_factor = 7000;
00453 
00454    input_var[22] = pS1_max/ECAL_norm_factor;
00455    input_var[23] = pS9_max/ECAL_norm_factor;
00456    input_var[24] = pS25_max/ECAL_norm_factor;
00457 
00458    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S1/ECAL_norm_factor = " << input_var[22];
00459    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S9/ECAL_norm_factor = " << input_var[23];
00460    LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: S25/ECAL_norm_factor = " << input_var[24] ;
00461 
00462    for(int i=0;i<EE_Indim;i++){
00463      if(input_var[i] > 1.0e+00) {
00464        valid_NNinput = false;
00465        break;
00466      }
00467    }
00468    
00469    LogTrace("EcalClusters") << " valid_NNinput = " << valid_NNinput ;
00470    
00471 
00472    return valid_NNinput;
00473 }
00474 
00475 
00476 //=====================================================================================
00477 // EndcapPiZeroDiscriminatorAlgo::calculateBarrelNNInputVariables(...), a method that calculates
00478 // the 12 barrel NN input
00479 // OUTPUT:
00480 // input_var[12] -> the 12 input to the barrel NN variables array
00481 //=====================================================================================
00482 
00483 void EndcapPiZeroDiscriminatorAlgo::calculateBarrelNNInputVariables(float et, double s1, double s9, double s25,
00484                                                                     double m2, double cee, double cep,double cpp,
00485                                                                     double s4, double s6, double ratio,
00486                                                                     double xcog, double ycog)
00487 {
00488   input_var = new float[EB_Indim];
00489 
00490   double lam, lam1, lam2;
00491 
00492   if(xcog < 0.) { 
00493     input_var[0] = -xcog/s25;
00494   } else { 
00495     input_var[0] = xcog/s25;
00496   }    
00497 
00498   input_var[1] = cee/0.0004;
00499 
00500   if(cpp<.001) {
00501     input_var[2] = cpp/.001;
00502   } else  { 
00503     input_var[2] = 0.;
00504   }
00505 
00506   if(s9!=0.) {
00507     input_var[3] = s1/s9; 
00508     input_var[8] = s6/s9; 
00509     input_var[10] = (m2+s1)/s9;
00510   }
00511   else {
00512     input_var[3] = 0.;
00513     input_var[8] = 0.;
00514     input_var[10] = 0.;
00515   }
00516 
00517   if(s25-s1>0.) {
00518     input_var[4] = (s9-s1)/(s25-s1); 
00519   } else {
00520     input_var[4] = 0.;
00521   }  
00522 
00523   if(s25>0.) {
00524     input_var[5] = s4/s25; 
00525   }  else {
00526     input_var[5] = 0.;
00527   } 
00528 
00529   if(ycog < 0.) {
00530     input_var[6] = -ycog/s25; 
00531   } else {
00532     input_var[6] = ycog/s25;
00533   } 
00534 
00535   input_var[7] = ratio;
00536 
00537   lam=sqrt((cee -cpp)*(cee -cpp)+4*cep*cep);
00538   lam1=(cee + cpp + lam)/2;
00539   lam2=(cee + cpp - lam)/2;
00540 
00541   if(lam1 == 0) {
00542     input_var[9] = .0;
00543   } else {
00544     input_var[9] = lam2/lam1;
00545   }
00546   if(s4!=0.) {
00547     input_var[11] = (m2+s1)/s4;
00548   } else {
00549     input_var[11] = 0.;
00550   }
00551 
00552 }
00553 
00554 
00555 //=====================================================================================
00556 // EndcapPiZeroDiscriminatorAlgo::GetNNOutput(...), a method that calculates the NNoutput
00557 // INPUTS: Super Cluster Energy
00558 // OUTPUTS : NNoutput
00559 //=====================================================================================
00560 float EndcapPiZeroDiscriminatorAlgo::GetNNOutput(float EE_Et)
00561 {
00562     Layers = EE_Layers; Indim = EE_Indim; Hidden = EE_Hidden; Outdim = EE_Outdim;  barrelstart = 0;  
00563 
00564     float nnout = -1;
00565 // Print the  NN input variables that are related to the Preshower + ECAL
00566 // ------------------------------------------------------------------------
00567     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo::GetNNoutput :nn_invar_presh = " ;
00568     
00569      for(int k1=0;k1<Indim;k1++) {
00570      LogTrace("EcalClusters") << input_var[k1] << " " ;
00571     
00572      }
00573      LogTrace("EcalClusters") << " " ;
00574     
00575      // select the appropriate Weigth file
00576      int sel_wfile;
00577      if(EE_Et<25.0)                     {sel_wfile = 0;}
00578      else if(EE_Et>=25.0 && EE_Et<35.0) {sel_wfile = 1;}
00579      else if(EE_Et>=35.0 && EE_Et<45.0) {sel_wfile = 2;}
00580      else if(EE_Et>=45.0 && EE_Et<55.0) {sel_wfile = 3;}
00581      else                               {sel_wfile = 4;} 
00582 
00583      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: Et_SC = " << EE_Et << " and I select Weight file Number = " << sel_wfile ;
00584 
00585      
00586      nnout = getNNoutput(sel_wfile); // calculate the nnoutput for the given ECAL object
00587      
00588      LogTrace("EcalClusters") <<  "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
00589      
00590     
00591    return nnout;
00592 }
00593 
00594 
00595 //=====================================================================================
00596 // EndcapPiZeroDiscriminatorAlgo::GetBarrelNNOutput(...), a method that calculates the barrel NNoutput
00597 // INPUTS: Super Cluster Energy
00598 // OUTPUTS : NNoutput
00599 //=====================================================================================
00600 float EndcapPiZeroDiscriminatorAlgo::GetBarrelNNOutput(float EB_Et)
00601 {
00602 
00603     Layers = EB_Layers; Indim = EB_Indim; Hidden = EB_Hidden; Outdim = EB_Outdim;  barrelstart = 1;  
00604 
00605     float nnout = -1;
00606 // Print the  NN input variables that are related to the ECAL Barrel
00607 // ------------------------------------------------------------------------
00608     LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo::GetBarrelNNoutput :nn_invar_presh = " ;
00609 
00610     for(int k1=0;k1<Indim;k1++) {
00611       LogTrace("EcalClusters") << input_var[k1] << " " ;
00612       
00613     }
00614      LogTrace("EcalClusters") << " " ;
00615      
00616      // select the appropriate Weigth file
00617      int sel_wfile;
00618      if(EB_Et<25.0)                     {sel_wfile = 0;}
00619      else if(EB_Et>=25.0 && EB_Et<35.0) {sel_wfile = 1;}
00620      else if(EB_Et>=35.0 && EB_Et<45.0) {sel_wfile = 2;}
00621      else if(EB_Et>=45.0 && EB_Et<55.0) {sel_wfile = 3;}
00622      else                               {sel_wfile = 4;}
00623      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: E_SC = " << EB_Et << " and I select Weight file Number = " << sel_wfile ;
00624     
00625      nnout = getNNoutput(sel_wfile); // calculate the nnoutput for the given ECAL object
00626      
00627      LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
00628     
00629    return nnout;
00630 }
00631