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
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());
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());
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
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) ) {
00084 for(int i=0;i<11;i++) {
00085 vout_stripE.push_back(-100.);
00086 }
00087 }
00088
00089
00090 road_2d.push_back(strip);
00091
00092
00093 EcalPreshowerNavigator navigator(strip, topology_p);
00094 navigator.setHome(strip);
00095
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
00101 RecHitsMap::iterator final_strip = rechits_map->end();
00102
00103
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)) {
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
00122
00123
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
00136 return vout_ElevenStrips_Energy;
00137
00138
00139 }
00140
00141
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
00153 if ( (candidate_tmp->first == lastID ) ||
00154 (candidate_it->second.energy() <= preshStripEnergyCut_ ) )
00155 {
00156 return false;
00157 }
00158
00159 return true;
00160 }
00161
00162
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
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
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
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
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 }
00228
00229 theESNav.home();
00230 }
00231
00232
00233
00234
00235
00236
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
00246
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') {
00255 fscanf(weights, "%d", &Layers);
00256 fscanf(weights, "%d", &Indim);
00257 fscanf(weights, "%d", &Hidden);
00258 fscanf(weights, "%d", &Outdim);
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') {
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'){
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'){
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'){
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
00299
00300
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
00348
00349
00350
00351
00352
00353
00354
00355
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
00366
00367
00368
00369
00370
00371
00372
00373
00374
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
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
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
00418
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
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
00478
00479
00480
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
00557
00558
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
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
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);
00587
00588 LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
00589
00590
00591 return nnout;
00592 }
00593
00594
00595
00596
00597
00598
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
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
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);
00626
00627 LogTrace("EcalClusters") << "EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
00628
00629 return nnout;
00630 }
00631