UnivRec.h

Go to the documentation of this file.

#ifndef __UnivRecNS_h_
#define __UnivRecNS_h_

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <fstream>

#define ExternalUse // Used within Offline

#ifdef ExternalUse
#include <tls/Atmosphere.h>
#include <tls/UnivParam.h>
#include <tls/UnivParamTime.h>
#else
#include <ToyShower.h>
#endif

#include <TVector3.h>
#include <tls/UnivCalibConstants.h>

#ifndef UNDEF
#define UNDEF -100
#endif

namespace UnivRecNS {

  const int debug=0;
 
  //--- Inherited from DetectorResponse.h ( change accordingly) 

  //-----------------------------
  // Detector geometry 
  //-----------------------------

  // WCD
  const double tankh = 120.;
  const double tankr = 180.;  

  // ASCII 2 m^2
  const double scinlength= 180.;
  const double scinwidth = 108.;
  const double scinheight= 1.;

  // AMIGA:  64 strips in 2 sides with 5.12 m^2 each 
  const double nMD=2.;
  const double MDlength= 400.;
  const double MDwidth = 128.;
  const double MDheight= 1.;  
  
 
  //-----------------------------
  //---- MC specific settings
  //-----------------------------
  const int nAoP=4;
  const double AreaOverPeak_Run[nAoP]={3.60,3.60,4.30,2.93};
  const double SmearingAngle=0.25; // Smearing angle for MC (Hybrid angular resolution)  

  //Saturation ( vem signal in a FADC bin )
  const double SaturationLevelVEM_WCD[2]={200.,2000.}; // [0] Current SDE [1] Upgraded SDE
  const double SaturationLevelVEM_Scin[2]={200.,800.}; // [0] Current SDE [1] Upgraded SDE
  //Saturation ( number of muons in 24 [ns] in a module with 64 scintillator bars )
  const double SaturationLevelMD=64;        

  //-----------------------------
  //----- Likelihood related
  //-----------------------------

  const int nQ=4;
  const double lnP_inf=-200.;

  const int nDetectorType=3;

  const unsigned int npar=9;
  const double unit[npar]={ 20.e2, 20.e2, 0.03, 0.20, 30.,  100., 2.e-3, 2.e-3 , 0.01 };
  const double low[npar]= {  -1500.e2,  -1500.e2, 17.0, 0.05,  400. , -500., -sin(25.*M_PI/180.), -sin(25.*M_PI/180.) , -0.5  };
  const double high[npar]={   1500.e2,   1500.e2, 22.5, 3.5 , 1200. ,  500.,  sin(25.*M_PI/180.),  sin(25.*M_PI/180.) ,  0.5 };
  const std::string ParName[npar]={ "xgcore_relative","ygcore_relative","logE", "Nmu" , "Xmax", "T0", "l", "m", "OffsetM_Mu" }; 
                                          
  const double edm_tolerance=0.1;

  const double InterGPSAccuracy[2]={10.,4.}; 
  const double ParameterizationRMSLimit=5.; 
  const double StartAccuracy[2]={6.25,2.0};  // ASCII/WCD  25/8 [ns] traces
  const double StartAccuracyAMIGA=2.0;       // MC -> 8 [ns] traces     Data -> 3.25 [ns] traces with pattern recognition

  //-----------------------------
  //----- Station selection
  //-----------------------------
  
  const bool UseStartTimeOnlyWhenSaturated=false;
  const double rStartCut[2]= { 0.e2,800.e2};    
  const double vemStartCut=3.;

  const double rQuantileCut[2]= { 0.e2,2000.e2};    
  const double vemTimeCut[nQ]={ 1.e10, 5., 5., 5. };      // [0] T1 [1] T10 [2] T50 [3] T90  ->   WCD [vem]  ASCII [mip]  AMIGA [nmu]  
  
  //---
  const double rSignalCut[2]= {0.e2,2500.e2};   //  WCD  ( ASCII and AMIGA are always used if WCD is flagged as useable )   
  const double vemSignalCut=3.;                          

  
  //-----------------------------
  //----- Event selection
  //-----------------------------

  const int nTimeStationsMin=4;
  const int nSignalStationsMin=3;


  //----------------------------------
  //  Shower Reconstruction results
  //----------------------------------

