ECLSplitterN1Module.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2016 - Belle II Collaboration                             *
 *                                                                        *
 * Main reconstruction splitter code for the nPhoton hypothesis.          *
 * Based on a connected region (CR) we look for local maxima and          *
 * create one shower for each local maximum (LM). In case of multiple     *
 * LM in one CR the energy is shared between the showers based on         *
 * their exponentially weighted distance in an iterative procedure. If    *
 * a CR has no LM the highest energetic digit in the CR is taken as LM.   *
 * The position is reconstructed using logarithmic weights for not too    *
 * small shower and linear weights otherwise ('lilo').                    *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Torben Ferber (ferber@physics.ubc.ca) (TF)               *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#pragma once
 
//STL
#include <vector>
 
// FRAMEWORK
#include <framework/core/Module.h>
#include <framework/gearbox/Unit.h>
#include <framework/datastore/StoreArray.h>
#include <framework/datastore/StoreObjPtr.h>
 
class TGraph2D;
class TFile;
class TH1F;
 
namespace Belle2 {
  class ECLCalDigit;
  class ECLConnectedRegion;
  class ECLShower;
  class ECLLocalMaximum;
  class EventLevelClusteringInfo;
 
  namespace ECL {
    class ECLNeighbours;
    class ECLGeometryPar;
  }
 
  class ECLSplitterN1Module : public Module {
 
  public:
    ECLSplitterN1Module();
 
    ~ECLSplitterN1Module();
 
    virtual void initialize() override;
 
    virtual void beginRun() override;
 
    virtual void event() override;
 
    virtual void endRun() override;
 
    virtual void terminate() override;
 
  private:
    // Module parameters:
 
    // Splitter
    double m_threshold; 
    double m_expConstant; 
    int m_maxIterations; 
    double m_shiftTolerance; 
    double m_minimumSharedEnergy; 
    int m_maxSplits; 
    const double c_molierRadius = 3.581 *
                                  Belle2::Unit::cm; 
    double m_cutDigitEnergyForEnergy; 
    double m_cutDigitTimeResidualForEnergy; 
    int m_useOptimalNumberOfDigitsForEnergy; 
    std::string m_fileBackgroundNormName; 
    std::string m_fileNOptimalFWDName; 
    std::string m_fileNOptimalBarrelName; 
    std::string m_fileNOptimalBWDName; 
    TFile* m_fileBackgroundNorm{nullptr}; 
    TFile* m_fileNOptimalFWD{nullptr}; 
    TFile* m_fileNOptimalBarrel{nullptr}; 
    TFile* m_fileNOptimalBWD{nullptr}; 
    TH1F* m_th1fBackgroundNorm{nullptr}; 
    TGraph2D* m_tg2dNOptimalFWD[13][9] {}; 
    TGraph2D* m_tg2dNOptimalBWD[10][9] {}; 
    TGraph2D* m_tg2dNOptimalBarrel{nullptr}; 
    // Position
    std::string m_positionMethod;  
    double m_liloParameterA; 
    double m_liloParameterB; 
    double m_liloParameterC; 
    std::vector<double> m_liloParameters; 
    // Background
    int m_fullBkgdCount; 
    const unsigned short c_nSectorCellIdBWD[10] = {9, 9, 6, 6, 6, 6, 6, 4, 4, 4}; 
    const unsigned short c_nSectorCellIdFWD[13] = {3, 3, 4, 4, 4, 6, 6, 6, 6, 6, 6, 9, 9}; 
    // Crystals per Ring
    const unsigned short c_crystalsPerRing[69] = {48, 48, 64, 64, 64, 96, 96, 96, 96, 96, 96, 144, 144, //FWD (13)
                                                  144, 144, 144, 144, 144, 144, 144,  // BARREL 20
                                                  144, 144, 144, 144, 144, 144, 144, 144, 144, 144, //30
                                                  144, 144, 144, 144, 144, 144, 144, 144, 144, 144, //40
                                                  144, 144, 144, 144, 144, 144, 144, 144, 144, 144, //50
                                                  144, 144, 144, 144, 144, 144, 144, 144, 144,//59
                                                  144, 144, 96, 96, 96, 96, 96, 64, 64, 64
                                                 }; 
    std::vector< int > m_StoreArrPosition;
 
    std::vector< int > m_StoreArrPositionLM;
 
    std::vector< int > m_cellIdInCR;
 
    ECL::ECLNeighbours* m_NeighbourMap9{nullptr}; 
    ECL::ECLNeighbours* m_NeighbourMap21{nullptr}; 
    StoreArray<ECLCalDigit> m_eclCalDigits;
 
    StoreArray<ECLConnectedRegion> m_eclConnectedRegions;
 
    StoreArray<ECLShower> m_eclShowers;
 
    StoreArray<ECLLocalMaximum> m_eclLocalMaximums;
 
    StoreObjPtr<EventLevelClusteringInfo> m_eventLevelClusteringInfo;
 
    virtual const char* eclCalDigitArrayName() const
    { return "ECLCalDigits" ; }
 
    virtual const char* eclConnectedRegionArrayName() const
    { return "ECLConnectedRegions" ; }
 
    virtual const char* eclLocalMaximumArrayName() const
    { return "ECLLocalMaximums" ; }
 
    virtual const char* eclShowerArrayName() const
    { return "ECLShowers" ; }
 
    virtual const char* eventLevelClusteringInfoName() const
    { return "EventLevelClusteringInfo" ; }
 
    ECL::ECLGeometryPar* m_geom{nullptr};
 
    void splitConnectedRegion(ECLConnectedRegion& aCR);
 
    int getNeighbourMap(const double energy, const double background);
 
    unsigned int getOptimalNumberOfDigits(const int cellid, const double energy, const double background);
 
    double getEnergySum(std::vector < std::pair<double, double> >& weighteddigits, const unsigned int n);
 
    double estimateEnergy(const int centerid);
 
  }; // end of ECLSplitterN1Module
 
 
  class ECLSplitterN1PureCsIModule : public ECLSplitterN1Module {
  public:
    virtual const char* eclCalDigitArrayName() const override
    { return "ECLCalDigitsPureCsI" ; }
 
    virtual const char* eclConnectedRegionArrayName() const override
    { return "ECLConnectedRegionsPureCsI" ; }
 
    virtual const char* eclLocalMaximumArrayName() const override
    { return "ECLLocalMaximumsPureCsI" ; }
 
    virtual const char* eclShowerArrayName() const override
    { return "ECLShowersPureCsI" ; }
 
    virtual const char* eventLevelClusteringInfoName() const override
    { return "EventLevelClusteringInfoPureCsI" ; }
 
  }; // end of ECLSplitterN1PureCsIModule
 
} // end of Belle2 namespace