RdBeamTimeOptimizer.h

Go to the documentation of this file.

#ifndef _RdBeamTimeOptimizer_h_
#define _RdBeamTimeOptimizer_h_

/* #################################################################

  The class evt/BeamQuantities.h and further support in the framework
  have been deleted with r34469 (29.10.21).

################################################################# */

#include <string>
#include <vector>

#include <fwk/VModule.h>
#include <utl/Vector.h>
#include <utl/Trace.h>
#include <utl/Triple.h>

#include <revt/REvent.h>
#include <evt/ShowerRRecData.h>
//#include <evt/BeamQuantities.h>
#include <revt/Station.h>
#include <revt/StationRecData.h>

#include <Minuit2/FCNBase.h>

#include "RdWaveModel.h"

namespace evt {
  class Event;
}

namespace revt {
  class Station;
}


// you should assign your module to a namespace

namespace RdBeamTimeOptimizer {

  enum gridtype_t { none = 0, polar = 1, carthesian = 2}; 
  
  class RdBeamTimeOptimizer : 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 RdBeamTimeOptimizer.cc)
    RdBeamTimeOptimizer();
    virtual ~RdBeamTimeOptimizer();

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

    utl::TraceD powertrace (revt::REvent& rEvent, evt::ShowerRRecData& rrec, const std::vector<double>& shifts) const;
    utl::TraceD crosscorr (revt::REvent& rEvent, evt::ShowerRRecData& rrec, const std::vector<double>& shifts) const;
    utl::TraceD xtrace (utl::TraceD ccTrace, double ccMean, utl::TraceD pwNormal);

    evt::BeamPeak findPeak (const utl::TraceD& trace);

  protected:
    RdWaveModel* waveModel;
    int startbin, stopbin; //< Configured start and stop time in bins.
    enum gridtype_t gridtype; /*< Configured type of ascii output grid.
                                  Meant for debug, but so far the only output method.*/
    std::string outfile; //< Configured output ascii file name. Prefixed by cc_, pw_ and x_.

    void delaysToShifts(revt::REvent& rEvent, utl::TimeStamp meanTime, std::vector<double> &times);
    void shiftsToDelays(revt::REvent& rEvent, utl::TimeStamp meanTime, std::vector<double> &times);
    
    std::vector<double> delaysFromRRec(revt::REvent& rEvent,  evt::ShowerRRecData& rrec);
    std::vector<double> shiftsFromRRec(revt::REvent& rEvent,  evt::ShowerRRecData& rrec);
    
  private:

    void writeRecData(revt::REvent &rEvent, std::vector<double> &shifts);
    void WriteASCIISpectrum(const revt::Station &,
                            const std::string &) const;
    
    void WriteTrace3D(const utl::TraceD& trace,
                      const std::string& filename,
                      const utl::Triple& coords) const;
    // /// Writes out a trace in an ASCII file with polar CS
    // void WriteTracePolar(const utl::TraceD &,
    //                   const std::string &,
    //                   const utl::Triple &) const;
    // /// Writes out a trace in an ASCII file with carthesian CS
    // void WriteTrace_XYZ(const utl::TraceD &,
    //                  const std::string &,
    //                  const utl::Triple &) const;

    void WriteASCIIPolar (const utl::TraceD &power,
                          const utl::TraceD &cc,
                          const utl::TraceD &x,
                          const evt::BeamPeak &pwPeak,
                          const evt::BeamPeak &ccPeak,
                          const evt::BeamPeak &xPeak,
                          const utl::Triple &coords) const;
    
    void WriteASCII_XYZ (const utl::TraceD &power,
                         const utl::TraceD &cc,
                         const utl::TraceD &x,
                         const evt::BeamPeak &pwPeak,
                         const evt::BeamPeak &ccPeak,
                         const evt::BeamPeak &xPeak,
                         const utl::Triple &coords) const;

    inline void efieldtopower (const revt::StationTimeSeries& in,
                               utl::TraceD& out,
                               int start, int stop) const;

    inline void efieldproduct (const revt::StationTimeSeries& in1,
                               const revt::StationTimeSeries& in2,
                               utl::TraceD& out,
                               int start, int stop) const;

    double peakTime;


    // 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("RdBeamTimeOptimizer", RdBeamTimeOptimizer);

  };

  class powerFCN : public ROOT::Minuit2::FCNBase {
  public:
    powerFCN(revt::REvent *rEventP,
             RdBeamTimeOptimizer *beam):rEventP(rEventP),beam(beam){}
    double operator()(const std::vector<double>& pars) const {
      utl::TraceD trace = beam->powertrace (*rEventP, *rrecP, pars);
      evt::BeamPeak peak = beam->findPeak (trace);
      
      return - (peak.GetHeight()  / peak.GetRMS() );
    }
    double Up() const { return 0.5; }
  protected:
    evt::ShowerRRecData *rrecP;
    revt::REvent *rEventP;
    RdBeamTimeOptimizer *beam;
  };


  class crossFCN : public ROOT::Minuit2::FCNBase {
  public:
    crossFCN(revt::REvent *rEventP,
             RdBeamTimeOptimizer *beam):rEventP(rEventP),beam(beam){}
    double operator()(const std::vector<double>& pars) const {
      utl::TraceD trace = beam->crosscorr (*rEventP, *rrecP, pars);
      evt::BeamPeak peak = beam->findPeak (trace);
      
      return - ( peak.GetHeight() / peak.GetRMS() );
    }
    double Up() const { return 0.5; }
  protected:
    evt::ShowerRRecData *rrecP;
    revt::REvent *rEventP;
    RdBeamTimeOptimizer *beam;
  };
}


#endif

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