CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Classes | Public Types | Public Member Functions | Public Attributes | Private Member Functions | Private Attributes
EcalDeadChannelRecoveryNN< DetIdT > Class Template Reference

#include <EcalDeadChannelRecoveryNN.h>

Classes

struct  MultiLayerPerceptronContext
 

Public Types

enum  CellID {
  CC = 0, UU = 1, DD = 2, LL = 3,
  LU = 4, LD = 5, RR = 6, RU = 7,
  RD = 8
}
 

Public Member Functions

 EcalDeadChannelRecoveryNN ()
 
double estimateEnergy (double *M3x3Input, double epsilon=0.0000001)
 
double makeNxNMatrice_RelDC (DetIdT itID, const EcalRecHitCollection &hit_collection, double *MNxN_RelDC, bool *AccFlag)
 
double makeNxNMatrice_RelMC (DetIdT itID, const EcalRecHitCollection &hit_collection, double *MNxN_RelMC, bool *AccFlag)
 
double recover (const DetIdT id, const EcalRecHitCollection &hit_collection, double Sum8Cut, bool *AcceptFlag)
 
double reorderMxNMatrix (EBDetId it, const std::vector< DetId > &window, const EcalRecHitCollection &hit_collection, double *MNxN, bool *AcceptFlag)
 
void setCaloTopology (const CaloTopology *topo)
 
 ~EcalDeadChannelRecoveryNN ()
 

Public Attributes

const int CellX [9]
 
const int CellY [9]
 

Private Member Functions

void load ()
 
void load_file (MultiLayerPerceptronContext &ctx, std::string fn)
 

Private Attributes

MultiLayerPerceptronContext ctx_ [9]
 
const CaloSubdetectorTopologytopology_
 

Detailed Description

template<typename DetIdT>
class EcalDeadChannelRecoveryNN< DetIdT >

Definition at line 19 of file EcalDeadChannelRecoveryNN.h.

Member Enumeration Documentation

template<typename DetIdT>
enum EcalDeadChannelRecoveryNN::CellID

Constructor & Destructor Documentation

template<typename T >
EcalDeadChannelRecoveryNN< T >::EcalDeadChannelRecoveryNN ( )

Definition at line 10 of file EcalDeadChannelRecoveryNN.cc.

References svgfig::load(), and NULL.

10  {
11  for (int id = 0; id < 9; ++id) {
12  ctx_[id].mlp = NULL;
13  }
14 
15  this->load();
16 }
#define NULL
Definition: scimark2.h:8
MultiLayerPerceptronContext ctx_[9]
template<typename T >
EcalDeadChannelRecoveryNN< T >::~EcalDeadChannelRecoveryNN ( )

Definition at line 19 of file EcalDeadChannelRecoveryNN.cc.

19  {
20  for (int id = 0; id < 9; ++id) {
21  if (ctx_[id].mlp) {
22  // @TODO segfaults for an uknown reason
23  // delete ctx[id].mlp;
24  // delete ctx[id].tree;
25  }
26  }
27 }
MultiLayerPerceptronContext ctx_[9]

Member Function Documentation

template<typename T >
double EcalDeadChannelRecoveryNN< T >::estimateEnergy ( double *  M3x3Input,
double  epsilon = 0.0000001 
)

Definition at line 165 of file EcalDeadChannelRecoveryNN.cc.

References Abs(), i, input, and combine::missing.

165  {
166  int missing[9];
167  int missing_index = 0;
168 
169  for (int i = 0; i < 9; i++) {
170  if (TMath::Abs(M3x3Input[i]) < epsilon) {
171  missing[missing_index++] = i;
172  } else {
173  // Generally the "dead" cells are allowed to have negative energies (since they will be estimated by the ANN anyway).
174  // But all the remaining "live" ones must have positive values otherwise the logarithm fails.
175 
176  if (M3x3Input[i] < 0.0) { return -2000000.0; }
177  }
178  }
179 
180  // Currently EXACTLY ONE AND ONLY ONE dead cell is corrected. Return -1000000.0 if zero DC's detected and -101.0 if more than one DC's exist.
181  int idxDC = -1 ;
182  if (missing_index == 0) { return -1000000.0; } // Zero DC's were detected
183  if (missing_index > 1) { return -1000001.0; } // More than one DC's detected.
184  if (missing_index == 1) { idxDC = missing[0]; }
185 
186  // Arrange inputs into an array of 8, excluding the dead cell;
187  int input_order[9] = { CC, RR, LL, UU, DD, RU, RD, LU, LD };
188  int input_index = 0;
189  Double_t input[8];
190 
191  for (int id : input_order) {
192  if (id == idxDC)
193  continue;
194 
195  input[input_index++] = TMath::Log(M3x3Input[id]);
196  }
197 
198  // Select the case to apply the appropriate NN and return the result.
199  M3x3Input[idxDC] = TMath::Exp(ctx_[idxDC].mlp->Evaluate(0, input));
200  return M3x3Input[idxDC];
201 }
dictionary missing
Definition: combine.py:4
int i
Definition: DBlmapReader.cc:9
static std::string const input
Definition: EdmProvDump.cc:44
T Abs(T a)
Definition: MathUtil.h:49
MultiLayerPerceptronContext ctx_[9]
const double epsilon
template<typename DetIdT>
void EcalDeadChannelRecoveryNN< DetIdT >::load ( )
private
template<typename T >
void EcalDeadChannelRecoveryNN< T >::load_file ( MultiLayerPerceptronContext ctx,
std::string  fn 
)
private

