//==================================================================
//
// TRphiWt.hh
//
// Finds the optimal weight for rphi hits in L2 using a simple
// straight line fit for L1 & L3, intersecting with L2 and fitting for
// a unit gaussian for the distribution of (pred-L2)/sigma.  Where
// sigma is a simple function of the weight.
//
// Aaron Dominguez (LBNL) May 17, 2002
//
//==================================================================
#if !defined (__CINT__) || defined (__MAKECINT__)

#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TGraphErrors.h"
#include "TClonesArray.h"
#include "TMinuit.h"

#include "Stntuple/loop/TStnModule.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"
#include "Stntuple/obj/TStnTrackLinkBlock.hh"
#include "Stntuple/obj/TSvxDataBlock.hh"
#include "Stntuple/obj/TSiIsectBlock.hh"
#include "Stntuple/obj/TSiStripBlock.hh"
#include <vector>
#include <string>
#include "Stntuple/obj/TSiAlign.hh"

#endif 

class TRPhiWt: public TStnModule {

public: 
  TRPhiWt(const char* name="RPhiWt", const char* title="Silicon Weights for RPhi", const char* filename="");
  ~TRPhiWt();

  static TRPhiWt *gRPhiWt;
  
  TTree *fT,*fC; // Tree of hits and a copy of tree after some cuts.  The copy fC is fit.
  struct rPhi_t {
    int w;
    double pt;
    double r,er; // radius of curvature from COT+SVX track fit and it's error
    double x[3]; // Global x-y coord rotated to vertical wedge 
    double y[3];
    int n[3]; // Cluster type
    int nst[3]; // Number of strips in cluster
    float q[3]; // Cluster charge
    double dx[3]; // Change in cluster centroid due to reclustering
  };
  TH1F *fR,*fR_full; // residual histogram
  TH1F *fS,*fS_full; // residual/fitted_error histogram
  TH2F *fPR, *fPS; // Residual and pull histogram as function of impact angle
  TH2F *fWR, *fWS; // Residual and pull histogram as function of wedge
  TH1F *fD; // Diff between perfect circle and approx at L2
  TH1F *fD0; // Diff between perfect circle and straight line at L2
  TGraphErrors *fG[4],*fW;

  // ****** Setters *******
  void SetNStripPerClus(int nType, const int *nmin, const int *nmax);
  void SetAnchor(int iType){fAnchorType=iType;}
  void LoadFile(bool load=true) {fLoadFile=load;}
  void PtCut(double minpt) {fMinPt=minpt;}

  void UseAlignment(const std::string &geoFile="./Stntuple/db/si/siGeo_100030_1_GOOD.txt");

  void FudgeClusters(float maxQ=127.0) {fMaxQ=maxQ;}
  void RecalculateClusterCentroid(int ihit, double &dx, float &q);

  // ****** overloaded methods of TStnModule ******
  int       BeginJob();
  int       BeginRun();
  int       Event   (int ientry);
  int       EndJob  ();
  
  // ****** other methods ******
  void BookHistograms();
  bool PassesTrackCuts(int itrk, int flag=0);
  int  ClusterType(int nstr);
  void GetSigmas(TMinuit *m, double *sigma, double *error=NULL);
  void AssignSigmas(const double *sigs, const int *n, double *layer_sigma);
  void FitWt();

  static void Likelihood(Int_t &npar, Double_t *gin, Double_t &f, 
		  Double_t *par, Int_t iflag);

  double a(const double *x, const double *y) const;
  double b(const double *x, const double *y) const;
  double sigma_a2(const double *x, const double *y, double s) const;
  double sigma_b2(const double *x, const double *y, double s) const;
  double cov_ab(const double *x, const double *y, double s) const;

  // Residual from straight line
  double resid(const double *x, const double *y) const;
  double sigma_resid2(const double *x, const double *y, double s) const;

  // Residual from circle
  double residC(const double *x, const double *y, double r) const;

  // Residual from linearization of circle
  double residL(const double *x, const double *y, double r) const;
  // Old error calculation was slightly incorrect and assumed equal errors for all layers
  double sigma_residL2_old(const double *x, const double *y, double r, double s) const;
  double sigma_residL2(const double *x, const double *y, double r, const double *s) const;

  // Helpers for linearization
  double A(const double *x, const double *y, double r) const;
  double sL(const double *x, const double *y, double r, double xL) const;

  double logGaussianZeroMean(double r, double sigma_r2);
  double logGaussian(double r, double sigma_r2, double mean);

  void circleFrom2PointsAndRadius(const double *x, const double *y, //input points
				  double r, // input radius
				  double &xc, double &yc) const; // output center of circle
  int CutTree(double mean, const double *sigma, double minpt) ; // Applies sigma cuts on tree
  // ****** private data ******
private: 
  TStnTrackBlock     *fTr;
  TStnTrackLinkBlock *fTl;
  TSvxDataBlock      *fHit;
  TSiIsectBlock      *fI;
  TSiStripBlock      *fSi;

  std::string fFilename;
  TFile *fFile;

  rPhi_t fp,fc;

  int fNStrTyp,*fNStrMin,*fNStrMax;
  bool fLoadFile;
  double fMinPt;
  int fAnchorType;

  TSiAlign *fGeo;

  float fMaxQ;

  ClassDef(TRPhiWt,0)
};