14 preshStripEnergyCut_(stripEnergyCut), preshSeededNstr_(nStripCut), pathToFiles_(path)
20 string file_pt[5] = {
"20",
"30",
"40",
"50",
"60"};
21 string file_barrel_pt[5] = {
"20",
"30",
"40",
"50",
"60"};
23 string nn_paterns_file =
"";
24 for(
int j=0;
j<Nfiles_EE;
j++) {
25 nn_paterns_file =
"endcapPiZeroDiscriminatorWeights_et"+file_pt[
j]+
".wts";
27 readWeightFile(WFile.fullPath().c_str());
34 for(
int i=0;
i<Indim*Hidden;
i++) I_H_Weight_all.push_back(I_H_Weight[
i]);
35 for(
int i=0;i<Hidden;i++) H_Thresh_all.push_back(H_Thresh[i]);
36 for(
int i=0;i<Outdim*Hidden;i++) H_O_Weight_all.push_back(H_O_Weight[i]);
37 for(
int i=0;i<Outdim;i++) O_Thresh_all.push_back(O_Thresh[i]);
40 for(
int k=0;
k<Nfiles_EB;
k++) {
41 nn_paterns_file =
"barrelPiZeroDiscriminatorWeights_et"+file_barrel_pt[
k]+
".wts";
43 readWeightFile(WFile.fullPath().c_str());
50 for(
int i=0;i<Indim*Hidden;i++) I_H_Weight_all.push_back(I_H_Weight[i]);
51 for(
int i=0;i<Hidden;i++) H_Thresh_all.push_back(H_Thresh[i]);
52 for(
int i=0;i<Outdim*Hidden;i++) H_O_Weight_all.push_back(H_O_Weight[i]);
53 for(
int i=0;i<Outdim;i++) O_Thresh_all.push_back(O_Thresh[i]);
65 vector<float> vout_stripE;
68 if ( rechits_map->size() == 0 ) {
69 edm::LogWarning(
"EndcapPiZeroDiscriminatorAlgo") <<
"RecHitsMap has size 0.";
75 vector<ESDetId> road_2d;
78 int plane = strip.
plane();
80 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPreshVectors: Preshower Seeded Algorithm - looking for clusters" <<
"n"<<
"findPreshVectors: Preshower is intersected at strip " << strip.
strip() <<
", at plane " << plane ;
84 for(
int i=0;i<11;i++) {
85 vout_stripE.push_back(-100.);
90 road_2d.push_back(strip);
94 navigator.setHome(strip);
98 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo:findPreshVectors: Total number of strips in the central road: " << road_2d.size() ;
101 RecHitsMap::iterator final_strip = rechits_map->end();
105 ESDetId last_stripID = final_strip->first;
108 vector<ESDetId>::iterator itID;
109 for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
110 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPreshVectors: ID = " << *itID ;
112 RecHitsMap::iterator strip_it = rechits_map->find(*itID);
114 E = strip_it->second.energy();
116 vout_stripE.push_back(E);
117 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPreshVectors: E = " << E ;
124 vector<float> vout_ElevenStrips_Energy;
125 for(
int i=0;i<11;i++)
127 vout_ElevenStrips_Energy.push_back(0.);
130 for(
unsigned int i=0;i<vout_stripE.size();i++)
132 vout_ElevenStrips_Energy[
i] = vout_stripE.at(i);
136 return vout_ElevenStrips_Energy;
144 RecHitsMap::iterator candidate_tmp = candidate_it;
147 if (candidate_tmp->first == lastID )
148 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: goodPi0Strip No such a strip in rechits_map " ;
150 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: goodPi0Strip Strip energy " << candidate_it->second.energy() <<
" is below threshold " ;
153 if ( (candidate_tmp->first == lastID ) ||
164 int plane, vector<ESDetId>& vout) {
165 if ( strip ==
ESDetId(0) )
return;
167 theESNav.setHome(strip);
168 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: starts from strip " << strip ;
173 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the East " ;
175 while ( ((next=theESNav.east()) !=
ESDetId(0) && next != strip) ) {
176 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: East: " << n_east <<
" current strip is " <<
next ;
178 vout.push_back(next);
184 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the West " ;
187 while ( ((next=theESNav.west()) !=
ESDetId(0) && next != strip )) {
188 LogTrace(
"EcalClusters") <<
"findPi0Road: West: " << n_west <<
" current strip is " <<
next ;
190 vout.push_back(next);
194 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 1-st plane is " << n_east+n_west ;
197 else if (plane == 2) {
200 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the North " ;
202 while ( ((next=theESNav.north()) !=
ESDetId(0) && next != strip) ) {
203 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: North: " << n_north <<
" current strip is " <<
next;
205 vout.push_back(next);
211 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: Go to the South " ;
214 while ( ((next=theESNav.south()) !=
ESDetId(0) && next != strip) ) {
215 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: South: " << n_south <<
" current strip is " <<
next ;
217 vout.push_back(next);
221 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: Total number of strips found in the road at 2-nd plane is " << n_south+n_north ;
225 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: findPi0Road: Wrong plane number, null cluster will be returned! " ;
244 bool checkinit=
false;
248 weights = fopen(Weights_file,
"r");
249 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: I opeded the Weights file = " << Weights_file ;
252 while( !feof(weights) ){
253 fscanf(weights,
"%s", line);
254 if (line[0] ==
'A') {
255 fscanf(weights,
"%d", &
Layers);
256 fscanf(weights,
"%d", &
Indim);
257 fscanf(weights,
"%d", &Hidden);
258 fscanf(weights,
"%d", &
Outdim);
267 }
else if (line[0] ==
'B') {
269 for (
int i = 0; i<
Indim; i++){
274 }
else if (line[0] ==
'C'){
276 for (
int i = 0; i<
Hidden; i++){
277 fscanf(weights,
"%f", &
H_Thresh[i]);
279 }
else if (line[0] ==
'D'){
281 for (
int i = 0; i<Hidden*
Outdim; i++){
284 }
else if (line[0] ==
'E'){
286 for (
int i = 0; i<
Outdim; i++){
287 fscanf(weights,
"%f", &
O_Thresh[i]);
291 else{
edm::LogError(
"EEPi0Discrim")<<
"EndcapPiZeroDiscriminatorAlgo: Not a Net file of Corrupted Net file " << endl;
310 I_SUM =
new float[
Hidden];
314 for(
int k1=0;k1<
Outdim;k1++) OUT[k1]=0.0;
318 for (
int i = 0; i<
Indim; i++){
323 for (
int o1 = 0; o1<
Outdim; o1++) {
328 for (
int o2 = 0; o2<
Outdim; o2++){
332 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: getNNoutput :: -> NNout = " << nnout ;
344 return( 1.0 / ( 1.0 +
exp(-2.0*SUM) ) );
358 float pS1_max,
float pS9_max,
float pS25_max,
int EScorr)
362 bool valid_NNinput =
true;
375 for(
int k = 0;
k<11;
k++) {
378 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " <<
k <<
" of X plane Do not exists" ;
383 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: Oops!!! Preshower Info for strip : " <<
k <<
" of Y plane Do not exists" ;
444 for(
int kk1=0;kk1<22;kk1++){
450 float ECAL_norm_factor = 500.;
451 if(pS25_max>500&&pS25_max<=1000) ECAL_norm_factor = 1000;
452 if(pS25_max>1000) ECAL_norm_factor = 7000;
454 input_var[22] = pS1_max/ECAL_norm_factor;
455 input_var[23] = pS9_max/ECAL_norm_factor;
456 input_var[24] = pS25_max/ECAL_norm_factor;
458 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: S1/ECAL_norm_factor = " <<
input_var[22];
459 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: S9/ECAL_norm_factor = " <<
input_var[23];
460 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: S25/ECAL_norm_factor = " <<
input_var[24] ;
464 valid_NNinput =
false;
469 LogTrace(
"EcalClusters") <<
" valid_NNinput = " << valid_NNinput ;
472 return valid_NNinput;
484 double m2,
double cee,
double cep,
double cpp,
485 double s4,
double s6,
double ratio,
486 double xcog,
double ycog)
490 double lam, lam1, lam2;
537 lam=
sqrt((cee -cpp)*(cee -cpp)+4*cep*cep);
538 lam1=(cee + cpp + lam)/2;
539 lam2=(cee + cpp - lam)/2;
567 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo::GetNNoutput :nn_invar_presh = " ;
569 for(
int k1=0;k1<
Indim;k1++) {
577 if(EE_Et<25.0) {sel_wfile = 0;}
578 else if(EE_Et>=25.0 && EE_Et<35.0) {sel_wfile = 1;}
579 else if(EE_Et>=35.0 && EE_Et<45.0) {sel_wfile = 2;}
580 else if(EE_Et>=45.0 && EE_Et<55.0) {sel_wfile = 3;}
581 else {sel_wfile = 4;}
583 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: Et_SC = " << EE_Et <<
" and I select Weight file Number = " << sel_wfile ;
588 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
608 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo::GetBarrelNNoutput :nn_invar_presh = " ;
610 for(
int k1=0;k1<
Indim;k1++) {
618 if(EB_Et<25.0) {sel_wfile = 0;}
619 else if(EB_Et>=25.0 && EB_Et<35.0) {sel_wfile = 1;}
620 else if(EB_Et>=35.0 && EB_Et<45.0) {sel_wfile = 2;}
621 else if(EB_Et>=45.0 && EB_Et<55.0) {sel_wfile = 3;}
622 else {sel_wfile = 4;}
623 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: E_SC = " << EB_Et <<
" and I select Weight file Number = " << sel_wfile ;
627 LogTrace(
"EcalClusters") <<
"EndcapPiZeroDiscriminatorAlgo: ===================> GetNNOutput : NNout = " << nnout ;
CaloNavigator< ESDetId > EcalPreshowerNavigator
void findPi0Road(ESDetId strip, EcalPreshowerNavigator &theESNav, int plane, std::vector< ESDetId > &vout)
float GetNNOutput(float EE_Et)
void 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)
float Activation_fun(float SUM)
std::vector< float > O_Thresh_all
std::map< DetId, EcalRecHit > RecHitsMap
EndcapPiZeroDiscriminatorAlgo()
double preshStripEnergyCut_
tuple path
else: Piece not in the list, fine.
std::vector< float > findPreshVector(ESDetId strip, RecHitsMap *rechits_map, CaloSubdetectorTopology *topology_p)
float getNNoutput(int sel_wfile)
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
std::vector< float > I_H_Weight_all
float GetBarrelNNOutput(float EB_Et)
bool goodPi0Strip(RecHitsMap::iterator candidate_it, ESDetId lastID)
std::vector< float > H_O_Weight_all
void readWeightFile(const char *WFile)
bool calculateNNInputVariables(std::vector< float > &vph1, std::vector< float > &vph2, float pS1_max, float pS9_max, float pS25_max, int EScorr)
std::vector< float > H_Thresh_all