Definition at line 42 of file EcalDeadChannelRecoveryNN.cc.

References edm::FileInPath::fullPath(), EcalDeadChannelRecoveryNN< DetIdT >::MultiLayerPerceptronContext::mlp, cmsHarvester::path, AlCaHLTBitMon_QueryRunRegistry::string, edmStreamStallGrapher::t, EcalDeadChannelRecoveryNN< DetIdT >::MultiLayerPerceptronContext::tmp, and EcalDeadChannelRecoveryNN< DetIdT >::MultiLayerPerceptronContext::tree.

Referenced by EcalDeadChannelRecoveryNN< EBDetId >::load(), and EcalDeadChannelRecoveryNN< EEDetId >::load().

42  {
44 
45  TTree *t = new TTree("t", "dummy MLP tree");
46  t->SetDirectory(0);
47 
48  t->Branch("z1", &(ctx.tmp[0]), "z1/D");
49  t->Branch("z2", &(ctx.tmp[1]), "z2/D");
50  t->Branch("z3", &(ctx.tmp[2]), "z3/D");
51  t->Branch("z4", &(ctx.tmp[3]), "z4/D");
52  t->Branch("z5", &(ctx.tmp[4]), "z5/D");
53  t->Branch("z6", &(ctx.tmp[5]), "z6/D");
54  t->Branch("z7", &(ctx.tmp[6]), "z7/D");
55  t->Branch("z8", &(ctx.tmp[7]), "z8/D");
56  t->Branch("zf", &(ctx.tmp[8]), "zf/D");
57 
58  ctx.tree = t;
59  ctx.mlp =
60  new TMultiLayerPerceptron("@z1,@z2,@z3,@z4,@z5,@z6,@z7,@z8:10:5:zf", t);
61  ctx.mlp->LoadWeights(path.c_str());
62 }
tuple path
else: Piece not in the list, fine.
std::string fullPath() const
Definition: FileInPath.cc:165
template<typename DetIdT>
double EcalDeadChannelRecoveryNN< DetIdT >::makeNxNMatrice_RelDC ( DetIdT  itID,
const EcalRecHitCollection hit_collection,
double *  MNxN_RelDC,
bool *  AccFlag 
)
template<typename DetIdT>
double EcalDeadChannelRecoveryNN< DetIdT >::makeNxNMatrice_RelMC ( DetIdT  itID,
const EcalRecHitCollection hit_collection,
double *  MNxN_RelMC,
bool *  AccFlag 
)

Definition at line 204 of file EcalDeadChannelRecoveryNN.cc.

References edm::SortedCollection< T, SORT >::end(), edm::false, edm::SortedCollection< T, SORT >::find(), and funct::true.

