/** CbmMuchDigitize.cxx
 *@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.
 **/


#include <cassert>
#include <vector>

// Includes from ROOT
#include "TClonesArray.h"
#include "TObjArray.h"
//#include "TRandom3.h"
#include "TDatabasePDG.h"

// Includes from base
#include "CbmRootManager.h"
#include "CbmRunAna.h"
#include "CbmRuntimeDb.h"
#include "CbmMCTrack.h"
#include "CbmMCPoint.h"

// Includes from MUCH
#include "CbmGeoMuchPar.h"
#include "CbmMuchDigi.h"
#include "CbmMuchDigiMatch.h"
#include "CbmMuchDigiPar.h"
#include "CbmMuchDigiScheme.h"
#include "CbmMuchDigitize.h"
#include "CbmMuchPoint.h"
#include "CbmMuchSector.h"
#include "CbmMuchStation.h"




// -----   Default constructor   ------------------------------------------
CbmMuchDigitize::CbmMuchDigitize() : CbmTask("MuchDigitize", 1) {
  fGeoPar      = NULL;
  fDigiPar     = NULL;
  fPoints      = NULL;
  fDigis       = NULL;
  fDigiMatches = NULL;
  fUseAvalanche = 0;
  fIntersPoints = new TClonesArray("Point", 8);
  fDigiScheme  = new CbmMuchDigiScheme();
  Reset();
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
CbmMuchDigitize::CbmMuchDigitize(Int_t iVerbose) 
  : CbmTask("MuchDigitize", iVerbose) { 
  fGeoPar      = NULL;
  fDigiPar     = NULL;
  fPoints      = NULL;
  fDigis       = NULL;
  fDigiMatches = NULL;
  fUseAvalanche = 0;
  fIntersPoints = new TClonesArray("Point", 8);
  fDigiScheme  = new CbmMuchDigiScheme();
  Reset();
}
// -------------------------------------------------------------------------



// -----   Constructor with name   -----------------------------------------
CbmMuchDigitize::CbmMuchDigitize(const char* name, Int_t iVerbose) 
  : CbmTask(name, iVerbose) { 
  fGeoPar      = NULL;
  fDigiPar     = NULL;
  fPoints      = NULL;
  fDigis       = NULL;
  fDigiMatches = NULL;
  fUseAvalanche = 0;
  fIntersPoints = new TClonesArray("Point", 8);
  fDigiScheme  = new CbmMuchDigiScheme();
  Reset();
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
CbmMuchDigitize::~CbmMuchDigitize() { 
  if ( fGeoPar)    delete fGeoPar;
  if ( fDigiPar)   delete fDigiPar;
  if ( fDigis ) {
    fDigis->Delete();
    delete fDigis;
  }
  if ( fDigiMatches ) {
    fDigiMatches->Delete();
    delete fDigiMatches;
  }
  if ( fIntersPoints ) {
    fIntersPoints->Delete();
    delete fIntersPoints;
  }
  if ( fDigiScheme ) delete fDigiScheme;
  fUseAvalanche = 0;
  Reset();
}
// -------------------------------------------------------------------------

// ------- Private method ExecSimple ---------------------------------------
Bool_t CbmMuchDigitize::ExecSimple(CbmMuchStation* station, CbmMuchPoint* point,
				   Int_t iPoint){

  Double_t xIn = point->GetXIn();
  Double_t yIn = point->GetYIn();
  Double_t xOut = point->GetXOut();
  Double_t yOut = point->GetYOut();
  
  // Get avalanche center position
  Double_t x0  = (xIn+xOut)/2;
  Double_t y0  = (yIn+yOut)/2;

  Int_t nSectors = station->GetNSectors();
  Int_t iChannel = -1;
  CbmMuchSector* sector = NULL;
  for (Int_t iSector=0; iSector<nSectors; iSector++) {
    sector = station->GetSector(iSector);
    iChannel = sector->GetChannel(x0, y0);
    if ( iChannel < 0 ) continue;
    break;
  }
  if ( iChannel < 0 ){
    fNOutside++;
    return kFALSE;
  }

  Int_t iDetectorId  = sector->GetDetectorId();
  Int_t stationNr = station->GetStationNr();
  Int_t sectorNr = sector->GetSectorNr();

  pair<Int_t, Int_t> aa(iDetectorId, iChannel);
  Int_t iDigi = -1;
  if ( fChannelMap.find(aa) == fChannelMap.end() ) {
    // Channel not yet active. Create new Digi and Match.
    iDigi = fNDigis;
    new ((*fDigis)[fNDigis]) CbmMuchDigi(stationNr, sectorNr, iChannel, point->GetTime());
    new ((*fDigiMatches)[fNDigis]) CbmMuchDigiMatch(iPoint);
    // Match pair (iDetectorId, iChannel) to index of the Digi.
    fChannelMap[aa] = iDigi;
    fNDigis++;
  }
  else {
    // Channel already active. 
    iDigi = fChannelMap[aa];
    CbmMuchDigi* digi = dynamic_cast<CbmMuchDigi*>(fDigis->At(iDigi));
    digi->AddTime(point->GetTime()); // add time info
    digi->AddDTime(); // add time resolution info
    CbmMuchDigiMatch* match = dynamic_cast<CbmMuchDigiMatch*>(fDigiMatches->At(iDigi));
    if ( match ) {
      match->AddPoint(iPoint);
      fNMulti++;
    }
  }
  return kTRUE;
}
// -------------------------------------------------------------------------

// ------- Private method ExecAdvanced -------------------------------------
Bool_t CbmMuchDigitize::ExecAdvanced(CbmMuchStation* station, CbmMuchPoint* point, 
				     Int_t iPoint){

  Int_t stationNr = station->GetStationNr();
  
  //**************  Simulate avalanche (begin) *********************************//
  
  // Get track length within the station
  Double_t xIn    = point->GetXIn();
  Double_t yIn    = point->GetYIn();
  Double_t xOut   = point->GetXOut();
  Double_t yOut   = point->GetYOut();
  Double_t deltaX = xOut-xIn;
  Double_t deltaY = yOut-yIn;
  Double_t lTrack = sqrt(deltaX*deltaX+
				deltaY*deltaY); // track length

  //********** Primary electrons from the track (begin) ************************//
  // Get particle's characteristics
  Int_t trackID = point->GetTrackID();
  if(trackID < 0) return kFALSE;
  CbmMCTrack* mcTrack = (CbmMCTrack*)fMCTracks->At(trackID);
  if(NULL == mcTrack) return kFALSE;
  Int_t pdgCode = mcTrack->GetPdgCode(); 

  // Reject funny particles
  if(pdgCode == 10010020 ||
     pdgCode == 10010030 ||
     pdgCode == 50000050 ||
     pdgCode == 50010051 ||
     pdgCode == 10020040) return kFALSE;
  
  TParticlePDG *particle = TDatabasePDG::Instance()->GetParticle(pdgCode);
  TVector3 momentum;                            // 3-momentum of the particle
  point->Momentum(momentum); 
  Double_t mom = momentum.Mag()*1e3;            // absolute momentum value [MeV/c]
  Double_t mom2 = mom*mom;                      // squared momentum of the particle
  Double_t mass = particle->Mass()*1e3;         // mass of the particle [MeV/c^2]
  Double_t mass2 = mass*mass;                   // squared mass of the particle
  Double_t Tkin = sqrt(mom2 + mass2) - mass;    // kinetic energy of the particle
  Tkin = Tkin*0.511/mass;                       // scaled kinetic energy with respect to electron
  Double_t logT = log(Tkin);
  Double_t sigma = e_sigma_n_e(logT);           // sigma for Landau distribution
  Double_t mpv = e_MPV_n_e(logT);               // most probable value for Landau distr.
  TRandom lanRnd;
  Int_t nElectrons = Int_t(lanRnd.Landau(mpv, sigma)); // number of prim. electrons per 0.5 cm gap
  nElectrons = Int_t(Double_t(nElectrons)*lTrack/0.5);           // number of prim.electrons for current track length  

//   cout << "DEBUG: nElectrons = " << nElectrons << "\t"
//        << "lTrack = " << lTrack << "cm\t"
//        << "Tkin = " << Tkin << "MeV\t"
//        << endl;
  //********** Primary electrons from the track (end)   ***********************//

  Double_t cosphi_tr = deltaX/lTrack;     // cos of track polar angle
  Double_t sinphi_tr = deltaY/lTrack;     // sin of track polar angle

  Int_t iChannel = -1;
  //  TRandom3 rnd;
  for(Int_t iElectron = 0; iElectron < nElectrons; iElectron++){ // loop over all primary electrons
    // Calculate coordinates of pr. electrons
    Double_t aL = fRnd.Rndm()*lTrack;
    Double_t xe = xIn + aL*cosphi_tr;
    Double_t ye = yIn + aL*sinphi_tr;
    // Calculate number of sec. electrons for each pr. electron
    Int_t nSecElectrons = GasGain();
    Double_t spotRad = 1e-5*nSecElectrons + 0.006; // sec. electrons' spot radius [cm]
    Double_t spotL = spotRad*2;                    // square for the spot

    //********** Calculate charge collected on pads (begin)  ********************//
    Int_t nSectors = station->GetNSectors();
    CbmMuchSector *sector = NULL;

    for(Int_t iSector=0; iSector<nSectors; iSector++){  // loop over all sectors
      sector = station->GetSector(iSector);
      Int_t iDetectorId  = sector->GetDetectorId();
      Int_t sectorNr     = sector->GetSectorNr();
      //      Double_t rot       = sector->GetRotation();
      Double_t xc        = sector->GetX0();
      Double_t yc        = sector->GetY0();
      Double_t lx        = sector->GetLx();
      Double_t ly        = sector->GetLy();
      Double_t dx        = sector->GetDx();
      Double_t dy        = sector->GetDy();
      Double_t sinrot    = sector->GetSinRot();
      Double_t cosrot    = sector->GetCosRot();

      // Search of spot-sector intersection
      Double_t padRad    = sqrt(dx*dx + dy*dy)/2; // pad ardius
      Double_t sectorRad = sqrt(lx*lx + ly*ly)/2; // sector radius
      Double_t distance  = sqrt((xe-xc)*(xe-xc) +   
				       (ye-yc)*(ye-yc)); // distance between sector center and pr. electron
      // Rough (fast) search
      if(distance > sectorRad + spotRad) continue;

      // Run detailed search, if rough one passed
      //      continue;
      Double_t area = 0.;
      TPolyLine secPolygon = GetPolygon(xc,yc,lx,ly,sinrot, cosrot);
      TPolyLine spotPolygon = GetPolygon(xe, ye, spotL, spotL, sinrot, cosrot);
      if(!PolygonsIntersect(secPolygon, spotPolygon, area)) continue;
      continue;

//       cout << "DEBUG: Detailed sector search passed:\n"
// 	   << "\tiSector = " << iSector + 1 << "\txc = " << xc << "\tyc = " << yc << "\t"
// 	   << "\tlx = " << lx << "\tly = " << ly << "\n"
// 	   << "DEBUG: Fired channels:\n"
// 	   << endl;
//       // Spot boundaries (global c.s.)
//       Double_t xMin = xe - spotL/2.;
//       Double_t yMin = ye - spotL/2.;
//       Double_t xMax = xe + spotL/2.;
//       Double_t yMax = ye + spotL/2.;
//       // Spot boundaries (sector internal c.s.)
//       Double_t xIntMin, yIntMin, xIntMax, yIntMax;
//       sector->IntCoord(xMin, yMin, xIntMin, yIntMin);
//       sector->IntCoord(xMax, yMax, xIntMax, yIntMax);
//       Int_t nColMin = (xIntMin < 0.) ? 0 : (Int_t)((xIntMin+1e-5)/dx);
//       Int_t nColMax = (xIntMin > lx) ? (Int_t)((lx+1e-5)/dx)
//                                    : (Int_t)((xIntMax+1e-5)/dx);
//       Int_t nRowMin = (yIntMin < 0.) ? 0 : (Int_t)((yIntMin+1e-5)/dy);
//       Int_t nRowMax = (yIntMin > ly) ? (Int_t)((ly+1e-5)/dy)
//                                    : (Int_t)((yIntMax+1e-5)/dy);
      // Fire the sector intersected with the spot
      Int_t nColumns = (Int_t)((lx + 1e-5)/dx);
      Int_t nRows = (Int_t)((ly + 1e-5)/dy);
      for(int iCol=0; iCol<nColumns; iCol++){
	for(int iRow=0; iRow<nRows; iRow++){
	  Double_t xint = (Double_t(iCol)+0.5 )*dx;
	  Double_t yint = (Double_t(iRow)+0.5 )*dy;
	  // Translation to centre of the sector
	  xint = xint - lx/2.;
	  yint = yint - ly/2.;  
	  // Rotation around the sector centre
	  Double_t x = xint * cosrot - yint * sinrot;
	  Double_t y = xint * sinrot + yint * cosrot;
	  // Translation into global c.s.
	  x = x + xc;
	  y = y + yc;

	  Double_t dist = TMath::Sqrt(TMath::Power(xe-x,2) +
				      TMath::Power(ye-y,2));   // distance between pad center and pr. electron
	  // Rough (fast) search
	  if(dist > spotRad + padRad) continue;
	  
	  // Fire pads intersected with the avalanche
	  TPolyLine padPolygon = GetPolygon(x, y, dx, dy, sinrot, cosrot);
	  if(!PolygonsIntersect(padPolygon, spotPolygon, area)) continue;
	  iChannel = iCol + iRow*nColumns;

	  // Create digis for fired pads
	  pair<Int_t, Int_t> aa(iDetectorId, iChannel);
	  Int_t iDigi = -1;
	  if ( fChannelMap.find(aa) == fChannelMap.end() ) {
	    // Channel not yet active. Create new Digi and Match.
	    iDigi = fNDigis;
	    new ((*fDigis)[fNDigis]) CbmMuchDigi(stationNr, sectorNr, iChannel, point->GetTime());
	    new ((*fDigiMatches)[fNDigis]) CbmMuchDigiMatch(iPoint);
	    // Match pair (iDetectorId, iChannel) to index of the Digi.
	    fChannelMap[aa] = iDigi;
	    fNDigis++;
	  }
	  else {
	    // Channel already active.
	    iDigi = fChannelMap[aa];
	    CbmMuchDigi* digi = dynamic_cast<CbmMuchDigi*>(fDigis->At(iDigi));
	    digi->AddTime(point->GetTime()); // add time info
	    digi->AddDTime(); // add time resolution info
	    CbmMuchDigiMatch* match = dynamic_cast<CbmMuchDigiMatch*>(fDigiMatches->At(iDigi));
	    if ( match ) {
	      match->AddPoint(iPoint);
	      fNMulti++;
	    }
	  }
	}
      }
    } // loop over all sectors
    //********** Calculate charge collected on pads (end)    ********************//
  } // loop over all primary electrons
  
  if ( iChannel < 0 ){
    //    fNOutside++;
    return kFALSE;
  }
  return kTRUE;
  
  //**************  Simulate avalanche (end) **********************************//

}
// -----   Public method Exec   --------------------------------------------
void CbmMuchDigitize::Exec(Option_t* opt) {
  // Reset all eventwise counters
  fTimer.Start();
  Reset();

  // Verbose screen output
  if ( fVerbose > 2 ) {
    cout << endl << "-I- " << fName << ": executing event" << endl;
    cout << "----------------------------------------------" << endl;
  }

  // Check for input arrays
  if( ! fMCTracks ) {
    cout << "-W- " << fName << "::Exec: No input array (MCTrack)" 
	 << endl;
    cout << "- " << fName << endl;
    return;
  }
  if ( ! fPoints ) {
    cerr << "-W- " << fName << "::Exec: No input array (MuchPoint) "
	 << endl;
    cout << "- " << fName << endl;
    return;
  }

  Int_t notUsable = 0; // DEBUG: counter for not usable points

  // Loop over all MuchPoints
  for (Int_t iPoint=0; iPoint<fPoints->GetEntriesFast(); iPoint++) {
    CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(iPoint);
    fNPoints++;

    // Get the station the point is in
    Int_t mcId = point->GetDetectorID();
    CbmMuchStation* station = fDigiScheme->GetStationByMcId(mcId);
    if ( ! station ) {
      if ( fVerbose > 2 ) cout << "MuchPoint " << iPoint 
			       << ", no digitization (MC volume "
			       << mcId << ")" << endl;
      fNFailed++;
      continue;
    }
    // Take only usable points
    if ( ! point->IsUsable() ) {
      notUsable++;
      continue;
    }

    // Produce Digis
    if(!fUseAvalanche){
      if(!ExecSimple(station, point, iPoint)) continue;
    }
    else {
      if(!ExecAdvanced(station, point, iPoint)) continue; 
    }
  }  // MuchPoint loop
  // Screen output
  fTimer.Stop();
  if ( ! fVerbose ) cout << "+ ";
  else              cout << "-I- ";
  cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8)
       << fixed << right << fTimer.RealTime()
       << " s, points " << fNPoints << ", failed " << fNFailed
       << ", not usable " << notUsable
       << ", outside " << fNOutside << ", multihits " << fNMulti
       << ", digis " << fNDigis << endl;
}
// -------------------------------------------------------------------------

 
// -----   Private method SetParContainers   -------------------------------
void CbmMuchDigitize::SetParContainers() {

  // Get run and runtime database
  CbmRunAna* run = CbmRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  CbmRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

  // Get MUCH geometry parameter container
  fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");

  // Get MUCH digitisation parameter container
  fDigiPar = (CbmMuchDigiPar*) db->getContainer("CbmMuchDigiPar");
 
}
// -------------------------------------------------------------------------



// -----   Private method Init   -------------------------------------------
InitStatus CbmMuchDigitize::Init() {

  CbmRootManager* ioman = CbmRootManager::Instance();
  if ( ! ioman ) Fatal("Init", "No CbmRootManager");

  // Get input array of MuchPoints
  fPoints = (TClonesArray*) ioman->GetObject("MuchPoint");

  // Get input array of MC tracks
  fMCTracks = (TClonesArray*)ioman->GetObject("MCTrack");

  // Register output array MuchDigi
  fDigis = new TClonesArray("CbmMuchDigi",1000);
  ioman->Register("MuchDigi", "Digital response in MUCH", fDigis, kTRUE);

  // Register output array MuchDigiMatches
  fDigiMatches = new TClonesArray("CbmMuchDigiMatch",1000);
  ioman->Register("MuchDigiMatch", "Digi Match in MUCH", fDigiMatches, kTRUE);

  // Build digitization scheme
  if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
    if      (fVerbose >= 1) fDigiScheme->Print();
    cout << "-I- " << fName << "::Init: "
	 << "MUCH digitization scheme succesfully initialised" << endl;
    cout << "    Stations: " << fDigiScheme->GetNStations() 
	 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: " 
	 << fDigiScheme->GetNChannels() << endl;

    if(!fUseAvalanche)
      cout << "Using simple algorithm for producing MUCH hits..." << endl;
    else
      cout << "Using advanced algorithm for producing MUCH hits..." << endl;
    return kSUCCESS;
  }
  return kERROR;
}
// -------------------------------------------------------------------------



