#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 CaloSubdetectorTopology *	topology_

Detailed Description

template<typename DetIdT>
class EcalDeadChannelRecoveryNN< DetIdT >

Definition at line 20 of file EcalDeadChannelRecoveryNN.h.

Member Enumeration Documentation

template<typename DetIdT>

enum EcalDeadChannelRecoveryNN::CellID

Enumerator
CC
UU
DD
LL
LU
LD
RR
RU
RD

Definition at line 36 of file EcalDeadChannelRecoveryNN.h.

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

Constructor & Destructor Documentation

template<typename T >

EcalDeadChannelRecoveryNN< T >::EcalDeadChannelRecoveryNN ( )

Definition at line 11 of file EcalDeadChannelRecoveryNN.cc.

References svgfig::load(), and NULL.

                                                         {
   for (int id = 0; id < 9; ++id) {
     ctx_[id].mlp = NULL;
   }
 
   this->load();
 }

template<typename T >

EcalDeadChannelRecoveryNN< T >::~EcalDeadChannelRecoveryNN ( )

Definition at line 20 of file EcalDeadChannelRecoveryNN.cc.

                                                          {
   for (int id = 0; id < 9; ++id) {
     if (ctx_[id].mlp) {
       // @TODO segfaults for an uknown reason
       // delete ctx[id].mlp;
       // delete ctx[id].tree;
     }
   }
 }

Member Function Documentation

template<typename T >

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

Definition at line 177 of file EcalDeadChannelRecoveryNN.cc.

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

                                                                                      {
   int missing[9];
   int missing_index = 0;
 
   for (int i = 0; i < 9; i++) {
     if (TMath::Abs(M3x3Input[i]) < epsilon) {
       missing[missing_index++] = i;
     } else {
       //  Generally the "dead" cells are allowed to have negative energies (since they will be estimated by the ANN anyway).
       //  But all the remaining "live" ones must have positive values otherwise the logarithm fails.
 
       if (M3x3Input[i] < 0.0) { return -2000000.0; }
     }
   }
 
   //  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.
   int idxDC = -1 ;
   if (missing_index == 0) { return -1000000.0; }    //  Zero DC's were detected
   if (missing_index  > 1) { return -1000001.0; }    //  More than one DC's detected.
   if (missing_index == 1) { idxDC = missing[0]; } 
 
   // Arrange inputs into an array of 8, excluding the dead cell;
   int input_order[9] = { CC, RR, LL, UU, DD, RU, RD, LU, LD };
   int input_index = 0;
   Double_t input[8];
 
   for (int id : input_order) {
     if (id == idxDC)
       continue;
 
     input[input_index++] = TMath::Log(M3x3Input[id]);
   }
 
   //  Select the case to apply the appropriate NN and return the result.
   M3x3Input[idxDC] = TMath::Exp(ctx_[idxDC].mlp->Evaluate(0, input));
   return M3x3Input[idxDC];
 }

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 43 of file EcalDeadChannelRecoveryNN.cc.