  typedef struct RecInfo {

    int id;
    double AugerDensity;    
    int iatm;

    double xgcore,xgcore_err;
    double ygcore,ygcore_err;      
    double zgcore;

    double logE, logE_err;

    double Xmax, Xmax_err;

    double Nmu, Nmu_err;

    double theta;
    double azi;

    double T0,T0_err;
    
    double logE_ref,logE_ref_err;
    double theta_ref,azi_ref;
    
    double OffsetM_Mu,OffsetM_Mu_err;

    double RA,Dec;
   
    int FitStatus[2]; 
    float MinEigenVal,edm;
    
    double lnP[3];  //Time/Signal/Constrains

    int nStationsWCD,nStationsScin,nStationsMD;    // non-accidental triggered stations above signal threshold ( Different threshold for Time/Signal rec. )
    int nTimeStations, nSignalStations;         

    int seed [3]; // The ids of the WCD detectors that form the triangle with the highest signal 
    TVector3 theDir_Ref;

    int FitStage; 

    //  Time reconstruction information
    bool AllTimesOK;           // All start time and time quantiles positive
    double TimeRejected[2];    // Rejected station  [0] Maximum standard deviation in time likelihood [1] Id of the station
    double TimeMaxdev[2];      // [0] Maximum standard deviation in time likelihood [1] Id of the station

    // Additional information 
    double rHot;
    int gDetHot;
    int nFCNcalls;
    bool IsSaturated[nDetectorType]; // WCD / ASCII / AMIGA

  } RecInfo_t;  

  //----------------------------------
  //  MC info  
  //----------------------------------
  typedef struct MCInfo{
  
    int im,ip,ie,ith,ish,iatm;

    double xgcore,ygcore,zgcore;
    double theta,azi;
    double logE,Xmax,Nmu,T0;

  } MCInfo_t;

  //----------------------------------
  //  Auger standard reconstruction info  
  //----------------------------------
  typedef struct AugerRecInfo{
  
    double xgcore_SD,ygcore_SD,zgcore_SD;
    double theta_SD,azi_SD;
    double logE_SD,logE_SD_err;
    bool isGoodSD; // (Data) relax T5

    double xgcore_FD,ygcore_FD,zgcore_FD;
    double theta_FD,azi_FD;
    double logE_FD,logE_FD_err;
    double Xmax_FD,Xmax_FD_err;
    bool isGoodFD; // (Data) fiducial cuts

    double RA, dec;
  } AugerRecInfo_t;

  //----------------------------------
  //  Rec Station
  //----------------------------------
  typedef struct RecStation {
    int type; // [0] WCD [1] Scin [2] MD
 
    int id;  

    double xg,yg,zg;
    double GPSnano;
    double AreaOverPeak;
    double DetectorArea; 

    bool mask_quantiles;
    double quantile[nQ],tquantile[2][nQ]; // [Resampled/Avg][iq]
    double tq_rec[nQ],tq_pred[nQ],tq_rms[nQ],tq_corr[nQ];
    bool UseGausPDF_quantiles[nQ];
    double lnP_quantiles;

    bool mask_starttime;
    double starttime;  
    double start_rec,start_pred,start_rms,start_corr;
    bool UseGausPDF_start;
    double lnP_start;

    bool mask_signal;
    double signalsize;
    double signalsize_pred,signalsize_rms;
    double lnP_signal;

    bool isSaturated;

    double r,psi,DX[4],DL[4],T0_CF[4],T0_PF,density;
  } RecStation_t; 



  //----------------------------------
  //  UnivRec 
  //----------------------------------

  class UnivRec  { 
    
  public:
    UnivRec();
    ~UnivRec() {}

    bool InitRec(  bool IsData, int RecType, int CalibOpt, int RecDetector /*all*/, int RecSys /*Only CalibOpt=-2*/ , int RecMixture /*Only CalibOpt=0,1*/ );
    bool InitRec_InternalFramework( bool RecInfill /*all*/, int RecDataSet, int year, std::string &FileName/*Only data*/, bool RecSDEUpgrade, int RecAoP /*Only MC */);
   
    void InitRecStation(RecStation *station);
    void ClearRecEvent();
#ifndef ExternalUse
    bool InitEvent(ToyShowerNS::ToyShower *fevt);
#else
    //#include "UnivRecCDAS.inc"
#endif
    bool InitRecParameters();

    float DistPoints(float x1, float y1, float x2, float y2);
    float NormalizeAngle(const float angle);

    void SetSPCoordinates(int idet);
    void SetAllSPCoordinates();

    bool GetRecEstimate(int itype);
    bool PlaneFit( bool CurvCorr, bool AltitudeCorr, bool RecSeed, int itype );
    bool PlaneFit(std::vector<double> fX,std::vector<double>fY,std::vector<double>fT,std::vector<double>fW);
    void GetRecSeed(int itype);