// -----   Private method ReInit   -----------------------------------------
InitStatus CbmMuchDigitize::ReInit() {
  // Clear digitisation scheme
  fDigiScheme->Clear();

  // Build new digitisation scheme
  if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
    if      (fVerbose >= 1) fDigiScheme->Print();
    cout << "-I- " << fName << "::Init: "
	 << "MUCH digitization scheme succesfully initialised" << endl;
    cout << "    Stations: " << fDigiScheme->GetNStations() 
	 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: " 
	 << fDigiScheme->GetNChannels() << endl;
    return kSUCCESS;
  }

  return kERROR;
}
// -------------------------------------------------------------------------



// -----   Private method Reset   ------------------------------------------
void CbmMuchDigitize::Reset() {
  fNPoints = fNFailed = fNOutside = fNMulti = fNDigis = 0;
  fChannelMap.clear();
  if ( fDigis ) fDigis->Clear();
  if ( fDigiMatches ) fDigiMatches->Clear();
}
// -------------------------------------------------------------------------



TPolyLine CbmMuchDigitize::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  ){
  Double_t a,b,d,cosphi,sinphi;
  xMin = 1000; xMax = -1000;
  yMin = 1000; yMax = -1000;
  
  
  d = sqrt(pow(xOut-xIn,2)+pow(yOut-yIn,2)); // distance between points
  Double_t x0  = (xIn+xOut)/2;       // x coordinate of the spread center
  Double_t y0  = (yIn+yOut)/2;       // y coordinate of the spread center
  a  = r + d/2;             // half-height of the spread
  b  = r;                   // half-width of the spread

  Double_t x[5]   , y[5];    // coordinates in the spread center r.f.
  Double_t xRot[5], yRot[5]; // coordinates in the global r.f.

  x[1]=-a; y[1]=+b; x[2]=+a; y[2]=+b;
  x[0]=-a; y[0]=-b; x[3]=+a; y[3]=-b;

  cosphi = (xOut-xIn)/d;     // cos of the spread angle in xy plane
  sinphi = (yOut-yIn)/d;     // sin of the spread angle in xy plane

  // rotate and translate the spread to the global r.f.
  for (Int_t i=0;i<4;i++){
    xRot[i] = x0 + x[i]*cosphi - y[i]*sinphi;
    yRot[i] = y0 + x[i]*sinphi + y[i]*cosphi;
    // Find boundaries of the avalanche
    if(xMin > xRot[i]) xMin = xRot[i];
    if(xMax < xRot[i]) xMax = xRot[i];
    if(yMin > yRot[i]) yMin = yRot[i];
    if(yMax < yRot[i]) yMax = yRot[i];
  }
  xRot[4]=xRot[0]; yRot[4]=yRot[0];

  TPolyLine pline(5, xRot, yRot);
  return pline;
}