References edm::FileInPath::fullPath(), EcalDeadChannelRecoveryNN< DetIdT >::MultiLayerPerceptronContext::mlp, eostools::move(), fed_dqm_sourceclient-live_cfg::path, AlCaHLTBitMon_QueryRunRegistry::string, t, EcalDeadChannelRecoveryNN< DetIdT >::MultiLayerPerceptronContext::tmp, and EcalDeadChannelRecoveryNN< DetIdT >::MultiLayerPerceptronContext::tree.

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

                                                                                            {
   std::string path = edm::FileInPath(fn).fullPath();
 
   auto t = std::make_unique<TTree>("t", "dummy MLP tree");
   t->SetDirectory(0);
 
   t->Branch("z1", &(ctx.tmp[0]), "z1/D");
   t->Branch("z2", &(ctx.tmp[1]), "z2/D");
   t->Branch("z3", &(ctx.tmp[2]), "z3/D");
   t->Branch("z4", &(ctx.tmp[3]), "z4/D");
   t->Branch("z5", &(ctx.tmp[4]), "z5/D");
   t->Branch("z6", &(ctx.tmp[5]), "z6/D");
   t->Branch("z7", &(ctx.tmp[6]), "z7/D");
   t->Branch("z8", &(ctx.tmp[7]), "z8/D");
   t->Branch("zf", &(ctx.tmp[8]), "zf/D");
 
   ctx.tree = std::move(t);
   //The TMultiLayerPerceptron will create a TEventList
   // and it will attach itself to the gDirectory
   // which will be some random TFile in the job
   // Need to reset gDirectory back to what it was to
   // avoid causing other code in other modules from
   // crashing.
   {
     auto resetGDirectory = [](TDirectory* oldDir) {gDirectory = oldDir;};
     std::unique_ptr<TDirectory, decltype(resetGDirectory)> guard(gDirectory, resetGDirectory);
     gDirectory = nullptr;
     ctx.mlp =
       std::make_unique<TMultiLayerPerceptron>("@z1,@z2,@z3,@z4,@z5,@z6,@z7,@z8:10:5:zf", ctx.tree.get());
   }
   ctx.mlp->LoadWeights(path.c_str());
 }

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 216 of file EcalDeadChannelRecoveryNN.cc.

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

                                                                                                                                                             {
   //  Since ANN corrects within a 3x3 window, the possible candidate 3x3 windows that contain 
   //  the "dead" crystal form a 5x5 window around it (totaly eight 3x3 windows overlapping).
   //  Get this 5x5 and locate the Max.Contain.Crystal within.
 
   //  Get the 5x5 window around the "dead" crystal -> vector "NxNaroundDC"
   std::vector<DetId> NxNaroundDC = topology_->getWindow(itID, 5, 5);
 
   DetIdT CellMax ;     //  Create a null DetId
   double EnergyMax = 0.0;
 
   //  Loop over all cells in the vector "NxNaroundDC", and for each cell find it's energy
   //  (from the EcalRecHits collection). Use this energy to detect the Max.Cont.Crystal.
   std::vector<DetId>::const_iterator theCells;
 
   for (theCells = NxNaroundDC.begin(); theCells != NxNaroundDC.end(); ++theCells) {
     DetIdT cell = DetIdT(*theCells);
 
     if (! cell.null()) {
       EcalRecHitCollection::const_iterator goS_it = hit_collection.find(cell);
       
       if ( goS_it !=  hit_collection.end() && goS_it->energy() >= EnergyMax ) {
         EnergyMax = goS_it->energy();
         CellMax = cell;
       }
     }
   }
 
   //  No Max.Cont.Crystal found, return back with no changes.
   if ( CellMax.null() ) { *AcceptFlag=false ; return 0.0 ; }
 
 #if 1
   // Not a smart check, because why not just get 4x4 matrix and have a guaranteed hit?
 
   //  Check that the "dead" crystal belongs to the 3x3 around  Max.Cont.Crystal
   bool dcIn3x3 = false ;
   
   std::vector<DetId> NxNaroundMaxCont = topology_->getWindow(CellMax,3,3);
   std::vector<DetId>::const_iterator testCell;
   for (testCell = NxNaroundMaxCont.begin(); testCell != NxNaroundMaxCont.end(); ++testCell) {
     if ( itID == DetIdT(*testCell) ) { dcIn3x3 = true ; } 
   }
   
   //  If the "dead" crystal is outside the 3x3 then do nothing.
   if (!dcIn3x3) { *AcceptFlag=false ; return 0.0 ; }
 #endif
   
   return makeNxNMatrice_RelDC(CellMax, hit_collection, MNxN_RelMC, AcceptFlag);
 }

template<typename T>

double EcalDeadChannelRecoveryNN< T >::recover	(	const T	id,
		const EcalRecHitCollection &	hit_collection,
		double	Sum8Cut,
		bool *	AcceptFlag
	)

Definition at line 107 of file EcalDeadChannelRecoveryNN.cc.

