CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
CorrectDeadChannelsClassic.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 "PositionCorrector.cxx"
#include "SplineCorrector.cxx"
#include "ExponCorrector.cxx"

Go to the source code of this file.

Functions

double CorrectDeadChannelsClassic (double *M11x11Input, const int DCeta)
 

Function Documentation

double CorrectDeadChannelsClassic ( double *  M11x11Input,
const int  DCeta 
)

Definition at line 40 of file CorrectDeadChannelsClassic.cc.

References PositionCorrector::CORRX(), PositionCorrector::CORRY(), gather_cfg::cout, epsilon, j, SplineCorrector::value(), and ExponCorrector::value().

Referenced by EcalDeadChannelRecoveryAlgos::correct().

40  {
41 
42  //std::cout<<"Inside Correction Function ...."<< std::endl;
43 
44  PositionCorrector *PosCorr = new PositionCorrector();
45  SplineCorrector *SplCorr = new SplineCorrector();
46  ExponCorrector *ExpCorr = new ExponCorrector();
47 
48 
49 
50 
51  double epsilon = 0.0000001;
52  float NEWx,NEWy;
53  float estimX,estimY,SUMlogFR;
54  float SUM24;
55  float crE[25];
56 
57 
58 
59  crE[0]=M11x11Input[40];
60  crE[1]=M11x11Input[51];
61  crE[2]=M11x11Input[62];
62  crE[3]=M11x11Input[73];
63  crE[4]=M11x11Input[84];
64 
65  crE[5]=M11x11Input[39];
66  crE[6]=M11x11Input[50];
67  crE[7]=M11x11Input[61];
68  crE[8]=M11x11Input[72];
69  crE[9]=M11x11Input[83];
70 
71  crE[10]=M11x11Input[38];
72  crE[11]=M11x11Input[49];
73  crE[12]=M11x11Input[60];
74  crE[13]=M11x11Input[71];
75  crE[14]=M11x11Input[82];
76 
77  crE[15]=M11x11Input[37];
78  crE[16]=M11x11Input[48];
79  crE[17]=M11x11Input[59];
80  crE[18]=M11x11Input[70];
81  crE[19]=M11x11Input[81];
82 
83  crE[20]=M11x11Input[36];
84  crE[21]=M11x11Input[47];
85  crE[22]=M11x11Input[58];
86  crE[23]=M11x11Input[69];
87  crE[24]=M11x11Input[80];
88 
89 
90  std::cout<<"Inside Correction Function : Dead Channel ETA" << DCeta << std::endl;
91  //revert crystal matrix because of negative eta
92  if(DCeta<0){
93 
94  crE[0]=M11x11Input[36];
95  crE[1]=M11x11Input[47];
96  crE[2]=M11x11Input[58];
97  crE[3]=M11x11Input[69];
98  crE[4]=M11x11Input[80];
99 
100  crE[5]=M11x11Input[37];
101  crE[6]=M11x11Input[48];
102  crE[7]=M11x11Input[59];
103  crE[8]=M11x11Input[70];
104  crE[9]=M11x11Input[81];
105 
106  crE[10]=M11x11Input[38];
107  crE[11]=M11x11Input[49];
108  crE[12]=M11x11Input[60];
109  crE[13]=M11x11Input[71];
110  crE[14]=M11x11Input[82];
111 
112  crE[15]=M11x11Input[39];
113  crE[16]=M11x11Input[50];
114  crE[17]=M11x11Input[61];
115  crE[18]=M11x11Input[72];
116  crE[19]=M11x11Input[83];
117 
118  crE[20]=M11x11Input[40];
119  crE[21]=M11x11Input[51];
120  crE[22]=M11x11Input[62];
121  crE[23]=M11x11Input[73];
122  crE[24]=M11x11Input[84];
123 
124 
125  }
126 
127 
128 
129 
130 
131  float SUMuu = 0.0;
132  float SUMu = 0.0;
133  float SUMdd = 0.0;
134  float SUMd = 0.0;
135  float SUMll = 0.0;
136  float SUMl = 0.0;
137  float SUMrr = 0.0;
138  float SUMr = 0.0;
139 
140  SUMuu = crE[4] + crE[9] + crE[14] + crE[19] + crE[24];
141  SUMu = crE[3] + crE[8] + crE[13] + crE[18] + crE[23];
142  SUMd = crE[1] + crE[6] + crE[11] + crE[16] + crE[21];
143  SUMdd = crE[0] + crE[5] + crE[10] + crE[15] + crE[20];
144 
145  SUMll = crE[0] + crE[1] + crE[2] + crE[3] + crE[4];
146  SUMl = crE[5] + crE[6] + crE[7] + crE[8] + crE[9];
147  SUMr = crE[15] + crE[16] + crE[17] + crE[18] + crE[19];
148  SUMrr = crE[20] + crE[21] + crE[22] + crE[23] + crE[24];
149 
150  /*
151  std::cout<<"================================================================="<<std::endl;
152  std::cout<<crE[4]<<setw(12)<<crE[9]<<setw(12)<<crE[14]<<setw(12)<<crE[19]<<setw(12)<<crE[24]<<std::endl;
153  std::cout<<crE[3]<<setw(12)<<crE[8]<<setw(12)<<crE[13]<<setw(12)<<crE[18]<<setw(12)<<crE[23]<<std::endl;
154  std::cout<<crE[2]<<setw(12)<<crE[7]<<setw(12)<<crE[12]<<setw(12)<<crE[17]<<setw(12)<<crE[22]<<std::endl;
155  std::cout<<crE[1]<<setw(12)<<crE[6]<<setw(12)<<crE[11]<<setw(12)<<crE[16]<<setw(12)<<crE[21]<<std::endl;
156  std::cout<<crE[0]<<setw(12)<<crE[5]<<setw(12)<<crE[10]<<setw(12)<<crE[15]<<setw(12)<<crE[20]<<std::endl;
157  std::cout<<"================================================================="<<std::endl;
158  */
159 
160 
162 
163  float XLOW[50],XHIG[50],YLOW[50],YHIG[50];
164  float CentX[50],CentY[50];
165  int NSUBS = 25;
166  for(int ix=0; ix<5 ; ix++){
167  for(int iy=0; iy<5 ; iy++){
168  int isub= ix*5+iy ;
169 
170  XLOW[isub]= -10.0 +float(ix)*4.0;
171  XHIG[isub]= XLOW[isub] + 4.0;
172  YLOW[isub]= -10.0 +float(iy)*4.0;;
173  YHIG[isub]= YLOW[isub] + 4.0;
174 
175  CentX[isub]=(XHIG[isub]+XLOW[isub])/2.0;
176  CentY[isub]=(YHIG[isub]+YLOW[isub])/2.0;
177  }
178  }
179 
180 
181 
182 
184 
185 
186 
187  //First Find the position of the Dead Channel in 3x3 matrix
188  int DeadCR = -1;
189  for(int ix=1; ix<4 ; ix++){
190  for(int iy=1; iy<4 ; iy++){
191  int idx = ix*5+iy;
192  if(fabs(crE[idx])<epsilon && DeadCR >0){std::cout<<" Problem 2 dead channels in sum9! Can not correct ... I return 0.0"<<std::endl; return 0.0;}
193  if(fabs(crE[idx])<epsilon && DeadCR==-1)DeadCR=idx;
194  }
195  }
196 
197 
198  //std::cout<<" THE DEAD CHANNEL IS : "<< DeadCR <<std::endl;
199  SUM24=0.0;
200  for(int j=0;j<25;j++)SUM24+=crE[j];
201  if(DeadCR == -1){std::cout<<" No Dead Channel in 3x3 ! I don't correct ... I return 0.0 ....look S25 = "<<SUM24<<std::endl; return 0.0;}
202 
203 
204 
205  // CR=12 must be the maximum in 5x5 !!!
206  int CisMax; CisMax=-1;
207  float MaxEin5x5; MaxEin5x5=0.0;
208  for(int ic=0;ic<24;ic++){
209  if(crE[ic]>MaxEin5x5){
210  MaxEin5x5 = crE[ic];
211  CisMax = ic ;
212  }
213  }
214  if(CisMax != 12){
215  std::cout<<"ERROR ----> Central has NOT the MAX energy in 5x5 ... I return 0.0"<<std::endl;
216  return 0.0;
217  }
218 
219 
220 
221  //AVOID BREM or OTHER DEPOSITIONS IN 5x5
222  //std::cout<<" CHECK FOR BREM or UNUSUAL PATTERNS ... " <<std::endl;
223  //std::cout<<" SUMuu="<<SUMuu<<",SUMu="<<SUMu<<",SUMd="<<SUMd <<",SUMdd="<<SUMdd <<",SUMll="<<SUMll <<",SUMl="<<SUMl <<",SUMrr="<<SUMrr <<",SUMr="<<SUMr <<std::endl;
224  if(DeadCR==6 && (SUMuu>SUMu || SUMrr>SUMr || SUMll>3.0*SUMr || SUMdd>3.0*SUMu)){std::cout<<"Unusual Pattern in 6 I return 0.0"<<std::endl; return 0.0;}
225  if(DeadCR==8 && (SUMdd>SUMd || SUMrr>SUMr || SUMll>3.0*SUMr || SUMuu>3.0*SUMd)){std::cout<<"Unusual Pattern in 8 I return 0.0"<<std::endl; return 0.0;}
226  if(DeadCR==16 && (SUMuu>SUMu || SUMll>SUMl || SUMrr>3.0*SUMl || SUMdd>3.0*SUMu)){std::cout<<"Unusual Pattern in 16 I return 0.0"<<std::endl; return 0.0;}
227  if(DeadCR==18 && (SUMdd>SUMd || SUMll>SUMl || SUMrr>3.0*SUMl || SUMuu>3.0*SUMd)){std::cout<<"Unusual Pattern in 18 I return 0.0"<<std::endl; return 0.0;}
228 
229  if(DeadCR==7 && (SUMuu>SUMu || SUMdd>SUMd || SUMrr>SUMr)){std::cout<<"Unusual Pattern in 7 I return 0.0"<<std::endl; return 0.0;}
230  if(DeadCR==17 && (SUMuu>SUMu || SUMdd>SUMd || SUMll>SUMl)){std::cout<<"Unusual Pattern in 17 I return 0.0"<<std::endl; return 0.0;}
231  if(DeadCR==11 && (SUMll>SUMl || SUMrr>SUMr || SUMuu>SUMu)){std::cout<<"Unusual Pattern in 11 I return 0.0"<<std::endl; return 0.0;}
232  if(DeadCR==13 && (SUMll>SUMl || SUMrr>SUMr || SUMdd>SUMd)){std::cout<<"Unusual Pattern in 13 I return 0.0"<<std::endl; return 0.0;}
233 
234 
235  //std::cout<<" NO BREM OR OTHER DEPOSITIONS IN THIS PATTERN --- I continue ... " <<std::endl;
236 
237  SUM24=0.0;
238  for(int j=0;j<25;j++)if(j!=DeadCR)SUM24+=crE[j];
239  // std::cout<<" THE Enorm is : "<< SUM24 <<std::endl;
240 
241 
242  //CHECK IF IT IS REALLY THE CORRECT DEAD CHANNEL
243  // This will be included in the next version
244 
245 
246 
249  // Estimate X,Y
252 
253  NEWx=0.0; NEWy=0.0;
254 
255 
256  estimX=0.0;
257  estimY=0.0;
258  SUMlogFR=0.0;
259  for(int ix=0; ix<5 ; ix++){
260  for(int iy=0; iy<5 ; iy++){
261  int idx = ix*5+iy;
262 
263  float xpos = 20.0*(float(ix)-2.0);
264  float ypos = 20.0*(float(iy)-2.0);
265 
266  if(idx != DeadCR){
267 
268  // weight definition
269  float w = crE[idx]/SUM24;
270 
271  if(w<0.0)w=0.0;
272  SUMlogFR = SUMlogFR + w;
273 
274  estimX = estimX + xpos*w;
275  estimY = estimY + ypos*w;
276  }
277 
278  } // iy loop
279  }// ix loop
280 
281  estimX = estimX/SUMlogFR;
282  estimY = estimY/SUMlogFR;
283 
284 
285  NEWx = PosCorr->CORRX(DeadCR,0,50,estimX);
286  NEWy = PosCorr->CORRY(DeadCR,0,50,estimY);
287 
288 
289  //std::cout<<" FINAL X,Y calculation: DC , estimX, estimY, NEWx, NEWy : "<<DeadCR<<" "<<estimX<<" "<<estimY<<" "<< NEWx<<" "<<NEWy<<std::endl;
290 
291 
292  if(DeadCR==7 && (estimX>7.0 || estimX< -2.0 || estimY>10.0 || estimY<-9.0) ){std::cout<<"DC=7 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
293  if(DeadCR==17 && (estimX>2.5 || estimX< -8.0 || estimY>10.0 || estimY<-8.0) ){std::cout<<"DC=17 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
294  if(DeadCR==11 && (estimX>8.0 || estimX<-10.0 || estimY> 9.0 || estimY<-4.0) ){std::cout<<"DC=11 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
295  if(DeadCR==13 && (estimX>8.0 || estimX<-10.0 || estimY> 4.5 || estimY<-8.0) ){std::cout<<"DC=13 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
296 
297  if(DeadCR==12 && (estimX>18.0 || estimX<-18.0 || estimY> 17.0 || estimY<-17.0) ){std::cout<<"DC=12 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
298 
299  if(DeadCR==6 && (estimX>8.0 || estimX< -9.0 || estimY>10.0 || estimY<-8.0) ){std::cout<<"DC=6 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
300  if(DeadCR==8 && (estimX>8.0 || estimX< -9.0 || estimY> 9.0 || estimY<-8.5) ){std::cout<<"DC=8 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
301  if(DeadCR==16 && (estimX>7.0 || estimX<-10.0 || estimY> 9.0 || estimY<-8.0) ){std::cout<<"DC=16 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
302  if(DeadCR==18 && (estimX>8.0 || estimX< -9.0 || estimY> 9.0 || estimY<-8.0) ){std::cout<<"DC=18 Position OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
303 
307 
308 
309 
310 
311 
312 
313 
314 
315 
316  // INFO from SPLINE
317  float RECOfrDcr = 0.0;
318  RECOfrDcr = SplCorr->value(DeadCR,0,50,NEWx,NEWy);
319 
320 
321 
322  // RESOLUTIONS PER AREA
323  for (int isub=0 ; isub<NSUBS ; isub++){
324  if( NEWx>XLOW[isub] && NEWx<XHIG[isub] && NEWy>YLOW[isub] && NEWy<YHIG[isub] ){
325 
326 
327  //TEST THE IDEA OF EXPs
328  if(DeadCR==6 && (isub==0 || isub ==1 || isub==5) ){
329  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
330  }
331  if(DeadCR==8 && (isub==4 || isub ==3 || isub==9) ){
332  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
333  }
334  if(DeadCR==16 && (isub==15 || isub ==21 || isub==20) ){
335  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
336  }
337  if(DeadCR==18 && (isub==19 || isub ==23 || isub==24) ){
338  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
339  }
340  if(DeadCR==7 || DeadCR==11 || DeadCR==12 || DeadCR==13 || DeadCR==17 ){
341  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
342  }
343 
344 
345  }
346  }
347 
348 
349  double ESTIMATED_ENERGY = RECOfrDcr*SUM24;
350  std::cout<<" THE ESTIMATED RECOfrDcr is : "<< RECOfrDcr <<std::endl;
351  std::cout<<" THE ESTIMATED ENERGY is : "<< ESTIMATED_ENERGY <<std::endl;
352 
353  if(RECOfrDcr<0.0){std::cout<<"NEGATIVE RECOfrDr I Return 0.0"<<std::endl; RECOfrDcr=0.0;}
354  if( (DeadCR==6 || DeadCR==8 || DeadCR==16 || DeadCR==18) && RECOfrDcr>0.20){std::cout<<"Fraction OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
355  if( (DeadCR==7 || DeadCR==17 || DeadCR==11 || DeadCR==13) && RECOfrDcr>1.00){std::cout<<"Fraction OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
356  if( DeadCR==12 && RECOfrDcr>5.00){std::cout<<"Fraction OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
357 
358 
360 // Using the findings above I build the position and the energy again!
362 
363 
364 
365  // THIS IS THE FINAL ANSWER
366  return ESTIMATED_ENERGY;
367 }
PositionCorrector::CORRX
double CORRX(int DeadCrystal, int DeadCrystalEta, int estimE, double estimX)
ExponCorrector::value
double value(int DeadCrystal, int DeadCrystalEta, int estimE, int subRegion, double estimX, double estimY)
SplineCorrector::value
double value(int DeadCrystal, int DeadCrystalEta, int estimE, double estimX, double estimY)
j
int j
Definition: DBlmapReader.cc:9
PositionCorrector::CORRY
double CORRY(int DeadCrystal, int DeadCrystalEta, int estimE, double estimY)
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
epsilon
const double epsilon
Definition: CosmicMuonParameters.h:70
Generated for CMSSW Reference Manual by doxygen 1.8.5