TPolyLine CbmMuchDigitize::GetPolygon(Double_t x0, Double_t y0, 
				      Double_t width, Double_t height, 
				      Double_t sinrot, Double_t cosrot){

  Double_t x[5], y[5];
  x[1]=-width/2; y[1]=+height/2;      x[2]=+width/2; y[2]=+height/2;
  x[0]=-width/2; y[0]=-height/2;      x[3]=+width/2; y[3]=-height/2;
  x[4]=x[0];  y[4]=y[0];
  for (Int_t i=0;i<5;i++) {
    Double_t xRot = x0+x[i]*cosrot-y[i]*sinrot;
    Double_t yRot = y0+y[i]*cosrot+x[i]*sinrot;
    x[i] = xRot;
    y[i] = yRot;
  }
  TPolyLine pline(5,x,y);
  return pline; 
}

Double_t CbmMuchDigitize::Angle2D(double x1, double y1, double x2, double y2){
    double dtheta,theta1,theta2;
   
  theta1 =TMath::ATan2(y1,x1);
  theta2 = TMath::ATan2(y2,x2);
  dtheta = theta2 - theta1;
  while (dtheta > TMath::Pi())
    dtheta -= TMath::TwoPi();
  while (dtheta < -TMath::Pi())
    dtheta += TMath::TwoPi();
								      
  return(dtheta);
}

