CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
CorrectDeadChannelsNN.cc File Reference
#include "TMultiLayerPerceptron.h"
#include <TTree.h>
#include <TCanvas.h>
#include <TGraph2D.h>
#include <TSystem.h>
#include <TMath.h>
#include <TGraphErrors.h>
#include <TProfile2D.h>
#include <TProfile.h>
#include <TPostScript.h>
#include <TLegend.h>
#include <TStyle.h>
#include <TRandom.h>
#include "getopt.h"
#include <map>
#include <iostream>
#include "string.h"
#include <sstream>
#include <unistd.h>
#include <iomanip>
#include <fstream>
#include <algorithm>
#include <time.h>
#include "Test_Central_1500.cxx"
#include "Test_Side_1500.cxx"
#include "Test_Corner_1500.cxx"
#include "TestPos_100.cxx"

Go to the source code of this file.

Functions

double CorrectDeadChannelsNN (double *M5x5Input)
 

Function Documentation

double CorrectDeadChannelsNN ( double *  M5x5Input)

Definition at line 41 of file CorrectDeadChannelsNN.cc.

References epsilon, i, recoMuon::in, j, position, Test_Side_1500::value(), Test_Central_1500::value(), Test_Corner_1500::value(), and TestPos_100::value().

Referenced by EcalDeadChannelRecoveryAlgos::correct().

