RdWaveFit.h

Go to the documentation of this file.

#ifndef _RdWaveFit_h_
#define _RdWaveFit_h_

#include <fwk/VModule.h>
#include <revt/Channel.h>
#include <revt/REvent.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/TimeStamp.h>
#include <utl/CoordinateSystem.h>
#include <utl/GeometryUtilities.h>
#include "Chi2ForPlaneWaveFit.h"
#include "Chi2ForSphericalWaveFit.h"
#include "Chi2ForSphericalWaveFitVarC.h"
#include "Chi2ForConicalWaveFit.h"

/*namespace evt {
 class Event;
 }*/

namespace revt {
  class REvent;
}

// you should assign your module to a namespace

namespace RdWaveFit {

  struct FitParameters {

    // !!! all angles in rad !!!
    double theta;
    double thetaError;
    double phi;
    double phiError;
    double r;   // Radius of a spherical wave fit
    double rError;
    double rho;  // angle for conical wave fit
    double rhoError;

    double c00;
    double c11;
    double c22;
    double c01;
    double c02;
    double c12;

    double gamma;       // gamma = c/c_0
    double gammaError;

    double chiSquareNDF;
    int status;
    int nCalls;
    double minChi;
    int NDF;
    bool corrected;
    int nStations;
    double recStage;    // RecStage defined in trunk/Framework/Shower/ShowerRRecData.cc

    // type defines which minimization mode is used: 0=migrad, 1=simplex, 2=combined, 3=scan
    // see Minuit2 Documentation for further information
    int type;
    FitParameters() :
          theta(0),
            thetaError(0),
            phi(0),
            phiError(0),
            r(0),
            rError(0),
            rho(0),
            rhoError(0),
            c00(0),
            c11(0),
            c22(0),
            c01(0),
            c02(0),
            c12(0),
            gamma(0),
            gammaError(0),
            chiSquareNDF(0),
            status(0),
            nCalls(0),
            minChi(0),
            NDF(0),
            corrected(0),
            nStations(0),
            recStage(0),
            type(0)

    {
    }
  };

  class RdWaveFit: public fwk::VModule {

  public:
    // If you define a constructor, you MUST provide an implementation
    // for it, otherwise the REGISTER_MODULE macro will not register
    // your module. (Here, implementation is in RdWaveFit.cc)
     RdWaveFit();
    //virtual ~RdWaveFit();

    fwk::VModule::ResultFlag Init();
    fwk::VModule::ResultFlag Run(evt::Event& event);
    fwk::VModule::ResultFlag Finish();

  private:
    fwk::VModule::ResultFlag PlaneWaveFit(revt::REvent& rEvent, FitParameters& fit, int type) const;

    fwk::VModule::ResultFlag SphericalWaveFit(revt::REvent& rEvent, FitParameters& fit, int type) const;

    fwk::VModule::ResultFlag SphericalWaveFitVarC(revt::REvent& rEvent, FitParameters& fit, int type) const;

    fwk::VModule::ResultFlag ConicalWaveFit(revt::REvent& rEvent, FitParameters& fit, int type) const;

    fwk::VModule::ResultFlag ScanPWF(revt::REvent& rEvent, FitParameters& sFitResults) const;

    fwk::VModule::ResultFlag ScanSWF(revt::REvent& rEvent, FitParameters& sFitResults) const;

    fwk::VModule::ResultFlag ScanSWFVarC(revt::REvent& rEvent, FitParameters& sFitResults) const;

    fwk::VModule::ResultFlag ScanCWF(revt::REvent& rEvent, FitParameters& sFitResults) const;

    fwk::VModule::ResultFlag CallPlaneWaveFit(revt::REvent& rEvent, FitParameters& sFitResults,
        bool reuseFit = false) const;

    fwk::VModule::ResultFlag CallSphericalWaveFit(revt::REvent& rEvent, FitParameters& sFitResults,
        bool reuseFit = false);

    fwk::VModule::ResultFlag CallSphericalWaveFitVarC(revt::REvent& rEvent, FitParameters& sFitResults, bool reuseFit =
        false);

    // ############################ Conical Wave Fit ###############################
    fwk::VModule::ResultFlag CallConicalWaveFit(revt::REvent& rEvent, FitParameters& sFitResults, bool reuseFit = false);
    // ########################### Conical Wave Fit End ############################

    fwk::VModule::ResultFlag ComputeBaryCenter(const revt::REvent& rEvent, utl::Point& baryCenter);

    const std::vector<double> GetTimeResiduals(const std::vector<utl::Vector> &AntennaPositions,
        const std::vector<double> &AntennaTimes, const utl::Vector& ShowerAxis,
        const unsigned int waveFrontModel, const double gamma = 1, const double rho = 0) const;

    utl::CoordinateSystemPtr Create(utl::Point arg1);

    enum InfoLevel {
      eNone = 0, eFinal = 1, eIntermediate = 2, eObscure = 3, eMinuit = 4
    };

    enum WaveFrontModel {
      ePlane = 0, eSpherical = 1, eSphericalVarC = 2, eConical = 3
    };

    // xml settings
    int fInfoLevel;
    bool fMinuitOutput;
    int fMinNumberOfStations;
    bool fAllowUnphysicalCosines;
    int fFitType;
    double fThetaMax;
    double fRMax;
    double fRhoMax;  // for conical fit
    int fNInitsTheta;
    int fNInitsPhi;
    bool fSphericalWaveFit;
    unsigned int fWaveFrontModel;
    unsigned int currentWaveFrontModel; // adjusts WaveFrontModel if number of stations is not sufficient

    //VarC Specific
    bool fSphericalWaveFitVarC;
    bool fComputeSigmaGamma;
    double fSigmaGamma;
    double fPenaltyTerm;
    double fThetaVarC;

    int fNumberOfStations;

    double fEventTimeNanoSecond;
    double fEventTimeSecond;
    // minuit globals
    utl::Point fReferenceCore;
    utl::CoordinateSystemPtr fLocalCS;

    Chi2ForPlaneWaveFit* ChiPWF;
    Chi2ForSphericalWaveFit* ChiSWF;
    Chi2ForSphericalWaveFitVarC* ChiSWFVarC;
    Chi2ForConicalWaveFit* ChiCWF;

    std::string fUsedDirection;
    std::string fUsedCorePosition;

    // The REGISTER_MODULE macro makes the framework aware of your module.
    // You need to put this macro at the very end of your module class definition.
    // Give your module a name (first argument) and specify its 
    // class name (second argument).  In the ModuleSequenceExample.xml file,
    // you can see this module being refered to by the name we give it here.

  REGISTER_MODULE("RdWaveFit", RdWaveFit);

  };

}

#endif

// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -k"
// End: