/** CbmMuchDigitize.h
 *@author Mikhail Ryzhinskiy <m.ryzhinskiy@gsi.de>
 *@since 19.03.07
 *@version 1.0
 **
 ** CBM task class for digitizing MUCH
 ** Task level RECO
 ** Produces objects of type CbmMuchDigi out of CbmMuchPoint.
 **/


#ifndef CBMMUCHDIGITIZE_H
#define CBMMUCHDIGITIZE_H 1

#include <set>
#include "TStopwatch.h"
#include "TPolyLine.h"
#include "TMath.h"
#include "TLine.h"
#include "CbmTask.h"
#include "TRandom3.h"

//using std::map;
//using std::pair;

class TClonesArray;
class CbmGeoMuchPar;
class CbmMuchDigiPar;
class CbmMuchStation;
class CbmMuchPoint;
class CbmMuchDigiScheme;

class Point;
class Segment;

class CbmMuchDigitize : public CbmTask
{

 public:

  /** Default constructor **/
  CbmMuchDigitize();


  /** Standard constructor **/
  CbmMuchDigitize(Int_t iVerbose);


  /** Constructor with name **/
  CbmMuchDigitize(const char* name, Int_t iVerbose);


  /** Destructor **/
  virtual ~CbmMuchDigitize();


  /** Execution **/
  virtual void Exec(Option_t* opt);

  /** Sets switch whether to use (=1) avalanche or not (=0).
   *  0 by default.
  **/
  void SetUseAvalanche(Int_t useAvalanche) { fUseAvalanche = useAvalanche; }


 private:

  CbmGeoMuchPar*     fGeoPar;        /** Geometry parameter container **/
  CbmMuchDigiPar*    fDigiPar;       /** Digitization parameter container **/
  CbmMuchDigiScheme* fDigiScheme;    /** Digitization scheme **/
  TClonesArray*      fPoints;        /** Input array of CbmMuchPoint **/
  TClonesArray*      fMCTracks;      /** Input array of MCTrack **/
  TClonesArray*      fDigis;         /** Output array of CbmMuchDigi **/
  TClonesArray*      fDigiMatches;   /** Output array of CbmMuchDigiMatches**/
  Int_t              fNPoints;       /** Total number of points **/ 
  Int_t              fNFailed;       /** Total number of points which digitization has failed **/
  Int_t              fNOutside;      /** Total number of points which was found outside a detector **/
  Int_t              fNMulti;        /** Total number of channels that was hitby several points **/
  Int_t              fNDigis;        /** Total number of created Digis **/
  TStopwatch         fTimer;         /** Timer **/
  TClonesArray*      fIntersPoints;  /** Intersection points array for internal use **/
  Int_t              fUseAvalanche;  /** Switch whether to use (=1) avalanches or not (=0). 0 by default. **/
  TRandom3           fRnd;

  /** Map of active channels (pair detectorId, channel number) 
   ** to index of MuchDigi **/
  std::map<std::pair<Int_t, Int_t>, Int_t> fChannelMap; //!


  /** Get parameter containers **/
  virtual void SetParContainers();

  /** Intialisation **/
  virtual InitStatus Init();

  /** Reinitialisation **/
  virtual InitStatus ReInit();



  /** Reset eventwise counters **/
  void Reset();

  /** Advanced digis production using avalanches **/
  Bool_t ExecAdvanced(CbmMuchStation* station, CbmMuchPoint* point, Int_t iPoint);

  /** Simple digis production without avalanches **/
  Bool_t ExecSimple(CbmMuchStation* station, CbmMuchPoint* point, Int_t iPoint);

  /** Builds a TPolyLine from given rectangle parameters 
   *@param x0     X of the center
   *@param y0     Y of the center
   *@param width  Width
   *@param height Height
   *@param sinrot Sin of rotation angle
   *@param cosrot Cos of rotation angle
   **/
  TPolyLine GetPolygon(Double_t x0, Double_t y0, 
		       Double_t width, Double_t height, 
		       Double_t sinrot, Double_t cosrot);

  /** Creates a polygon for the avalanche and calculates its boundaries**/
  TPolyLine CreateSpread(Double_t xIn, Double_t yIn,
			 Double_t xOut, Double_t yOut, Double_t r, 
			 Double_t& xMin, Double_t& yMin, 
			 Double_t& xMax, Double_t& yMax);

  /** Return the angle between two vectors on a plane
   * The angle is from vector 1 to vector 2, positive anticlockwise
   * The result is between -pi -> pi   
   **/
  Double_t Angle2D(double x1, double y1, double x2, double y2);
  