41  {
42 
43 
44  //From 07 May 07 we have as input a vector of :
45  //So now, we have a vector which is ordered around the Maximum containement and which contains a dead channel as:
46  //Filling of the vector : NxNaroundDC with N==11 Typo are possible...
47  // 060 is Maximum containement which is in +/- 2 from DC
48  //
49  // 120 119 118 117 116 115 114 113 112 111 110
50  // 109 108 107 106 105 104 103 102 101 100 099
51  // 098 097 096 095 094 093 092 091 090 089 088
52  // 087 086 085 084 083 082 081 080 079 078 077
53  // 076 075 074 073 072 071 070 069 068 067 066
54  // 065 064 063 062 061 060 059 058 057 056 055
55  // 054 053 052 051 050 049 048 047 046 045 044
56  // 043 042 041 040 039 038 037 036 035 034 033
57  // 032 031 030 029 028 027 026 025 024 023 022
58  // 021 020 019 018 017 016 015 014 013 012 011
59  // 010 009 008 007 006 005 004 003 002 001 000
61 
62  //Conversion between input and input need by NN:
63  double M5x5[25];
64  M5x5[0]=M5x5Input[60];
65  M5x5[1]=M5x5Input[71];
66  M5x5[2]=M5x5Input[49];
67  M5x5[3]=M5x5Input[61];
68  M5x5[4]=M5x5Input[72];
69  M5x5[5]=M5x5Input[50];
70  M5x5[6]=M5x5Input[59];
71  M5x5[7]=M5x5Input[70];
72  M5x5[8]=M5x5Input[48];
73  M5x5[9]=M5x5Input[82];
74  M5x5[10]=M5x5Input[83];
75  M5x5[11]=M5x5Input[81];
76  M5x5[12]=M5x5Input[38];
77  M5x5[13]=M5x5Input[39];
78  M5x5[14]=M5x5Input[37];
79  M5x5[15]=M5x5Input[58];
80  M5x5[16]=M5x5Input[69];
81  M5x5[17]=M5x5Input[80];
82  M5x5[18]=M5x5Input[47];
83  M5x5[19]=M5x5Input[36];
84  M5x5[20]=M5x5Input[62];
85  M5x5[21]=M5x5Input[73];
86  M5x5[22]=M5x5Input[84];
87  M5x5[23]=M5x5Input[51];
88  M5x5[24]=M5x5Input[40];
89 
90 
91  Test_Central_1500 *NNCentral =new Test_Central_1500();
92  Test_Side_1500 *NNAdj =new Test_Side_1500();
93  Test_Corner_1500 *NNCorner =new Test_Corner_1500();
94 
95  double in[6]; //Input of neurals networks
96  double NN1=-1;//Output Neural Net
97 
98  double epsilon = 0.0000001;
99 
100  //variables
101  bool Adjacent=false;
102  bool Corner=false;
103  bool Central=false;
104  double SUM8=-1;
105  double logX8=-1;
106  double logY8=-1;
107  double SUM24=-1;
108 
109  //First Find the position of the Dead Channel in 3x3 matrix
110  int IndDeadCha=-1;
111  for(int i =0 ; i< 9;i++){
112  if(fabs(M5x5[i])<epsilon && IndDeadCha >0){/*cout<<" Problem 2 dead channels in sum9! Can not correct "<<std::endl;*/ return 0.0;}
113  if(fabs(M5x5[i])<epsilon && IndDeadCha==-1)IndDeadCha=i;
114  }
115 
116 
117 
118 // std::cout<<" THE DEAD CHANNEL IS : " << IndDeadCha <<std::endl;
119 
120 // std::cout<<"========================="<<std::endl;
121 // std::cout<<M5x5[22]<<" "<<M5x5[10]<<" "<<M5x5[9] <<" "<<M5x5[11]<<" "<<M5x5[17]<<std::endl;
122 // std::cout<<M5x5[21]<<" "<<M5x5[4] <<" "<<M5x5[1] <<" "<<M5x5[7] <<" "<<M5x5[16]<<std::endl;
123 // std::cout<<M5x5[20]<<" "<<M5x5[3] <<" "<<M5x5[0] <<" "<<M5x5[6] <<" "<<M5x5[15]<<std::endl;
124 // std::cout<<M5x5[23]<<" "<<M5x5[5] <<" "<<M5x5[2] <<" "<<M5x5[8] <<" "<<M5x5[18]<<std::endl;
125 // std::cout<<M5x5[24]<<" "<<M5x5[13]<<" "<<M5x5[12]<<" "<<M5x5[14]<<" "<<M5x5[19]<<std::endl;
126 // std::cout<<"========================="<<std::endl;
127 
128 
129 // int isOK; isOK=0;
130 // if(IndDeadCha == 0)isOK=1;
131 // if(isOK == 0 )return 1000.0;
132 
133 
134 
135  int lineX1[3];
136  int lineX2[3];
137  int lineY1[3];
138  int lineY2[3];
139  int voisin1,voisin2,voisin3,voisin4;
140 
141  switch (IndDeadCha){
142  case 0:
143  lineX1[0]=3;lineX1[1]=4;lineX1[2]=5;
144  lineX2[0]=6;lineX2[1]=7;lineX2[2]=8;
145  lineY1[0]=1;lineY1[1]=4;lineY1[2]=7;
146  lineY2[0]=2;lineY2[1]=8;lineY2[2]=5;
147  voisin1=1;voisin2=2;voisin3=3;voisin4=6;
148  break;
149  case 1:
150  lineY1[0]=3;lineY1[1]=4;lineY1[2]=5;
151  lineY2[0]=6;lineY2[1]=7;lineY2[2]=8;
152  lineX2[0]=3;lineX2[1]=0;lineX2[2]=6;
153  lineX1[0]=2;lineX1[1]=8;lineX1[2]=5;
154  voisin1=0;voisin2=4;voisin3=7;voisin4=9;
155  break;
156  case 2:
157  lineY2[0]=3;lineY2[1]=4;lineY2[2]=5;
158  lineY1[0]=6;lineY1[1]=7;lineY1[2]=8;
159  lineX1[0]=1;lineX1[1]=4;lineX1[2]=7;
160  lineX2[0]=3;lineX2[1]=0;lineX2[2]=6;
161  voisin1=0;voisin2=5;voisin3=8;voisin4=12;
162  break;
163  case 3:
164  lineX2[0]=0;lineX2[1]=1;lineX2[2]=2;
165  lineX1[0]=6;lineX1[1]=7;lineX1[2]=8;
166  lineY1[0]=1;lineY1[1]=4;lineY1[2]=7;
167  lineY2[0]=2;lineY2[1]=8;lineY2[2]=5;
168  voisin1=0;voisin2=20;voisin3=4;voisin4=5;
169  break;
170  case 4:
171  lineX1[0]=0;lineX1[1]=1;lineX1[2]=2;
172  lineX2[0]=6;lineX2[1]=7;lineX2[2]=8;
173  lineY1[0]=0;lineY1[1]=3;lineY1[2]=6;
174  lineY2[0]=2;lineY2[1]=8;lineY2[2]=5;
175  voisin1=1;voisin2=10;voisin3=3;voisin4=21;
176  break;
177  case 5:
178  lineY1[0]=0;lineY1[1]=1;lineY1[2]=2;
179  lineY2[0]=6;lineY2[1]=7;lineY2[2]=8;
180  lineX2[0]=4;lineX2[1]=1;lineX2[2]=7;
181  lineX1[0]=0;lineX1[1]=3;lineX1[2]=6;
182  voisin1=13;voisin2=2;voisin3=3;voisin4=23;
183  break;
184  case 6:
185  lineX1[0]=3;lineX1[1]=4;lineX1[2]=5;
186  lineX2[0]=0;lineX2[1]=1;lineX2[2]=2;
187  lineY1[0]=1;lineY1[1]=4;lineY1[2]=7;
188  lineY2[0]=5;lineY2[1]=2;lineY2[2]=8;
189  voisin1=15;voisin2=0;voisin3=7;voisin4=8;
190  break;
191  case 7:
192  lineY2[0]=3;lineY2[1]=4;lineY2[2]=5;
193  lineY1[0]=0;lineY1[1]=1;lineY1[2]=2;
194  lineX1[0]=0;lineX1[1]=3;lineX1[2]=6;
195  lineX2[0]=2;lineX2[1]=8;lineX2[2]=5;
196  voisin1=11;voisin2=1;voisin3=6;voisin4=16;
197  break;
198  case 8:
199  lineX2[0]=3;lineX2[1]=4;lineX2[2]=5;
200  lineX1[0]=0;lineX1[1]=1;lineX1[2]=2;
201  lineY2[0]=1;lineY2[1]=4;lineY2[2]=7;
202  lineY1[0]=0;lineY1[1]=3;lineY1[2]=6;
203  voisin1=6;voisin2=2;voisin3=18;voisin4=14;
204  break;
205  default:
206  // std::cout<<" Error, Dead Channel to far from main containement one, no correction is applied"<<std::endl;
207  return 0.0;
208  break;
209  }//end switch
210 
211  float XL8=-50;
212  float XR8=-50;
213  float YL8=-50;
214  float YR8=-50;
215 
216  XL8=0;
217  XR8=0;
218  for(int j=0;j<3;j++){XL8+=M5x5[lineX1[j]];}
219  for(int j=0;j<3;j++){XR8+=M5x5[lineX2[j]];}
220  YL8=0;
221  YR8=0;
222  for(int j=0;j<3;j++)YL8+=M5x5[lineY1[j]];
223  for(int j=0;j<3;j++)YR8+=M5x5[lineY2[j]];
224 
225  SUM8 =0;
226  for(int j=0;j<9;j++)if(j!=IndDeadCha)SUM8+=M5x5[j];
227 
228  float XR24=XR8+M5x5[11]+M5x5[14]+M5x5[15]+M5x5[16]+M5x5[17]+M5x5[18]+M5x5[19];
229  float YR24=YR8+M5x5[18]+M5x5[23]+M5x5[24]+M5x5[13]+M5x5[12]+M5x5[14]+M5x5[19];
230 
231  float sum24 = 0.;
232  for(int j=0;j<25;j++)if(j!=IndDeadCha)sum24+=M5x5[j];
233 
234  if(XR8 > 0 && XL8>0 && YL8 >0 && YR8>0 && SUM8>0 && XR24>0 &&YR24>0){
235  logX8=TMath::Log(XL8/XR8);
236  logY8=TMath::Log(YL8/YR8);
237  SUM24=sum24;
238  //Added 15 Janvier 2007 to get ride of energy dependance!
239  SUM24 = SUM8/SUM24;
240 
241  Adjacent=false;
242  Corner=false;
243  Central=false;
244 
245  //Il faut trouver le canal d'energie maximal, et par rapport a ce canal regarde l'energie dans les adjacent pour voir si nous sommes sur le central.
246  float maxi;
247  int IndMax;
248  if(IndDeadCha!=0){
249  maxi=M5x5[0];
250  IndMax=0;
251  }else{
252  maxi=M5x5[1];
253  IndMax=1;
254  }
255  for(int j=IndMax;j<9;j++){
256  if(j!=IndDeadCha){
257  if(M5x5[j] > maxi){IndMax=j;maxi=M5x5[j];}
258  }
259  }
260  float Secmaxi=0;
261  int IndSecMax=0;
262  for(int j=IndSecMax;j<9;j++){
263  if(j!=IndDeadCha && j!=IndMax){
264  if(M5x5[j] > Secmaxi){IndSecMax=j;Secmaxi=M5x5[j];}
265  }
266  }
267 
268  in[0]=logX8;
269  in[1]=logY8;
270  in[2]=SUM24;
271 
272  //Define Adjacent/Central/Side with a Neural Net
274  float indpos = position->value(0,in[0],in[1],in[2],(SUM8-maxi)/sum24,(maxi-Secmaxi)/maxi,M5x5[voisin1]/SUM8,M5x5[voisin2]/SUM8,M5x5[voisin3]/SUM8,M5x5[voisin4]/SUM8);
275  if( indpos < 1.5){
276  Central=true;
277  Adjacent=false;
278  Corner=false;
279  }else{
280  if(indpos <2.5){
281  Central=false;
282  Adjacent=true;
283  Corner=false;
284  }else{
285  Central=false;
286  Adjacent=false;
287  Corner=true;
288  }
289  }
290  delete position;
291 
292 
293  if(logX8!=-50 && logY8!=-50 && SUM8>0 && XR24>0 &&YR24>0){
294  if(Adjacent)NN1 =NNAdj->value(0,in[0],in[1],in[2]);
295  if(Central)NN1 =NNCentral->value(0,in[0],in[1],in[2]);
296  if(Corner)NN1 =NNCorner->value(0,in[0],in[1],in[2]);
297  }else{
298  NN1=0;
299  }
300  }
301 
302  return NN1*SUM8;
303 }
int i
Definition: DBlmapReader.cc:9
double value(int index, double in0, double in1, double in2)
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
double value(int index, double in0, double in1, double in2)
double value(int index, double in0, double in1, double in2)
int j
Definition: DBlmapReader.cc:9
const double epsilon
double value(int index, double in0, double in1, double in2, double in3, double in4, double in5, double in6, double in7, double in8)