References edm::false, and i.

 {
   // use the correct probability matrix
   typedef class EcalCrystalMatrixProbality<T> P;
 
   double NewEnergy = 0.0;
 
   double NewEnergy_RelMC = 0.0;
   double NewEnergy_RelDC = 0.0;
 
   double MNxN_RelMC[9] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ;
   double MNxN_RelDC[9] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ;
 
   double sum8 = 0.0;
 
   double sum8_RelMC = makeNxNMatrice_RelMC(id, hit_collection, MNxN_RelMC,AcceptFlag);
   double sum8_RelDC = makeNxNMatrice_RelDC(id, hit_collection, MNxN_RelDC,AcceptFlag);
 
   //  Only if "AcceptFlag" is true call the ANN
   if ( *AcceptFlag ) {
     if (sum8_RelDC > Sum8Cut && sum8_RelMC > Sum8Cut) {
       NewEnergy_RelMC = estimateEnergy(MNxN_RelMC);
       NewEnergy_RelDC = estimateEnergy(MNxN_RelDC);
       
       //  Matrices "MNxN_RelMC" and "MNxN_RelDC" have now the full set of energies, the original ones plus 
       //  whatever "estimates" by the ANN for the "dead" xtal. Use those full matrices and calculate probabilities.
       //  
       double SumMNxN_RelMC =  MNxN_RelMC[LU] + MNxN_RelMC[UU] + MNxN_RelMC[RU] + 
                               MNxN_RelMC[LL] + MNxN_RelMC[CC] + MNxN_RelMC[RR] + 
                               MNxN_RelMC[LD] + MNxN_RelMC[DD] + MNxN_RelMC[RD] ;
       
       double frMNxN_RelMC[9];  for (int i=0; i<9; i++) { frMNxN_RelMC[i] = MNxN_RelMC[i] / SumMNxN_RelMC ; }
       
       double prMNxN_RelMC  =  P::Diagonal(  frMNxN_RelMC[LU] ) * P::UpDown(  frMNxN_RelMC[UU] ) * P::Diagonal(  frMNxN_RelMC[RU] ) * 
                               P::ReftRight( frMNxN_RelMC[LL] ) * P::Central( frMNxN_RelMC[CC] ) * P::ReftRight( frMNxN_RelMC[RR] ) * 
                               P::Diagonal(  frMNxN_RelMC[LD] ) * P::UpDown(  frMNxN_RelMC[DD] ) * P::Diagonal(  frMNxN_RelMC[RD] ) ;
       
       double SumMNxN_RelDC =  MNxN_RelDC[LU] + MNxN_RelDC[UU] + MNxN_RelDC[RU] + 
                               MNxN_RelDC[LL] + MNxN_RelDC[CC] + MNxN_RelDC[RR] + 
                               MNxN_RelDC[LD] + MNxN_RelDC[DD] + MNxN_RelDC[RD] ;
       
       double frMNxN_RelDC[9];  for (int i=0; i<9; i++) { frMNxN_RelDC[i] = MNxN_RelDC[i] / SumMNxN_RelDC ; }
       
       double prMNxN_RelDC  =  P::Diagonal(  frMNxN_RelDC[LU] ) * P::UpDown(  frMNxN_RelDC[UU] ) * P::Diagonal(  frMNxN_RelDC[RU] ) * 
                               P::ReftRight( frMNxN_RelDC[LL] ) * P::Central( frMNxN_RelDC[CC] ) * P::ReftRight( frMNxN_RelDC[RR] ) * 
                               P::Diagonal(  frMNxN_RelDC[LD] ) * P::UpDown(  frMNxN_RelDC[DD] ) * P::Diagonal(  frMNxN_RelDC[RD] ) ;
       
       if ( prMNxN_RelDC > prMNxN_RelMC )  { NewEnergy = NewEnergy_RelDC ; sum8 = sum8_RelDC ; } 
       if ( prMNxN_RelDC <= prMNxN_RelMC ) { NewEnergy = NewEnergy_RelMC ; sum8 = sum8_RelMC ; } 
       
       
       //  If the return value of "CorrectDeadChannelsNN" is one of the followin negative values then
       //  it corresponds to an error condition. See "CorrectDeadChannelsNN.cc" for possible values.
       if ( NewEnergy == -1000000.0 ||
            NewEnergy == -1000001.0 ||
            NewEnergy == -2000000.0 ||
            NewEnergy == -3000000.0 ||
            NewEnergy == -3000001.0 ) { *AcceptFlag=false ; NewEnergy = 0.0 ; }             
     }
   }
   
   // Protect against non physical high values
   // From the distribution of (max.cont.xtal / Sum8) we get as limit 5 (hard) and 10 (softer)
   // Choose 10 as highest possible energy to be assigned to the dead channel under any scenario.
   if ( NewEnergy > 10.0 * sum8 ) { *AcceptFlag=false ; NewEnergy = 0.0 ; }
   
   return NewEnergy;
 }

template<typename DetIdT>

double EcalDeadChannelRecoveryNN< DetIdT >::reorderMxNMatrix	(	EBDetId	it,
		const std::vector< DetId > &	window,
		const EcalRecHitCollection &	hit_collection,
		double *	MNxN,
		bool *	AcceptFlag
	)