    void PrintRecInfo( bool PrintResiduals);
    void WriteTextFile(FILE *fp);
    bool RWRecinfo(FILE *fp, bool RWStations, bool readFlag);
    void GetChi2(double *Chi2, int *nChi2);
    void GetResiduals( std::vector<double> &mean, std::vector<double> &rms, std::vector<double> &rec, int optY, std::vector<double> &fx, int optX,bool OnlyNotSaturated );
    void GetResiduals( std::vector<double> &mean, std::vector<double> &rms, std::vector<double> &rec, int optY, std::vector<double> &fx, int optX );
    void GetMeanRMS_AoP(double &mean, double &rms);
    bool ScanRange(int iparX, int iparY, std::vector<double> par, std::vector<int> parStatus,  double *rangeX, double *rangeY, double *d, 
                   double RotAngle, double errDef, double InitStep, double &edm);
    bool ScanRange(int iparX, int iparY, std::vector<double> par, std::vector<int> parStatus,  double *rangeX, double *rangeY, double *d, 
                   double RotAngle, double errDef, double &edm);
    void ScanLikelihood( int iparX, int iparY, int N, std::vector<std::pair<double,double> > &parXY,  std::vector<double> &lnP,  
                         std::pair<double,double> &parXY_min, double &lnP_min);
    bool IsGoodRec();
    bool IsGoodRec(bool CutBadYears);
    bool CheckEnergyCut(bool PrintOut);

    void SetT0FromHot(bool UseT0_relative, double T0_relative );

    double GetTimeLikelihood();
    double GetTimeLikelihood(bool UseAvgQuantiles);
    double GetSignalLikelihood();

    bool Check_DX_DL( RecStation* );
    bool SetFitStage(int FitStage);
    bool StationSelection();
    bool TimeSignalOK();
    bool RejectTimeOutliers();
    void SetOptPDF();
 
    void SaveFCNParameters(std::vector<double> par, std::vector<int> status );
    int  InitFCNParameters(std::vector<double> &par, std::vector<double> &epar , std::vector<int> &status ,std::vector<bool> &hasLimits);

    void Unrotate(double l_p, double m_p, double &theta, double &azi);
    int GetErrors(std::vector<double> par,  std::vector<double> &epar, std::vector<int> parStatus  );
    double GetLikelihood(std::vector<double> par, std::vector<int> status);
      
    //---
    RecStation_t *GetRecStation(int idet) { return &fstation[idet]; }
    RecInfo_t *GetRecInfo()            { return &recinfo; }
    MCInfo_t *GetMCInfo()              { return &MCinfo; }
    AugerRecInfo_t *GetAugerRecInfo()  { return &augerinfo; }

    unsigned int GetGPSsec() { return GPSsec;   }
    unsigned int GetnRecStations() { return nRecStations;}
    bool isUpgrade() { return gRecSDE_Upgrade;}
    bool isMC()      { return gRecMC; }
    void GetXmaxRange(double &Xmax0, double &Xmax1){ Xmax0=XmaxLow; Xmax1=XmaxHigh;}
    void SetXmaxRange(double Xmax0, double Xmax1){XmaxLow=Xmax0; XmaxHigh=Xmax1;}
    bool IsTimeReco(){ return IsTimeRec; }

    //private:
    std::vector<RecStation_t> fstation;
    int nRecStations;
    RecInfo_t   recinfo;
    AugerRecInfo_t augerinfo;
    MCInfo_t    MCinfo;

    int gRecType;
    bool gRecMC,gRecInfill;
    bool gRecSDE_Upgrade;
    int gRecAoP;
    int gCalibOpt,gRecMixture; // Assumption done in the MC reconstruction 
    bool gDetectorTypeMask[nDetectorType];
    int gRecDataSet; //[0] SD standard/infill [1] Golden standard/infill [2] Golden standard/infill (no FOV) [3] SD high energy
                   
    bool IsTimeRec, IsGoldenRec;

    int gRecSys; 

    unsigned int GPSsec;
    double XmaxLow,XmaxHigh;

    UnivCalibConstantsNS::UnivCalibConstants *fcalib;
    UnivParamNS::UnivParam *fsignal[nDetectorType];
    UnivParamTimeNS::UnivParamTime *ftime[nDetectorType];
    AtmosphereNS::Atmosphere *fatm;
  };

}

#endif // __UnivRecNS_h_