dipoleHiRes.h

#ifndef _DIPOLE_H_
#define _DIPOLE_H_

#include "healpixmap.h"
#include "events.h"
#include "TH1F.h"
#include "TF1.h"

#define WITHSINE

/*!
  The method for detecting the presence of a dipole source model will be based upon comparison between the real 
  data and a large quantity of events generated by Monte Carlo simulation program. In order to measure the 
  presence of a dipole in the events sample, a binning technique that considers the event counts for the full 
  range of opening angles from the center of each proposed dipole distribution is used. The method is described 
  in : HiRes Collaboration, Astroparticle Physics 21 (2004) 111-123
*/

//! Look for a dipolar modulation in the data set.
class TDipoleHiRes
{
 public:
  //! Constructor.
  TDipoleHiRes(const vector<TEvent>& events, double lDipole, double bDipole, double dCosTheta = 0.05);
  
  //! Destructor.
  ~TDipoleHiRes();
  
  //! Set fExtension. 
  void SetExtension(string ext) {fExtension = ext;}

  /*! 
    Computes the cosinus opening angle distribution. For each event \f$\vec{u_i}\f$, the cosinus opening angle 
    \f$ \cos\theta = \vec{u_i}.\vec{u}_{dipole} \f$ is computed. The #fHistoEvents histogram of #fNbins. 
   */ 
  void ComputeEventsOpeningAngle();

  /*!
    Computes the background cosinus opening angle distribution using Monte Carlo simulation. The simulated data 
    possess the zenith angle given by thVal and pthVal. The #fHistoCov histogram of #fNbins is then filled.
   */
  void ComputeCoverageOpeningAngle(double longitude, double latitude, const THealpixMap& map, unsigned int nSimu, const vector<double>& thVal, const vector<double>& pthVal, string utcFile = "", string jdFile = "", string globalFile = "");

  //! The true cosinus opening angle is obtained from the ratio of #fHistoEvents and #fHistoCov. 
  void ComputeTrueOpeningAngle();

  //! Linear fit of the normalized bin count.
  void fitDipole();

  //! Draw #fHistoEvents.
  void DrawEvents() const;
  
  //! Draw #fHistoCov.
  void DrawCov() const;

 //! Draw #fHistoEvents and #fHistoCov in the same canvas.
  void DrawBoth() const;

  //! Draw #fHistoTrue.
  void DrawTrue() const;

  //! Draw #fHistoTrue and #fFitFcn in the same canvas.
  void DrawTrueFit() const;

  //! Returns #fNbins.
  int GetNbins() const {return fNbins;}

  //! Returns #fFitParameters.
  vector<double> GetFitParameters() const {return fFitParameters;} 

  //! Returns #fFitParametersErrors.
  vector<double> GetFitParametersErrors() const {return fFitParametersErrors;}

  //! Returns the number of counts in each bin of #fHistoEvents.
  vector<double> GetEventsCounts() const;

  //! Returns the number of counts in each bin of #fHistoCov.
  vector<double> GetCoverageCounts() const;

  //! Returns the normalized count in each bin of #fHistoTrue.
  vector<double> GetTrueValues() const;

  //! Returns #fCosThetaBins
  vector<double> GetBinsCenter() const {return fCosThetaBins;}

private:
  //! Computes #fCosTheta.
  void ComputeOpeningAngle(const vector<TEvent>& events);

  //! Initializes #fNbins, #fCosThetaBins, #fHistoEvents, #fHistoCov and #fHistoTrue.
  void InitHisto();

  //! Events (real or simulated).
  const vector<TEvent>& fEvents;

  //! Number of bins of #fHistoEvents, #fHistoCov and #fHistoTrue.
  unsigned int fNbins;

  //! Size of the bins of #fHistoEvents, #fHistoCov and #fHistoTrue.
  double fdCosTheta;

  //! Opening angle between the arrival direction of an event and the center of the proposed dipole source model.
  double * fCosTheta;

  //! Events cosinus opening angle distribution.
  TH1F * fHistoEvents;

  //! Background cosinus opening angle distribution.
  TH1F * fHistoCov;

  //! True correlation.
  TH1F * fHistoTrue;

  //! Bins center of #fHistoEvents, #fHistoCov and # fHistoTrue.
  vector<double> fCosThetaBins;

  //! Number of events.
  unsigned int fNevents;

  //! Galactic longitude of the dipole.
  double fLdipole;

  //! Galactic Latitude of the dipole.
  double fBdipole;

  //! Function used for the fit.
  TF1* fFitFcn;

  //! Parameters of the fit.
  vector<double> fFitParameters;
  
  //! Error on the parameters of the fit.
  vector<double> fFitParametersErrors;

  //! Figure extension.
  string fExtension;
};

//! Linear function.
double linearFunction(double *x, double *par);

#endif