Bool_t CbmMuchDigitize::InsidePolygon(TPolyLine polygon, Point p){
  double angle=0.;
  Point p1;
  Point p2;
  Double_t *x = polygon.GetX();
  Double_t *y = polygon.GetY();
  Int_t n = polygon.GetN();
       
  for (int i=0;i<n;i++) {
    p1.SetX(x[i] - p.GetX());
    p1.SetY(y[i] - p.GetY());
    p2.SetX(x[(i+1)%n] - p.GetX());
    p2.SetY(y[(i+1)%n] - p.GetY());
    angle += Angle2D(p1.GetX(),p1.GetY(),p2.GetX(),p2.GetY());
  }
			

  if (TMath::Abs(angle) < TMath::Pi()) return kFALSE;
  else return kTRUE;
}

Double_t CbmMuchDigitize::Area(map<Double_t, Point*> pointsMap){
  if((int)pointsMap.size()==0) return 0.;
  map<Double_t, Point*>::iterator it ;  
  map<Double_t, Point*>::iterator itnext;  
  
  Double_t area = 0.0;

  for(itnext=it=pointsMap.begin(), itnext++; it!=pointsMap.end(); it++, itnext++ ){
    if(itnext==pointsMap.end()) itnext = pointsMap.begin();
    Double_t x = it->second->GetX();
    Double_t y = it->second->GetY();
    Double_t xnext = itnext->second->GetX();
    Double_t ynext = itnext->second->GetY();
    area += x*ynext - xnext*y;
  }
//   if(it) delete it;
//   if(itnext) delete itnext;
  return TMath::Abs(area)/2.;
}

