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(), ExponCorrector::value(), and w().

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  int NSUBS = 25;
165  for(int ix=0; ix<5 ; ix++){
166  for(int iy=0; iy<5 ; iy++){
167  int isub= ix*5+iy ;
168 
169  XLOW[isub]= -10.0 +float(ix)*4.0;
170  XHIG[isub]= XLOW[isub] + 4.0;
171  YLOW[isub]= -10.0 +float(iy)*4.0;;
172  YHIG[isub]= YLOW[isub] + 4.0;
173 
174  }
175  }
176 
177 
178 
179 
181 
182 
183 
184  //First Find the position of the Dead Channel in 3x3 matrix
185  int DeadCR = -1;
186  for(int ix=1; ix<4 ; ix++){
187  for(int iy=1; iy<4 ; iy++){
188  int idx = ix*5+iy;
189  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;}
190  if(fabs(crE[idx])<epsilon && DeadCR==-1)DeadCR=idx;
191  }
192  }
193 
194 
195  //std::cout<<" THE DEAD CHANNEL IS : "<< DeadCR <<std::endl;
196  SUM24=0.0;
197  for(int j=0;j<25;j++)SUM24+=crE[j];
198  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;}
199 
200 
201 
202  // CR=12 must be the maximum in 5x5 !!!
203  int CisMax; CisMax=-1;
204  float MaxEin5x5; MaxEin5x5=0.0;
205  for(int ic=0;ic<24;ic++){
206  if(crE[ic]>MaxEin5x5){
207  MaxEin5x5 = crE[ic];
208  CisMax = ic ;
209  }
210  }
211  if(CisMax != 12){
212  std::cout<<"ERROR ----> Central has NOT the MAX energy in 5x5 ... I return 0.0"<<std::endl;
213  return 0.0;
214  }
215 
216 
217 
218  //AVOID BREM or OTHER DEPOSITIONS IN 5x5
219  //std::cout<<" CHECK FOR BREM or UNUSUAL PATTERNS ... " <<std::endl;
220  //std::cout<<" SUMuu="<<SUMuu<<",SUMu="<<SUMu<<",SUMd="<<SUMd <<",SUMdd="<<SUMdd <<",SUMll="<<SUMll <<",SUMl="<<SUMl <<",SUMrr="<<SUMrr <<",SUMr="<<SUMr <<std::endl;
221  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;}
222  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;}
223  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;}
224  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;}
225 
226  if(DeadCR==7 && (SUMuu>SUMu || SUMdd>SUMd || SUMrr>SUMr)){std::cout<<"Unusual Pattern in 7 I return 0.0"<<std::endl; return 0.0;}
227  if(DeadCR==17 && (SUMuu>SUMu || SUMdd>SUMd || SUMll>SUMl)){std::cout<<"Unusual Pattern in 17 I return 0.0"<<std::endl; return 0.0;}
228  if(DeadCR==11 && (SUMll>SUMl || SUMrr>SUMr || SUMuu>SUMu)){std::cout<<"Unusual Pattern in 11 I return 0.0"<<std::endl; return 0.0;}
229  if(DeadCR==13 && (SUMll>SUMl || SUMrr>SUMr || SUMdd>SUMd)){std::cout<<"Unusual Pattern in 13 I return 0.0"<<std::endl; return 0.0;}
230 
231 
232  //std::cout<<" NO BREM OR OTHER DEPOSITIONS IN THIS PATTERN --- I continue ... " <<std::endl;
233 
234  SUM24=0.0;
235  for(int j=0;j<25;j++)if(j!=DeadCR)SUM24+=crE[j];
236  // std::cout<<" THE Enorm is : "<< SUM24 <<std::endl;
237 
238 
239  //CHECK IF IT IS REALLY THE CORRECT DEAD CHANNEL
240  // This will be included in the next version
241 
242 
243 
246  // Estimate X,Y
249 
250  NEWx=0.0; NEWy=0.0;
251 
252 
253  estimX=0.0;
254  estimY=0.0;
255  SUMlogFR=0.0;
256  for(int ix=0; ix<5 ; ix++){
257  for(int iy=0; iy<5 ; iy++){
258  int idx = ix*5+iy;
259 
260  float xpos = 20.0*(float(ix)-2.0);
261  float ypos = 20.0*(float(iy)-2.0);
262 
263  if(idx != DeadCR){
264 
265  // weight definition
266  float w = crE[idx]/SUM24;
267 
268  if(w<0.0)w=0.0;
269  SUMlogFR = SUMlogFR + w;
270 
271  estimX = estimX + xpos*w;
272  estimY = estimY + ypos*w;
273  }
274 
275  } // iy loop
276  }// ix loop
277 
278  estimX = estimX/SUMlogFR;
279  estimY = estimY/SUMlogFR;
280 
281 
282  NEWx = PosCorr->CORRX(DeadCR,0,50,estimX);
283  NEWy = PosCorr->CORRY(DeadCR,0,50,estimY);
284 
285 
286  //std::cout<<" FINAL X,Y calculation: DC , estimX, estimY, NEWx, NEWy : "<<DeadCR<<" "<<estimX<<" "<<estimY<<" "<< NEWx<<" "<<NEWy<<std::endl;
287 
288 
289  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;}
290  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;}
291  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;}
292  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;}
293 
294  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;}
295 
296  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;}
297  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;}
298  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;}
299  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;}
300 
304 
305 
306 
307 
308 
309 
310 
311 
312 
313  // INFO from SPLINE
314  float RECOfrDcr = 0.0;
315  RECOfrDcr = SplCorr->value(DeadCR,0,50,NEWx,NEWy);
316 
317 
318 
319  // RESOLUTIONS PER AREA
320  for (int isub=0 ; isub<NSUBS ; isub++){
321  if( NEWx>XLOW[isub] && NEWx<XHIG[isub] && NEWy>YLOW[isub] && NEWy<YHIG[isub] ){
322 
323 
324  //TEST THE IDEA OF EXPs
325  if(DeadCR==6 && (isub==0 || isub ==1 || isub==5) ){
326  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
327  }
328  if(DeadCR==8 && (isub==4 || isub ==3 || isub==9) ){
329  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
330  }
331  if(DeadCR==16 && (isub==15 || isub ==21 || isub==20) ){
332  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
333  }
334  if(DeadCR==18 && (isub==19 || isub ==23 || isub==24) ){
335  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
336  }
337  if(DeadCR==7 || DeadCR==11 || DeadCR==12 || DeadCR==13 || DeadCR==17 ){
338  RECOfrDcr = ExpCorr->value(DeadCR,0,50,isub,NEWx,NEWy);
339  }
340 
341 
342  }
343  }
344 
345 
346  double ESTIMATED_ENERGY = RECOfrDcr*SUM24;
347  std::cout<<" THE ESTIMATED RECOfrDcr is : "<< RECOfrDcr <<std::endl;
348  std::cout<<" THE ESTIMATED ENERGY is : "<< ESTIMATED_ENERGY <<std::endl;
349 
350  if(RECOfrDcr<0.0){std::cout<<"NEGATIVE RECOfrDr I Return 0.0"<<std::endl; RECOfrDcr=0.0;}
351  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;}
352  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;}
353  if( DeadCR==12 && RECOfrDcr>5.00){std::cout<<"Fraction OUT of LIMIT I return 0.0"<<std::endl; return 0.0;}
354 
355 
357 // Using the findings above I build the position and the energy again!
359 
360 
361 
362  // THIS IS THE FINAL ANSWER
363  return ESTIMATED_ENERGY;
364 }
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:121
epsilon
const double epsilon
Definition: CosmicMuonParameters.h:70
w
T w() const
Definition: newBasic3DVector.h:234
Generated for CMSSW Reference Manual by doxygen 1.8.5