  /** Determines whether the given point is inside the polygon 
   *@param polygon   A polygon
   *@param p         A 2D point
   **/
  Bool_t InsidePolygon(TPolyLine polygon, Point p);

  /** Area of the polygon 
   *@param pointsMap  Map of polygon vertices sorted by polar angle 
   **/
  Double_t Area(std::map<Double_t, Point*> pointsMap);

  /** Determines whether two segments are intersected
   *@param s1  First segment 
   *@param s2  Second segment
   *@param iX  X of intersection point (if exists)
   *@param iY  Y of intersection point (if exists)
   **/
  Bool_t SegmentsIntersect(Segment s1, Segment s2, Double_t& iX, Double_t& iY);

  /** Sorts given points by polar angle with respect to the given center
   *@param points  Array of 2D points to sort
   *@param cp      Central 2D point 
   **/
  std::map<Double_t, Point*> SortByPhi(TClonesArray *points, Point cp);

  /** Determines whether the two given polygons intersect and calculates the
   *  area of intersection.
   *@param polygon1  Fisrt polygon
   *@param polygon2  Second polygon
   *@param area      Intersection area
   **/
  Bool_t PolygonsIntersect(TPolyLine polygon1, TPolyLine polygon2, 
			   Double_t& area);

  /*****************************************************************
   * Function returns a random number distributed according
   * exponential law, which reproduces the gas gain fluctuation
   *@author Volodia Nikulin
   *@since 14/04/2007
   ****************************************************************/
  inline Int_t GasGain();

  /** Parametrization of sigma for electrons for Lndau distribution 
   *@param logT  Logarithm of electron kinetic energy (energy in [MeV])
   **/
  inline double   e_sigma_n_e(double &logT);

  /** Parametrization of most probable value for electrons for Lndau distribution 
   *@param logT  Logarithm of electron kinetic energy (energy in [MeV])
   **/
  inline double   e_MPV_n_e(double &logT);

  ClassDef(CbmMuchDigitize,1)
};
#endif




#ifndef POINT_H
#define POINT_H 1
/** A 2D point inherited from TObject **/
class Point : public TObject 
{
 public:
  Point() {}
  Point(Double_t x, Double_t y) {
    fX = x; fY = y;
  }
  Double_t GetX() { return fX; }
  Double_t GetY() { return fY; }
  void SetX(Double_t x) { fX = x; }
  void SetY(Double_t y) { fY = y; }  

 private:
  Double_t fX;
  Double_t fY;
  ClassDef(Point,1)


};
#endif


#ifndef SEGMENT_H
#define SEGMENT_H 1
/** Holds a segment object **/
class Segment : public TLine 
{
 public:
  
  Segment(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
    :TLine(x1,y1,x2,y2){
    
    // Get a, b, c
    a = y1 - y2;
    b = x2 - x1;
    c = x1*y2 - x2*y1;
  }
  
  Bool_t IsHorizontal(){
    if(TMath::Abs(fY1 - fY2) <=1e-5) return kTRUE;
    return kFALSE;
  }
  
  Bool_t IsVertical(){
    if(TMath::Abs(fX1 - fX2) <= 1e-5) return kTRUE;
    return kFALSE;
  }
  
  // Checks whether the segment contains the given point
  Bool_t ContainsPoint(Double_t x, Double_t y){
    if(TMath::Abs(a*x + b*y + c) > 1e-5) 
      return kFALSE;
    
    if(IsVertical()){
      if(fY1 < fY2) return (y>=fY1 && y<=fY2) ? kTRUE : kFALSE;
      else return (y>=fY2 && y<=fY1) ? kTRUE : kFALSE;
    }
    else if(IsHorizontal()){
      if(fX1 < fX2) return (x>=fX1 && x<=fX2) ? kTRUE : kFALSE;
      else return (x>=fX2 && x<=fX1) ? kTRUE : kFALSE;
    }
    else{
      if(fX1 < fX2) return (x>=fX1 && x<=fX2) ? kTRUE : kFALSE;
      else return (x>=fX2 && x<=fX1) ? kTRUE : kFALSE;
    }
    return kFALSE;
    
  }

  Double_t GetA() { return a; }
  Double_t GetB() { return b; }
  Double_t GetC() { return c; }


 private:
  Double_t a;
  Double_t b;
  Double_t c;  

  ClassDef(Segment,1)

};
#endif