Bool_t CbmMuchDigitize::SegmentsIntersect(Segment s1, Segment s2, 
					  Double_t& iX, Double_t& iY){
  Double_t det = s1.GetA()*s2.GetB() - s1.GetB()*s2.GetA();
  if(TMath::Abs(det) <= 1e-5) return kFALSE;

  Double_t x = -(s1.GetC()*s2.GetB() - s2.GetC()*s1.GetB())/det;
  Double_t y = -(s1.GetA()*s2.GetC() - s2.GetA()*s1.GetC())/det;
  
  if(s1.ContainsPoint(x, y) && s2.ContainsPoint(x, y)){
    iX = x; iY = y;
    return kTRUE;
  }
  return kFALSE;
}

map<Double_t, Point*> CbmMuchDigitize::SortByPhi(TClonesArray *points, Point cp){ 
  map<Double_t, Point*> pointsMap;

  for(int i=0;i<points->GetEntriesFast();i++){
    Point *p = dynamic_cast<Point*>(points->At(i));
    Double_t y = p->GetY() - cp.GetY();
    Double_t x = p->GetX() - cp.GetX();
    Double_t phi = TMath::ATan2(y,x);
    if(y<0) phi += TMath::TwoPi();
    pointsMap[phi] = p;
  }
  return pointsMap;
}

