#include <CaloCellNeighborsAverageCorr.h>

Inheritance diagram for CaloCellNeighborsAverageCorr:

Collaboration diagram for CaloCellNeighborsAverageCorr:

Public Member Functions
	CaloCellNeighborsAverageCorr (const std::string &type, const std::string &name, const IInterface *parent)

virtual	~CaloCellNeighborsAverageCorr ()

virtual StatusCode	initialize () override
	initialize method More...

virtual StatusCode	process (CaloCellContainer *theCellContainer, const EventContext &ctx) const override
	process calo cell collection to apply corrections More...

Private Attributes
const CaloCell_ID *	m_calo_id

const TileID *	m_tile_id

bool	m_testMode

bool	m_skipDeadFeb

bool	m_skipDeadLAr

bool	m_skipDeadDrawer

bool	m_skipDeadTile

Detailed Description

Definition at line 19 of file CaloCellNeighborsAverageCorr.h.

Constructor & Destructor Documentation

◆ CaloCellNeighborsAverageCorr()

CaloCellNeighborsAverageCorr::CaloCellNeighborsAverageCorr	(	const std::string &	type,
		const std::string &	name,
		const IInterface *	parent
	)

Definition at line 28 of file CaloCellNeighborsAverageCorr.cxx.

   :base_class(type, name, parent),
    m_calo_id(nullptr),
    m_tile_id(nullptr)
 {
   declareProperty("testMode",m_testMode=false,"test mode");
   declareProperty("skipDeadFeb",    m_skipDeadFeb=true,     "Skip dead LAr Febs (default = true)");
   declareProperty("skipDeadLAr",    m_skipDeadLAr=false,    "Skip all dead LAr cells (default = false)");
   declareProperty("skipDeadDrawer", m_skipDeadDrawer=false, "Skip dead Tile drawers (default = false)");
   declareProperty("skipDeadTile",   m_skipDeadTile=true,    "Skip all dead Tile cells (default = true)");
 }

◆ ~CaloCellNeighborsAverageCorr()

virtual CaloCellNeighborsAverageCorr::~CaloCellNeighborsAverageCorr ( )

inlinevirtual

Definition at line 29 of file CaloCellNeighborsAverageCorr.h.

29 {};

Member Function Documentation

◆ initialize()

StatusCode CaloCellNeighborsAverageCorr::initialize ( )

overridevirtual

initialize method

Definition at line 46 of file CaloCellNeighborsAverageCorr.cxx.

 {
   ATH_MSG_INFO ( " in CaloCellNeighborsAverageCorr::initialize() " );
  
   if (m_testMode)
     ATH_MSG_DEBUG ( " Running in Test Mode. " );
  
   // don't need to check status of the drawer if all dead cells are ignored
   if (m_skipDeadTile && m_skipDeadDrawer) {
     ATH_MSG_WARNING ( "Please, don't set skipDeadDrawer and skipDeadTile to True at the same time" );
     ATH_MSG_WARNING ( "Setting skipDeadDrawer=False" );
     m_skipDeadDrawer=false; 
   }
   
   ATH_MSG_INFO ( "Skip Dead Feb    = " << ((m_skipDeadFeb)?"true":"false")    );
   ATH_MSG_INFO ( "Skip Dead LAr    = " << ((m_skipDeadLAr)?"true":"false")    );
   ATH_MSG_INFO ( "Skip Dead Drawer = " << ((m_skipDeadDrawer)?"true":"false") );
   ATH_MSG_INFO ( "Skip Dead Tile   = " << ((m_skipDeadTile)?"true":"false")   );
  
   
   ATH_CHECK( detStore()->retrieve (m_calo_id, "CaloCell_ID") );
   m_tile_id   = m_calo_id->tile_idHelper();
  
   ATH_MSG_INFO ( "CaloCellNeighborsAverageCorr initialize() end" );
   
   return StatusCode::SUCCESS;
 }

◆ process()

StatusCode CaloCellNeighborsAverageCorr::process	(	CaloCellContainer *	theCellContainer,
		const EventContext &	ctx
	)		const

overridevirtual

process calo cell collection to apply corrections