204  {
205  // Since ANN corrects within a 3x3 window, the possible candidate 3x3 windows that contain
206  // the "dead" crystal form a 5x5 window around it (totaly eight 3x3 windows overlapping).
207  // Get this 5x5 and locate the Max.Contain.Crystal within.
208 
209  // Get the 5x5 window around the "dead" crystal -> vector "NxNaroundDC"
210  std::vector<DetId> NxNaroundDC = topology_->getWindow(itID, 5, 5);
211 
212  DetIdT CellMax ; // Create a null DetId
213  double EnergyMax = 0.0;
214 
215  // Loop over all cells in the vector "NxNaroundDC", and for each cell find it's energy
216  // (from the EcalRecHits collection). Use this energy to detect the Max.Cont.Crystal.
217  std::vector<DetId>::const_iterator theCells;
218 
219  for (theCells = NxNaroundDC.begin(); theCells != NxNaroundDC.end(); ++theCells) {
220  DetIdT cell = DetIdT(*theCells);
221 
222  if (! cell.null()) {
223  EcalRecHitCollection::const_iterator goS_it = hit_collection.find(cell);
224 
225  if ( goS_it != hit_collection.end() && goS_it->energy() >= EnergyMax ) {
226  EnergyMax = goS_it->energy();
227  CellMax = cell;
228  }
229  }
230  }
231 
232  // No Max.Cont.Crystal found, return back with no changes.
233  if ( CellMax.null() ) { *AcceptFlag=false ; return 0.0 ; }
234 
235 #if 1
236  // Not a smart check, because why not just get 4x4 matrix and have a guaranteed hit?
237 
238  // Check that the "dead" crystal belongs to the 3x3 around Max.Cont.Crystal
239  bool dcIn3x3 = false ;
240 
241  std::vector<DetId> NxNaroundMaxCont = topology_->getWindow(CellMax,3,3);
242  std::vector<DetId>::const_iterator testCell;
243  for (testCell = NxNaroundMaxCont.begin(); testCell != NxNaroundMaxCont.end(); ++testCell) {
244  if ( itID == DetIdT(*testCell) ) { dcIn3x3 = true ; }
245  }
246 
247  // If the "dead" crystal is outside the 3x3 then do nothing.
248  if (!dcIn3x3) { *AcceptFlag=false ; return 0.0 ; }
249 #endif
250 
251  return makeNxNMatrice_RelDC(CellMax, hit_collection, MNxN_RelMC, AcceptFlag);
252 }
std::vector< EcalRecHit >::const_iterator const_iterator
double makeNxNMatrice_RelDC(DetIdT itID, const EcalRecHitCollection &hit_collection, double *MNxN_RelDC, bool *AccFlag)
const CaloSubdetectorTopology * topology_
const_iterator end() const
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
iterator find(key_type k)
volatile std::atomic< bool > shutdown_flag false
template<typename T>
double EcalDeadChannelRecoveryNN< T >::recover ( const T  id,
const EcalRecHitCollection hit_collection,
double  Sum8Cut,
bool *  AcceptFlag 
)

Definition at line 95 of file EcalDeadChannelRecoveryNN.cc.

References edm::false, and i.