Bool_t CbmMuchDigitize::PolygonsIntersect(TPolyLine polygon1, TPolyLine polygon2, 
					  Double_t& area){

  // Find intersection points
  fIntersPoints->Clear();
  Double_t *x1 = polygon1.GetX();
  Double_t *y1 = polygon1.GetY();
  Double_t *x2 = polygon2.GetX();
  Double_t *y2 = polygon2.GetY();
  Point p;

  Int_t nIntersPoints = 0;
  for(int i = 0; i<polygon1.GetN()-1; i++){
    Segment s1(x1[i],y1[i],x1[i+1],y1[i+1]);
    for(int j=0; j<polygon2.GetN()-1; j++){      
      Segment s2(x2[j],y2[j],x2[j+1],y2[j+1]);
      Double_t fx, fy;
      if(SegmentsIntersect(s1,s2,fx,fy))
	new ((*fIntersPoints)[nIntersPoints++]) Point(fx, fy);            
    }
  }

  // Find if the vertices of one rectangle are contained in another one
  for(int i = 0; i<polygon1.GetN()-1; i++){
    p = Point(x2[i], y2[i]);
    if(InsidePolygon(polygon1, p))
      new ((*fIntersPoints)[nIntersPoints++]) Point(p);
  }
  for(int i = 0; i<polygon2.GetN()-1; i++){
    p = Point(x1[i], y1[i]);
    if(InsidePolygon(polygon2, p))
      new ((*fIntersPoints)[nIntersPoints++]) Point(p);
  }

  if(fIntersPoints->GetEntriesFast() == 0) return kFALSE;
  

  //   Find center  
  Point cp(0.,0.);
  int n = fIntersPoints->GetEntriesFast();
  if(n>0){
    for(int i =0;i<n;i++){
      Point* pt = (Point*) fIntersPoints->At(i);
      cp.SetX( cp.GetX() + pt->GetX() );
      cp.SetY( cp.GetY() + pt->GetY() );
    }
    cp.SetX( cp.GetX() / n);
    cp.SetY( cp.GetY() / n);
  }
 
  // Sort points
  map<Double_t, Point*> pointsMap = SortByPhi(fIntersPoints, cp);

  // Calculate area
  area = Area(pointsMap);

  // Free memory
  if(pointsMap.size()>0) pointsMap.clear();
  

  return kTRUE;
}

