MdLDFFinder.h

Go to the documentation of this file.

#ifndef _MdLDFFinderAG_MdLDFFinder_h_
#define _MdLDFFinderAG_MdLDFFinder_h_

#include <utl/config.h> // before all
#include <utl/Point.h>
#include <utl/Vector.h>
#include <fwk/VModule.h>
//
#include <evt/ShowerRecData.h>
#include <mevt/MEvent.h>
#include <mevt/ModuleRecData.h>
#include <mdet/MDetector.h>
//
#include <boost/utility.hpp>
#include <utility>
//
#include "VMinMethodFunctor.h"

namespace mdet {
   class Counter;
}

namespace evt {
   class ShowerMRecData;
}

namespace MdLDFFinderAG {

  class VLDFFunctor;
  class VMinMethodFunctor;

  class MdLDFFinder : public boost::noncopyable,
                      public fwk::VModule {
  public:
    enum LDFType {
      eKascadeGrande, 
      eUndef_I,       
      eUndef_II       
    };
    static const char* const kLDFTags[];

    enum MinMethodType {
      eLikelihood, 
      eLikelihood2,       
      eLikelihood3       
    };
    static const char* const kMinMethodTags[];

    MdLDFFinder();
    ~MdLDFFinder();
    VModule::ResultFlag Init();
    VModule::ResultFlag Run(evt::Event& e);
    VModule::ResultFlag Finish();



  private:

   /*
   *   \brief The MD reconstructor state
   *   This class was designed to keep the same state after
   *   processing each event. This implies that its members
   *   do not refer to invidual events but to the
   *   module configuration. This design prevents passing
   *   unwanted data between events.
   *
   */

    CounterList SelectCandidateCounters(mevt::MEvent&);
    CounterList SelectSilentCounters(mevt::MEvent&);

    double CalculateRho0Seed( const CounterList counters, const utl::Point rCore, const utl::Vector rAxis) const;
    double CalculateBeta( const double theta) const;

// From Pierre Billoir Thoughts....(same as SdReconstruction)
// - if there are enough stations (especially at distances around ReferenceDistance), a fit
// with a variable LDF slope (beta) is attempted; the conditions are:
// * at least fCountersToFixBeta counters, with one of the following requirements:
// * at least 2 within (fLowLimToRelaxBeta, fUppLimToRelaxBeta), with a difference at least fCounterSpacings1 in between
// * at least 3 within (fLowLimToRelaxBeta, fUppLimToRelaxBeta), with a maximal difference at fCounterSpacings2
// * at least 4 within (fLowLimToRelaxBeta, fUppLimToRelaxBeta), with a maximal difference at fCounterSpacings3

    bool FixBetaFlag( const CounterList, const utl::Point&, const utl::Vector& ) const;

/*
* The fit residuals
* The residuals of the muon LDF fit are calculated for each candidate module that is not saturated.
* The fit residuals are the measured density minus the density according the LDF. All densities
* are referred to the shower plane.
*/
    void SetLdfResiduals(CounterList counters, const std::vector<double> parameters,
                         const utl::Point rCore, const utl::Vector rAxis) const;

/*
The fit chi2
The chi2 of the fit is calculated from candidates modules that are not saturated.
This chi2 should not be used for any serious analysis as it is just a raw approximation. For a better
estimator of the goodness of the fit use instead the value of the likelihood minimum.
For each module the chi2 is calculated from the number of measured muons and the number of muons expected
from the LDF. A Poisson variance equal to the number of muons from the LDF is assumed. However this approximation
fails for the limits of few and many muons. In the case of few muons the approximation is invalid because Poisson
cannot be approximated by a Gaussian. For many muons the problem is that the variance of the segmented counter is
greater than the Poisson variance (see Appendix A of GAP-2014-37).
*/
    double CalculateChi2(const CounterList counters, const std::vector<double> parameters,
                         const utl::Point rCore, const utl::Vector rAxis) const;

    size_t CalculateDegreesOfFreedom(const CounterList, const std::vector<bool>) const;

    void FillTabulatedFunction(evt::ShowerMRecData&, const std::vector<double>& fPars,
                               const std::vector<double>& fErrPars);

    // Fill reconstructed shower with parameters
    void FillReconstructedParameters( evt::ShowerMRecData&, const ROOT::Minuit2::MnUserParameterState& , double );

    // Helper method to get the index of a free fit parameter skipping the fixed parameters
    // It is required to select elements in the covariance matrix provided by Minuit
    unsigned int GetFreeParIndex(const std::vector<ROOT::Minuit2::MinuitParameter>, unsigned int) const;

    // Filling the SP distances of the mevt::ModuleRecData objects
    void FillModulesShowerPlaneDistances(const evt::ShowerSRecData&, const mdet::MDetector&, mevt::MEvent&);

    bool fFixBeta;  // Beta flag setting from the configuration. Input to all events
    unsigned int fCountersToFixBeta;


    double fLowLimToRelaxBeta;
    double fUppLimToRelaxBeta;

    double fCounterSpacings1;
    double fCounterSpacings2;
    double fCounterSpacings3;
    double fA0;
    double fA1;
    double fA2;

    double fReferenceDistance;
    bool fUseSilent;
    size_t fSilentLimit;

    bool fUseSaturated;
    size_t fSaturationLimit;

    bool fFixCoreFromSd;

    LDFType fLDFType;
    MinMethodType fMinMethodType;
    std::unique_ptr<VLDFFunctor> fLDF;
    std::unique_ptr<VMinMethodFunctor> fMinimizer;

    utl::CoordinateSystemPtr siteCS;
                const mdet::MDetector* mDetector;

    REGISTER_MODULE("MdLDFFinderAG", MdLDFFinder);
  };
}

#endif // _MdLDFFinderAG_MdLDFFinder_h_