#ifndef LCDCALRECON_H
#define LCDCALRECON_H

#include "TMath.h"
#include "TRandom.h"
#include "LCDRecModule.h"  
#include "LCDGetParameters.h"
#include "LCDMcPart.h"
#include "LCDDetectorVolume.h"
#include "LCDSmearTraj.h"
#include "LCDCluster.h"

//=========================================================
class LCDCalRecon : public LCDRecModule {
 private:
  //Int_t       nPart;         // Total number of MC particles in event
  //TObjArray   smears;        // "Clusters" -- actually of class CalSmeared 
  TObjArray*  ClsList;       // "Clusters" 
  // Ultimately should learn how to use TCloneArray for smears
  
  static Double_t RanExpDecay(TRandom *pRand, Double_t mean = 1.0)
    {return -TMath::Log(pRand->Rndm()) * mean;}
  
  Int_t   Swim(LCDSmearTraj* pTraj, Bool_t hadron);
  
  LCDGetParameters*  gp;  
  
  Int_t  Smearable(LCDMcPart*);
  
  TRandom m_rand;  // Use Root for randoms, gaussian or uniform
  
  // Assume all relevant volumes are either barrels or endcaps, in any
  // case cylinders
  
  LCDDetectorVolume* pVolumes;       // Array of detector volumes
  LCDDetectorVolume* pInnerVolume;   // Innermost volume about IP
  Int_t        nVolumes;       // Number of detector volumes in array

  Double_t  merge_siz_EM   ; // Cut Value for merging. by cm (0.8 cm).
  Double_t  merge_siz_HAD  ; // Default value = 1.3 cm.
  Double_t  merge_angle_EM ; // Cut Value for merging. by radian (0.0200 rad)
  Double_t  merge_angle_HAD; // Default value = 0.0600 rad
  
  Int_t m_cluster_merge; // =1 With cluster merging (default) 
                         // =0 Without cluster merging 

  void Init();

 public:
  LCDCalRecon();
  LCDCalRecon(LCDGetParameters* getparam);
  ~LCDCalRecon(){}
  
  void cleanup();
  void doit(LCDEvent* event);
  void spew(FILE* ofile)const;
  
  void SetDetectorParameters(LCDGetParameters* getparam);
  void Smear(Double_t Epart,LCDSmearTraj& traj, Bool_t hadron, 
	     LCDGetParameters* gp, LCDCluster* kal);
  void doMerge();
  void AddCluster(LCDCluster* kal1, LCDCluster* kal2);
  
  void SetMergeSizeEM(Double_t a)   { merge_siz_EM    = a; }
  void SetMergeSizeHAD(Double_t a)  { merge_siz_HAD   = a; }
  void SetMergeAngleEM(Double_t a)  { merge_angle_EM  = a; }
  void SetMergeAngleHAD(Double_t a) { merge_angle_HAD = a; }
  
  void SetNoClusterMerge() {m_cluster_merge = 0;}

  TVector3 ToPolar(const TVector3& pRect);
  TVector3 Rotate(const TVector3& inVec, const TVector3& angleVec);
  TVector3 RotateTheta(const Double_t theta, const TVector3& iVec);
  TVector3 RotatePhi(const Double_t phi, const TVector3& iVec);
  
  ClassDef(LCDCalRecon,0)    // Manager class for FastMC clustering
};

// energy smearing is done according to the formula
//   sigma = E-perfect * sqrt(a**2/E-perfect + b**2)
// Magnitude of transverse position smear is
//   sigma_x = sqrt(a**2/E-perfect + b**2)
//   phi (with axis = momentum angle) is uniformly distributed in (0, 2pi)
#endif