Int_t CbmMuchDigitize::GasGain(){
  const Double_t q_mean  = -1.e4;  // mean gas gain, arbitrary value
  return (int) (q_mean*TMath::Log(1-fRnd.Rndm()));
}


Double_t CbmMuchDigitize::e_sigma_n_e(Double_t &logT){
  if (logT<-3.21888) logT=-3.21888; 
  if (logT>9.21034) logT=9.21034; 
  const Int_t n = 7;
  Double_t    val=0;
  Double_t    arg=1;
  Double_t    p[n]={4.28999
		  , -0.179917
		  , 0.479522
		  , -0.152795
		  , 0.0247173
		  , -0.00206182
		  , 6.82384e-05 };
  for (Int_t i=0; i<n; i++){
    val=val+arg*p[i];
    arg = arg*logT;
  }
  return val;
}

Double_t CbmMuchDigitize::e_MPV_n_e(Double_t &logT){
  if (logT<-3.21888) logT=-3.21888;
  if (logT>9.21034) logT=9.21034;
  const Int_t n = 7;
  Double_t    val=0;
  Double_t    arg=1;
  Double_t    p[n]={15.7884
		  , -0.721914
		  , 2.47486
		  , -0.951817
		  , 0.183168
		  , -0.0177078
		  , 0.000665668 };
  for (Int_t i=0; i<n; i++){
    val=val+arg*p[i];
    arg = arg*logT;
  }
  return val;
}

ClassImp(CbmMuchDigitize)
ClassImp(Segment)
ClassImp(Point)