10#include <ecl/modules/eclLocalMaximumFinder/ECLLocalMaximumFinderModule.h>
13#include <ecl/dataobjects/ECLCalDigit.h>
14#include <ecl/dataobjects/ECLConnectedRegion.h>
15#include <ecl/dataobjects/ECLDigit.h>
16#include <ecl/dataobjects/ECLElementNumbers.h>
17#include <ecl/dataobjects/ECLHit.h>
18#include <ecl/dataobjects/ECLLocalMaximum.h>
19#include <ecl/geometry/ECLGeometryPar.h>
20#include <ecl/geometry/ECLNeighbours.h>
27#include <framework/gearbox/Const.h>
28#include <framework/logging/Logger.h>
29#include <mdst/dataobjects/MCParticle.h>
45 m_mcParticles(mcParticleArrayName()),
46 m_eclHits(eclHitArrayName()),
47 m_eclDigits(eclDigitArrayName()),
48 m_eclCalDigits(eclCalDigitArrayName()),
49 m_eclConnectedRegions(eclConnectedRegionArrayName()),
50 m_eclLocalMaximums(eclLocalMaximumArrayName())
61 "Run in training mode (i.e. fill file with MVA input variables and determine MC truth of LM.).", 0);
62 addParam(
"outfileName",
m_outfileName,
"Output file name for training file.", std::string(
"ECLLocalMaximumFinderOutput.root"));
63 addParam(
"method",
m_method,
"Method to determine the LM (cut, none, fastbdt).", std::string(
"none"));
65 addParam(
"cutOffset",
m_cutOffset,
"Cut method specific: Offset. (BaBar: 2.5, high eff: 1.40)", 2.5);
66 addParam(
"cutSlope",
m_cutSlope,
"Cut method specific: Slope. (BaBar: 0.5, high eff: 3.0)", 0.5);
78 B2DEBUG(200,
"ECLLocalMaximumFinderModule::initialize()");
89 B2WARNING(
"ECLLocalMaximumFinderModule::initialize: Energy threshold too small, resetting to " <<
c_minEnergyCut <<
" GeV");
94 B2WARNING(
"ECLLocalMaximumFinderModule::initialize: Matching fraction must be above 1e-6, input: " <<
m_truthFraction <<
95 " will be reset to 1e-6.");
100 B2WARNING(
"ECLLocalMaximumFinderModule::initialize: Matching fraction must be below 1.0, input: " <<
m_truthFraction <<
101 " will be reset to 1.0.");
118 const int cBufferLength = 500;
119 char tmpBuffer[cBufferLength];
120 int tmpBufferLength = sprintf(tmpBuffer,
"%s",
m_outfileName.c_str());
122 if (tmpBufferLength >= cBufferLength) {
123 B2FATAL(
"ECLLocalMaximumFinderModule::initialize(): Output training filename length too long!");
126 m_outfile =
new TFile(tmpBuffer,
"RECREATE");
127 m_tree =
new TTree(
"locmax",
"locmax");
168 B2DEBUG(200,
"ECLLocalMaximumFinderModule::event()");
177 std::vector< double > vNeighourEnergies;
182 const int crId = aCR.getCRId();
191 std::fill_n(vNeighourEnergies.begin(), vNeighourEnergies.size(),
198 int neighbourCount = 0;
200 if (neighbourId == aECLCalDigit.getCellId())
continue;
204 double energyNeighbour = 0.0;
207 vNeighourEnergies[neighbourCount] = energyNeighbour;
210 vNeighourEnergies[neighbourCount] = 0.0;
214 if (energyNeighbour > aECLCalDigit.getEnergy()) {
223 for (
unsigned int npos = 0; npos < vNeighourEnergies.size(); ++npos) {
224 if (vNeighourEnergies[npos] >= 0) {
230 m_time = aECLCalDigit.getTime();
236 m_energy =
static_cast < float >(aECLCalDigit.getEnergy());
237 m_cellId =
static_cast < float >(aECLCalDigit.getCellId());
240 m_CRId =
static_cast < float >(crId);
241 m_LMId =
static_cast < float >(iLM);
252 auto relatedParticlePairs = aECLCalDigit.getRelationsWith<
MCParticle>();
254 for (
unsigned int irel = 0; irel < relatedParticlePairs.size(); irel++) {
255 const auto particle = relatedParticlePairs.object(irel);
256 const double weight = relatedParticlePairs.weight(irel);
259 int motherindex = -1;
265 B2DEBUG(175,
" -> digt energy: " << aECLCalDigit.getEnergy());
298 ", aECLCalDigit.getEnergy(): " << aECLCalDigit.getEnergy());
300 B2DEBUG(200,
"(m_maxNeighbourEnergy - m_cutRatioCorrection) / (aECLCalDigit.getEnergy() - m_cutRatioCorrection)" <<
322 int highestEnergyCellId = -1;
323 double highestEnergy = 0.0;
327 if (aECLCalDigit.getEnergy() > highestEnergy) {
328 highestEnergyCellId = aECLCalDigit.getCellId();
329 highestEnergy = aECLCalDigit.getEnergy();
340 uint16_t nLMPerRegion[3] = {};
342 const int iCellId = aLM.getCellId();
359 B2DEBUG(200,
"ECLLocalMaximumFinderModule::endRun()");
365 B2DEBUG(200,
"ECLLocalMaximumFinderModule::terminate()");
385 B2DEBUG(175,
"ECLLocalMaximumFinderModule::makeLocalMaximum(): local maximum cellid: " << cellId);
388 aLocalMaximum->setLMId(lmId);
391 aLocalMaximum->setCellId(cellId);
406 for (
unsigned int i = 0; i < 10; ++i) {
407 for (
unsigned int j = 0; j < 5; ++j) {
440 int index = particle.getArrayIndex();
451 pi0index =
m_mcParticles[index]->getMother()->getArrayIndex();
465 pi0arrayindex = pi0index;
472 const double theta = vertex.Theta();
474 if (theta > 0.555015 and theta < 2.26369) {
475 double radius = vertex.Perp();
477 }
else if (theta <= 0.555015) {
479 }
else if (theta >= 2.26369) {
489 for (
int i = 0; i < 5; i++) {
507 for (
unsigned int i = 0; i < 10; ++i) {
508 for (
unsigned int j = 0; j < 5; ++j) {
537 }
else if (abs(motherpdg) ==
Const::muon.getPDGCode()) {
static const ParticleType neutron
neutron particle
static const ParticleType pi0
neutral pion particle
static const ChargedStable muon
muon particle
static const ChargedStable pion
charged pion particle
static const ParticleType Klong
K^0_L particle.
static const ChargedStable proton
proton particle
static const ChargedStable kaon
charged kaon particle
static const ParticleType photon
photon particle
static const ChargedStable electron
electron particle
Class to store calibrated ECLDigits: ECLCalDigits.
Class to store connected regions (CRs)
TTree * m_tree
tree that contain information for MVA training
double m_cutRatioCorrection
correction for nominator and denominator of the ratio.
double m_signalEnergy[10][5]
total energy per MC matching type of this digit
void getEnteringMother(const MCParticle &particle, int &pdg, int &arrayindex, int &pi0arrayindex)
Get enterging mother of this particle.
StoreArray< ECLConnectedRegion > m_eclConnectedRegions
Store array: ECLConnectedRegion.
float m_timeFitFailed
failed fit
float m_energyRatioNeighbour[c_nMaxNeighbours]
energy ratio of neighbour 0..9 to center
double m_cutOffset
cut offset
double m_totalSignalEnergy
total energy of this digit
float m_phiId
local maximum center theta Id
double m_energyCut
energy cut for seed
virtual void initialize() override
Initialize.
StoreArray< ECLLocalMaximum > m_eclLocalMaximums
Store array: ECLLocalMaximum.
float m_cellId
local maximum center cell Id
virtual void event() override
Event.
virtual const char * eclLocalMaximumArrayName() const
Name to be used for default option: ECLLocalMaximums.
int getIdPosition(const int type, const int id)
Get Id position in the vector.
float m_energy
Variables to monitor the MVA training.
virtual void endRun() override
End run.
std::string m_method
Method to find the local maximum.
virtual void terminate() override
Terminate (close ROOT files here if you have opened any).
ECL::ECLGeometryPar * m_geom
Geometry.
std::string m_outfileName
file name prefix of the training output file
ECL::ECLNeighbours * m_neighbourMap
Neighbour maps.
void makeLocalMaximum(const ECLConnectedRegion &aCR, const int cellId, const int lmId)
Make local maximum for a given connected region.
void resetClassifierVariables()
Reset Classifier Variables.
float m_thetaId
local maximum center theta Id
float m_maxNeighbourEnergy
highest energy of all neighbours
ECLLocalMaximumFinderModule()
Constructor.
bool isEnteringECL(const B2Vector3D &vec)
Check if particle is produced outside of the ECL.
TFile * m_outfile
Output training files and trees.
float m_timeResolution
time resolution
virtual void beginRun() override
Begin.
void addToSignalEnergy(int motherpdg, int motherindex, int pi0index, double weight)
Add energy to vector.
float m_targetpi0index
target array index
virtual const char * eclConnectedRegionArrayName() const
Name to be used for default option: ECLConnectedRegions.
float m_target
MC truth variables.
static const unsigned c_nMaxNeighbours
Variables (possibly) used for MVA classification.
StoreArray< MCParticle > m_mcParticles
Store array: MCParticle.
virtual const char * eclCalDigitArrayName() const
Name to be used for default or PureCsI option: ECLCalDigits.
float m_targetindex
target array index
float m_maxEnergyRatio
Highest energetic neighbour energy divided by LM energy.
int m_signalId[10][5]
total energy per MC matching type of this digit
const double c_minEnergyCut
Minimum LM energy.
std::vector< int > m_StoreArrPosition
vector (ECLElementNumbers::c_NCrystals + 1 entries) with cell id to store array positions
StoreObjPtr< EventLevelClusteringInfo > m_eventLevelClusteringInfo
EventLevelClusteringInfo.
void resetTrainingVariables()
Reset Debug Variables.
void getMax(int &maxtype, int &maxpos)
Get the highest energy deposition particle type.
StoreArray< ECLCalDigit > m_eclCalDigits
Store array: ECLCalDigit.
double m_cutSlope
cut slope.
int m_isTrainingMode
training mode for MVA methods (i.e.
double m_truthFraction
MC truth fraction.
float m_nNeighbours10
Variables (possibly) used for cut classification.
virtual ~ECLLocalMaximumFinderModule()
Destructor.
Class to store local maxima (LM)
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
int GetPhiID()
Get Phi Id.
void Mapping(int cid)
Mapping theta, phi Id.
int GetThetaID()
Get Theta Id.
Class to get the neighbours for a given cell id.
const std::vector< short int > & getNeighbours(short int cid) const
Return the neighbours for a given cell ID.
A Class to store the Monte Carlo particle information.
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
static const double MeV
[megaelectronvolt]
static const double cm
Standard units with the value = 1.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
bool isForward(int cellId)
Check whether the crystal is in forward ECL.
bool isBarrel(int cellId)
Check whether the crystal is in barrel ECL.
const int c_NCrystals
Number of crystals.
bool isBackward(int cellId)
Check whether the crystal is in backward ECL.
Abstract base class for different kinds of events.