96 {
97  // use the correct probability matrix
98  typedef class EcalCrystalMatrixProbality<T> P;
99 
100  double NewEnergy = 0.0;
101 
102  double NewEnergy_RelMC = 0.0;
103  double NewEnergy_RelDC = 0.0;
104 
105  double MNxN_RelMC[9] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ;
106  double MNxN_RelDC[9] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ;
107 
108  double sum8 = 0.0;
109 
110  double sum8_RelMC = makeNxNMatrice_RelMC(id, hit_collection, MNxN_RelMC,AcceptFlag);
111  double sum8_RelDC = makeNxNMatrice_RelDC(id, hit_collection, MNxN_RelDC,AcceptFlag);
112 
113  // Only if "AcceptFlag" is true call the ANN
114  if ( *AcceptFlag ) {
115  if (sum8_RelDC > Sum8Cut && sum8_RelMC > Sum8Cut) {
116  NewEnergy_RelMC = estimateEnergy(MNxN_RelMC);
117  NewEnergy_RelDC = estimateEnergy(MNxN_RelDC);
118 
119  // Matrices "MNxN_RelMC" and "MNxN_RelDC" have now the full set of energies, the original ones plus
120  // whatever "estimates" by the ANN for the "dead" xtal. Use those full matrices and calculate probabilities.
121  //
122  double SumMNxN_RelMC = MNxN_RelMC[LU] + MNxN_RelMC[UU] + MNxN_RelMC[RU] +
123  MNxN_RelMC[LL] + MNxN_RelMC[CC] + MNxN_RelMC[RR] +
124  MNxN_RelMC[LD] + MNxN_RelMC[DD] + MNxN_RelMC[RD] ;
125 
126  double frMNxN_RelMC[9]; for (int i=0; i<9; i++) { frMNxN_RelMC[i] = MNxN_RelMC[i] / SumMNxN_RelMC ; }
127 
128  double prMNxN_RelMC = P::Diagonal( frMNxN_RelMC[LU] ) * P::UpDown( frMNxN_RelMC[UU] ) * P::Diagonal( frMNxN_RelMC[RU] ) *
129  P::ReftRight( frMNxN_RelMC[LL] ) * P::Central( frMNxN_RelMC[CC] ) * P::ReftRight( frMNxN_RelMC[RR] ) *
130  P::Diagonal( frMNxN_RelMC[LD] ) * P::UpDown( frMNxN_RelMC[DD] ) * P::Diagonal( frMNxN_RelMC[RD] ) ;
131 
132  double SumMNxN_RelDC = MNxN_RelDC[LU] + MNxN_RelDC[UU] + MNxN_RelDC[RU] +
133  MNxN_RelDC[LL] + MNxN_RelDC[CC] + MNxN_RelDC[RR] +
134  MNxN_RelDC[LD] + MNxN_RelDC[DD] + MNxN_RelDC[RD] ;
135 
136  double frMNxN_RelDC[9]; for (int i=0; i<9; i++) { frMNxN_RelDC[i] = MNxN_RelDC[i] / SumMNxN_RelDC ; }
137 
138  double prMNxN_RelDC = P::Diagonal( frMNxN_RelDC[LU] ) * P::UpDown( frMNxN_RelDC[UU] ) * P::Diagonal( frMNxN_RelDC[RU] ) *
139  P::ReftRight( frMNxN_RelDC[LL] ) * P::Central( frMNxN_RelDC[CC] ) * P::ReftRight( frMNxN_RelDC[RR] ) *
140  P::Diagonal( frMNxN_RelDC[LD] ) * P::UpDown( frMNxN_RelDC[DD] ) * P::Diagonal( frMNxN_RelDC[RD] ) ;
141 
142  if ( prMNxN_RelDC > prMNxN_RelMC ) { NewEnergy = NewEnergy_RelDC ; sum8 = sum8_RelDC ; }
143  if ( prMNxN_RelDC <= prMNxN_RelMC ) { NewEnergy = NewEnergy_RelMC ; sum8 = sum8_RelMC ; }
144 
145 
146  // If the return value of "CorrectDeadChannelsNN" is one of the followin negative values then
147  // it corresponds to an error condition. See "CorrectDeadChannelsNN.cc" for possible values.
148  if ( NewEnergy == -1000000.0 ||
149  NewEnergy == -1000001.0 ||
150  NewEnergy == -2000000.0 ||
151  NewEnergy == -3000000.0 ||
152  NewEnergy == -3000001.0 ) { *AcceptFlag=false ; NewEnergy = 0.0 ; }
153  }
154  }
155 
156  // Protect against non physical high values
157  // From the distribution of (max.cont.xtal / Sum8) we get as limit 5 (hard) and 10 (softer)
158  // Choose 10 as highest possible energy to be assigned to the dead channel under any scenario.
159  if ( NewEnergy > 10.0 * sum8 ) { *AcceptFlag=false ; NewEnergy = 0.0 ; }
160 
161  return NewEnergy;
162 }
int i
Definition: DBlmapReader.cc:9
double makeNxNMatrice_RelMC(DetIdT itID, const EcalRecHitCollection &hit_collection, double *MNxN_RelMC, bool *AccFlag)
#define P
double makeNxNMatrice_RelDC(DetIdT itID, const EcalRecHitCollection &hit_collection, double *MNxN_RelDC, bool *AccFlag)
double estimateEnergy(double *M3x3Input, double epsilon=0.0000001)
volatile std::atomic< bool > shutdown_flag false
long double T
template<typename DetIdT>
double EcalDeadChannelRecoveryNN< DetIdT >::reorderMxNMatrix ( EBDetId  it,
const std::vector< DetId > &  window,
const EcalRecHitCollection hit_collection,
double *  MNxN,
bool *  AcceptFlag 
)
template<typename DetIdT>
void EcalDeadChannelRecoveryNN< DetIdT >::setCaloTopology ( const CaloTopology topo)

Member Data Documentation

template<typename DetIdT>
const int EcalDeadChannelRecoveryNN< DetIdT >::CellX[9]
Initial value:
= { 0, 0, 0 , -1, -1, -1 ,
1, 1, 1 }

Definition at line 51 of file EcalDeadChannelRecoveryNN.h.

template<typename DetIdT>
const int EcalDeadChannelRecoveryNN< DetIdT >::CellY[9]
Initial value:
= { 0, -1, 1 , 0, -1, 1 ,
0, -1, 1 }

Definition at line 54 of file EcalDeadChannelRecoveryNN.h.

template<typename DetIdT>
MultiLayerPerceptronContext EcalDeadChannelRecoveryNN< DetIdT >::ctx_[9]
private
template<typename DetIdT>
const CaloSubdetectorTopology* EcalDeadChannelRecoveryNN< DetIdT >::topology_
private