Definition at line 77 of file CaloCellNeighborsAverageCorr.cxx.

 {
   ATH_MSG_VERBOSE ( " in process " );
  
   std::vector<IdentifierHash> theNeighbors;
   theNeighbors.reserve(22);
   
   std::map<const CaloCell*,float> cannedUncorrectedEnergies;
   if (m_testMode) { //fill in the map of uncorrected energies, to preserve it for later use.
     CaloCellContainer::const_iterator itrCell_test=theCont->begin();
     for (; itrCell_test!=theCont->end();++itrCell_test){
       const CaloCell* pointerToCell=*itrCell_test; //yes, this is an iterator that when dereferenced returns a pointer, not a CaloCell directly.
       cannedUncorrectedEnergies[pointerToCell] = pointerToCell->energy();
       //std::cout << "GEORGIOS DEBUGGING For CaloCell* = " << pointerToCell << " we have energy = " << cannedUncorrectedEnergies[pointerToCell] << std::endl;
     }
   }
  
   int nGoodCellsPerDrawer[4][64];
  
   if (m_skipDeadDrawer) {
     memset(nGoodCellsPerDrawer,0,sizeof(nGoodCellsPerDrawer));
       
     size_t cellItr = theCont->indexFirstCellCalo(CaloCell_ID::TILE);
     size_t lastCell = theCont->indexLastCellCalo(CaloCell_ID::TILE);
  
     for (; cellItr <= lastCell; ++cellItr) {
       const CaloCell* CCell = (*theCont)[cellItr];
       if (CCell && ! CCell->badcell() ) {
         Identifier id = CCell->ID();
         int part = std::min(m_tile_id->section(id),2)+m_tile_id->side(id);
         int module = m_tile_id->module(id);
         ++nGoodCellsPerDrawer[part][module];
       }
     }
     if (msgLvl(MSG::VERBOSE)) {
       //                       1-1   2-1   1+1   2+1
       const char * name[4] = {"LBC","EBC","LBA","EBA"};
       for (int part=0; part<4; ++part) {
         for (int module=0; module<64; ++module) {
           //log << MSG::VERBOSE << name[part] << module+1 << " nCells " 
           //    << nGoodCellsPerDrawer[part][module]
           //    << ( ( nGoodCellsPerDrawer[part][module] < 2 ) ? " bad drawer" : "")
           //    << endmsg;
           if ( nGoodCellsPerDrawer[part][module] < 2 ) {
             ATH_MSG_VERBOSE ( "Bad drawer " << name[part] << module+1  );
           }
         }
       }
     }
   }
   
   CaloCellContainer::iterator itrCell=theCont->begin();
   for (; itrCell!=theCont->end();++itrCell){
     CaloCell * aCell=*itrCell;
     
     if (aCell->badcell() || m_testMode) {
  
       const CaloDetDescrElement* caloDDE = aCell->caloDDE();
       if (!caloDDE) continue;
       int theCellSubCalo = caloDDE->getSubCalo();
  
       bool isTile = (theCellSubCalo == CaloCell_ID::TILE);
  
       if (isTile) {
  
         if (m_skipDeadTile) {
           ATH_MSG_VERBOSE ( " skipping Tile hash " << m_calo_id->calo_cell_hash(aCell->ID()) );
           continue;
         } else if (m_skipDeadDrawer) {
  
           Identifier id = aCell->ID();
           int part = std::min(m_tile_id->section(id),2)+m_tile_id->side(id);
           int module = m_tile_id->module(id);
           if ( nGoodCellsPerDrawer[part][module] < 2 ) {
             ATH_MSG_VERBOSE ( " skipping Tile drawer hash " << m_calo_id->calo_cell_hash(aCell->ID()) );
             continue;
           }
         }
         
       } else {
  
         if (m_skipDeadLAr) {
           ATH_MSG_VERBOSE ( " skipping LAr hash " << m_calo_id->calo_cell_hash(aCell->ID()) );
           continue;
         } else if (m_skipDeadFeb && ((aCell->provenance() & 0x0200) == 0x0200) ) {
           ATH_MSG_VERBOSE ( " skipping LAr Feb hash " << m_calo_id->calo_cell_hash(aCell->ID()) );
           continue;
         }
       }
       
       //inspired from http://alxr.usatlas.bnl.gov/lxr-stb3/source/atlas/Calorimeter/CaloRec/src/CaloTopoClusterMaker.cxx#585
       //We need a IdentifierHash to pass as input to the get_neighbors().
       //   IdentifierHash theCellHashID = theCell->getID(); //this doesn't work (any more?)
       Identifier theCellID = aCell->ID();
       const float oldE=aCell->energy();
       IdentifierHash theCellHashID = m_calo_id->calo_cell_hash(theCellID); 
       
       //Find now the neighbors around theCell, and store them in theNeighbors vector.
       m_calo_id->get_neighbours(theCellHashID,LArNeighbours::all2D,theNeighbors);
       
       //std::cout << "Found " << theNeighbors.size() << " neighbors." << std::endl;
       
       //first let's find the volume of theCell we are correcting.
       float volumeOfTheCell=0;
       if (caloDDE) volumeOfTheCell = caloDDE->volume();
       //     std::cout << " volumeOfTheCell " << volumeOfTheCell << std::endl;
       if (volumeOfTheCell==0) continue;
       //int theCellSampling = caloDDE->getSampling();
       //      if (theCellSubCalo == CaloCell_ID::TILE) {
       //      std::cout << "theCell = " << theCellSubCalo << "  " << theCellSampling  << " is at eta=" << aCell->eta() << " phi=" << aCell->phi() << " E=" << aCell->energy() << " V=" << volumeOfTheCell*1e-6 << "D=" << aCell->energy() / volumeOfTheCell * 1e6 << std::endl;
       //      }
     
  
       const Identifier theCellRegion=m_calo_id->region_id(theCellID);
       //loop through neighbors, and calculate average energy density of guys who have a legitimate volume (namely >0).
       //float totalEnergyDensity=0;
       //float legitimateNeighbors=0;
  
       float goodNeighborEnergyDensitySum=0;
       unsigned goodNeighbors=0;
  
       float betterNeighborEnergyDensitySum=0;
       unsigned betterNeighbors=0;
  
       ATH_MSG_VERBOSE("Cell " << theCellID.get_identifier32().get_compact() << " has " << theNeighbors.size() << " neighbors");
       for (unsigned int iN=0;iN<theNeighbors.size();iN++) {
     const CaloCell* thisNeighbor = theCont->findCell(theNeighbors[iN]);
     //if there is a hardware problem in real data, then the above pointer may be NULL.  In that case the geometric neighbor is absent from the container, so we should move on to the next geometric neighbor.
     if (!thisNeighbor) continue;
     const CaloDetDescrElement* thisNeighborDDE = (thisNeighbor)->caloDDE();     
     float thisEnergy = thisNeighbor->energy();
     if (m_testMode) {
       thisEnergy = cannedUncorrectedEnergies[thisNeighbor];
       //      std::cout << "GEORGIOS DEBUGGING Retrieving cannedUncorrectedEnergies[" << thisNeighbor << "] = " << cannedUncorrectedEnergies[thisNeighbor] << std::endl;
     }
     if (thisNeighbor->badcell()) {
       ATH_MSG_VERBOSE("Ignoring neighbor " << thisNeighbor->ID().get_identifier32().get_compact() << " because of it's bad cell status");
       continue;
     }
     //int thisNeighborSubCalo = thisNeighborDDE->getSubCalo();
     //int thisNeighborSampling = thisNeighborDDE->getSampling();
     //if (thisNeighborSubCalo != theCellSubCalo) continue;
     //if (thisNeighborSampling != theCellSampling) continue; //if the quality of the cell is very different, it's a wrong idea that dE/dV would be similar.
  
  
     float thisVolume = thisNeighborDDE->volume();
     if (thisVolume <= 0) continue;
  
     //A good neightbor if we arrive at this point
     goodNeighbors++;
     float thisEnergyDensity= thisEnergy / thisVolume;
     goodNeighborEnergyDensitySum += thisEnergyDensity;
  
  
     const Identifier thisNeighborRegion=m_calo_id->region_id(thisNeighbor->ID());
     //Better neighbors are in the same region
     if (thisNeighborRegion==theCellRegion) {
       ATH_MSG_VERBOSE("Neighbor is in different region " << thisNeighborRegion.get_identifier32().get_compact() << "/" << theCellRegion.get_identifier32().get_compact());
       betterNeighbors++;
       betterNeighborEnergyDensitySum+=thisEnergyDensity;
     }
        
     //  if (theCellSubCalo == CaloCell_ID::TILE && theCellSampling == CaloCell_ID::TileBar0)
     //  std::cout << "Neighbor " << iN << " : " << thisNeighborSubCalo << " , " << thisNeighborSampling << " : " << thisNeighbor->eta() << " , " << thisNeighbor->phi() << " E=" << thisNeighbor->energy() << " V=" << thisVolume*1e-6 << " D=" << thisEnergyDensity*1e6 << std::endl;
  
       } //end loop over neighbors
  
  
       unsigned nNeighbors=0;
       float neighborEnergyDensitySum=0;
  
       if (betterNeighbors>=2) {//Have at least 2 good neighbors in the same region
     nNeighbors=betterNeighbors;
     neighborEnergyDensitySum=betterNeighborEnergyDensitySum;
       }
       else {//No good neighbors in the same region
     ATH_MSG_DEBUG("Cell " <<  theCellID.get_identifier32().get_compact() << ": Did not find enough good neighbors in the same region. Will use neighbors from different regions");
     nNeighbors=goodNeighbors;
     neighborEnergyDensitySum=goodNeighborEnergyDensitySum;
       }
  
       if(nNeighbors <= 0) {
     ATH_MSG_DEBUG("Did not get any suitable neighbor for cell " << theCellID.get_identifier32().get_compact() << ". Not corrected.");
     continue;
       }
  
       const float averageEnergyDensity = neighborEnergyDensitySum/nNeighbors;
       
       //now use the average energy density to make a prediction for the energy of theCell
       const float predictedEnergy = averageEnergyDensity * volumeOfTheCell;
       aCell->setEnergy(predictedEnergy);
       ATH_MSG_VERBOSE ( " correcting " << ((isTile)?"Tile":"LAr") <<  " id " << theCellID.get_identifier32().get_compact() << " Eold=" 
             <<oldE << "Enew=" << predictedEnergy << ", used " << nNeighbors <<  " neighbors" );
     }  // end of if(badcell)
   }  // end loop over cells
   return StatusCode::SUCCESS;
 }