// ----------------------------------------------------------------------------
// $Id: LCDFMCCalRecon.cxx,v 1.1 2001/05/10 17:55:55 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDFMCCalRecon.cxx,v $
// Revision 1.1  2001/05/10 17:55:55  toshi
// Changes C++ name convention from *.C to *.cxx to allow Windows usage.
//
//
#include <iostream.h>

#include "LCDFMCCalRecon.h"
#include "TVector3.h"

#include "LCDEDeposit.h"

#define ELEC_TYPE      11  
#define GAMMA_TYPE     22   
#define PIPLUS_TYPE   211 
#define K0SHORT_TYPE  310
#define K0LONG_TYPE   130
#define KPLUS_TYPE    321 
#define PROTON_TYPE  2212
#define NEUTRON_TYPE 2112
#define LAMBDA_TYPE  3122

#define NSMEARABLE  9

#define SMEARABLE_NOT  0
#define SMEARABLE_EM   1
#define SMEARABLE_HAD  2

#define VERY_SMALL 0.0000001
#define MCSTAT_FINALSTATE  1

//______________________________________________________________________
//                                                                       
// LCDFMCCalRecon
//                                                                
// LCDFMCCalRecon is the over-all manager for Fast MC calorimetry.  It must
// be instantiated for each detector (or change to smearing parameters)
// but need only be cleared between events belonging to the same
// detector

ClassImp(LCDFMCCalRecon)   

//------------------------------------------------------------------------
 LCDFMCCalRecon::LCDFMCCalRecon()  { 
  // Constructor.   Primary job is to assemble various parameters
  // describing detector geometry and resolution.
  Init();
}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
 LCDFMCCalRecon::LCDFMCCalRecon(Int_t iseed)  { 
  // Constructor.   Primary job is to assemble various parameters
  // describing detector geometry and resolution.
  Init();
  m_rand.SetSeed(iseed);
}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
 LCDFMCCalRecon::LCDFMCCalRecon(LCDGetParameters* getparam)  { 
  // Constructor.   Primary job is to assemble various parameters
  // describing detector geometry and resolution.
  Init();
  SetDetectorParameters(getparam);
}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
 LCDFMCCalRecon::LCDFMCCalRecon(LCDGetParameters* getparam, Int_t iseed)  { 
  // Constructor.   Primary job is to assemble various parameters
  // describing detector geometry and resolution.
  Init();
  m_rand.SetSeed(iseed);
  SetDetectorParameters(getparam);
}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
 void LCDFMCCalRecon::Init() {
  ClsList= 0;
  // Merge sizes (default)
  merge_siz_EM   =0.8;
  merge_siz_HAD  =1.3;
  merge_angle_EM =0.0200;
  merge_angle_HAD=0.0600;

  Int_t iseed=(Int_t)&m_rand;
  m_rand.SetSeed(iseed);

  m_cluster_merge = kTRUE;  // Default.. do cluster merging
}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
 void LCDFMCCalRecon::SetDetectorParameters(LCDGetParameters* Getparam)  {
  gp     = Getparam;  

  nVolumes     = gp->GetNVolumes();
  if (nVolumes == 0) {
    Int_t ok = gp->SetDetectorVolumes();
    if (ok) {
      nVolumes     = gp->GetNVolumes();
      pInnerVolume = gp->GetInnerVolume();
      pVolumes     = gp->GetDetectorVolumes();
    } else  { // make it clear we have no detector
      nVolumes     = 0;
      pInnerVolume = 0;
      pVolumes     = 0;
      printf(" Error in CalRecon :: Cannot define detector volumes");
    }
  }else{
    pInnerVolume = gp->GetInnerVolume();
    pVolumes     = gp->GetDetectorVolumes();
  }

  // Set cluster merging sizes
  merge_siz_EM   = gp->GetEMMergeSize();
  merge_siz_HAD  = gp->GetHADMergeSize();
  merge_angle_EM = gp->GetEMMergeAngle();
  merge_angle_HAD= gp->GetHADMergeAngle();

}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
void 
 LCDFMCCalRecon::Cleanup(LCDEvent* event) {
  if (event) {
    event->ClusterLst()->Delete();
  }
}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
void 
 LCDFMCCalRecon::Doit(LCDEvent* event) {
  // Process an event.  For each MC particle, check to see whether it's
  // "smearable".  If so, propagate, then smear energy and transverse
  // position.  

  LCDSwimTraj  traj;
  if (!nVolumes) return;  // no detector

  ClsList = event->ClusterLst();
  ClsList = event->ClusterLst();
  ClsList->Delete();

  //----------> create clusters 
  Int_t nParts = event->MCparticles()->GetEntries();
  Int_t nkal =0;

  LCDMcPart*     pPart;
  Int_t          canSmear;
  Bool_t         hadron;
  TVector3       pos0;
  TVector3       mom0;
  Int_t          ifSwim;
  LCDCluster*    cluster;
  TLorentzVector p4mc;
  TVector3       pPos;
  Double_t       energy;
  Double_t       theta;
  Double_t       r;
  Double_t       r_cal;
  Double_t       z_cal;

  for (Int_t i = 0 ; i< nParts ; i++) {
    pPart = (LCDMcPart*)event->MCparticles()->At(i);
    
    canSmear = Smearable(pPart);
 
    if( canSmear ) {
      hadron = (canSmear == SMEARABLE_HAD);
      
      pos0=*(pPart->GetPositionPtr());
      mom0=pPart->Get4MomentumPtr()->Vect();
      Int_t    charge = (Int_t)pPart->GetCharge();
      traj.SetParticleParams(pos0, mom0, charge);
      
      ifSwim = Swim(&traj, hadron);

      if ( ifSwim == 1 ) {

	cluster = new((*event->ClusterLst())[nkal++]) LCDCluster();

	p4mc=*(pPart->Get4MomentumPtr());
	Smear(p4mc.E(),&traj,hadron,cluster);

	pPos=*(traj.GetNewPosPtr());
	energy = p4mc.E();
	theta  = pPos.Theta();
	r      = pPos.Mag();
	
	r_cal = r*TMath::Sin(theta);
	z_cal = TMath::Abs(r_cal*TMath::Cos(theta)/TMath::Sin(theta));
	
	if( r_cal> gp->GetEMBarrelInnerR() &&
	    r_cal< gp->GetHADBarrelInnerR()&& 
	    z_cal< gp->GetEMBarrelOuterZ()   )   { //EM  Barrel 
	  cluster->SetBarEnd(0); cluster->SetSystem(0);
	} else if(r_cal> gp->GetHADBarrelInnerR()){ //HAD Barrel 	   
	  cluster->SetBarEnd(0); cluster->SetSystem(1);
        } else if(r_cal< gp->GetEMBarrelInnerR() && 
		  z_cal< gp->GetEMBarrelOuterZ() ){ //EM  endcap
	  cluster->SetBarEnd(1); cluster->SetSystem(0);
	} else if(r_cal< gp->GetHADBarrelInnerR()&& 
		  z_cal> gp->GetEMBarrelOuterZ() ){ //HAD endcap
	  cluster->SetBarEnd(1); cluster->SetSystem(1);
	}
	
	cluster->AddParticle(i,energy);

      }
    }
  }

  // Ultimately may have another stage where distances between first-pass
  // smears are examined to see if any are close enough (w.r.t. supplied
  // parameters) to be combined into a single smear
  if (m_cluster_merge) {
    DoMerge();
  }
}


 Int_t LCDFMCCalRecon::Smearable(LCDMcPart* pPart) {
  // Check a) that we have a final-state particle b) particle type.  
  // Must be electron/positron, photon, or an acceptable hadron. 
  // c) energy.  Must be non-negligible.

  static Int_t willSmear[NSMEARABLE] = {
    ELEC_TYPE,
    GAMMA_TYPE,
    PIPLUS_TYPE,
    K0SHORT_TYPE,
    K0LONG_TYPE,
    KPLUS_TYPE,
    PROTON_TYPE,
    NEUTRON_TYPE,
    LAMBDA_TYPE};

  Bool_t ok = kFALSE;

  // Check status
  Int_t mcStat = pPart->GetStatus();
  if (mcStat != MCSTAT_FINALSTATE) return SMEARABLE_NOT;

  Int_t partType = TMath::Abs(pPart->GetParticleID());

  for (int ix = 0; ix < NSMEARABLE; ix++) {
    if (partType == willSmear[ix]) {
      ok = kTRUE;
      break;
    }
  }
  if (!ok) return SMEARABLE_NOT;

  // Energy is very small -- forGet it
  if (pPart->Get4MomentumPtr()->E() < VERY_SMALL) return SMEARABLE_NOT;

  return  ((partType == ELEC_TYPE) || (partType == GAMMA_TYPE)) ?
    SMEARABLE_EM  : SMEARABLE_HAD;
}

 Int_t  LCDFMCCalRecon::Swim(LCDSwimTraj* pTraj, Bool_t hadron) {
  // Swim particle with initial conditions described in pTraj through
  // the detector.  Return true if particle showers in the detector
  // (where the coil counts as part of the detector)

  LCDSwimTraj::PropagateRet status;
  Int_t              charge = pTraj->GetCharge();
  LCDDetectorVolume* pCur = (LCDDetectorVolume *) pInnerVolume;
  Double_t           decrS;     // measured in cm
  Double_t           decrPath;
  TVector3           pPos(0.,0.,0.);
  TVector3           pMom(pTraj->GetIniMom());
  LCDDetectorVolume* pNextVolume;

  Int_t           nVolume = 0;       // debug info
  // traj has been initialized with position, momentum, charge

  // First swim through empty central volume
  if ((pCur->field != 0.0) && (charge != 0)) { // helix
    status = pTraj->HelixToCyl(pCur, &decrS);
  } else {  // straight line
    status = pTraj->LineToCyl(pCur, &decrS);
  }

  if ((status == LCDSwimTraj::INTERSECT_NONE) ||
      (status == LCDSwimTraj::INTERSECT_ERROR)   ) return 0;

  // Check that we didn't exit to the air volume next to the beampipe.
  pPos = pTraj->GetNewPos();

  if ((status == LCDSwimTraj::INTERSECT_ZPLUS) || 
      (status == LCDSwimTraj::INTERSECT_ZMINUS)) {
    if (pCur->pNextZ == 0) return 0;
    pNextVolume = (LCDDetectorVolume *) pCur->pNextZ;
    if (pPos.Perp() < pNextVolume->outerR) return 0;
  }

  pTraj->Update();

  // Throw an interaction length or rad. length depending on particle type
  // measured in interaction/rad len
  Double_t           remPath = RanExpDecay();

  Double_t newR=0;
  Double_t newZ=0;
  Bool_t helix;
  while (kTRUE) {
    nVolume++;

    switch(status) {   // find new volume
    case LCDSwimTraj::INTERSECT_ZMINUS:
    case LCDSwimTraj::INTERSECT_ZPLUS:
      pCur = (LCDDetectorVolume *) pCur->pNextZ;
      if (!pCur) return 0;
      // In the special case in which we're coming from an endcap air
      // volume headed in a z-direction, there will typically be more 
      // than one adjacent volume.  pNextZ pointer gives us innermost.
      // Check this against our current radius.  If necessary, follow
      // pNextR pointers until we find an acceptable volume
      newR=pTraj->GetNewPos().Perp();
      while (newR > pCur->outerR) {
	pCur = (LCDDetectorVolume *) pCur->pNextR;
	if (!pCur) return 0;
      }
      break;
    case LCDSwimTraj::INTERSECT_CYLIN:
      // One extra thing to check here as well.  In case we go inwards, 
      // there may be more than one component next to the original. 
      // pPrevR points
      // to the one with smallest z.  If this isn't right (based on
      // new position), follow its pNextZ pointer.      
      newZ = TMath::Abs(pTraj->GetNewPos().Z());
      pCur = (LCDDetectorVolume *) pCur->pPrevR;
      if (!pCur) return 0;
      while (newZ > pCur->outerZ) {
	pCur = (LCDDetectorVolume *) pCur->pNextZ;
	if (!pCur) return 0;
      }
      break;
    case LCDSwimTraj::INTERSECT_CYLOUT:
      pCur = (LCDDetectorVolume *) pCur->pNextR;
      break;
    default:
      return 0;
    }

    if (!pCur) return 0;  // Went outside detector

    // Propagate to boundary of volume
    helix = (pCur->field != 0.0) && (charge != 0);
    if (helix) {
      status = pTraj->HelixToCyl(pCur, &decrS);
    } else {  // straight line
      status = pTraj->LineToCyl(pCur, &decrS);
    }

    if ((status == LCDSwimTraj::INTERSECT_NONE) ||
	(status == LCDSwimTraj::INTERSECT_ERROR)  ) return 0;

    decrPath = hadron ? pCur->iLenPerCm : pCur->radLenPerCm;
    decrPath *= decrS;

    if (decrPath < remPath) { // keep going
      remPath -= decrPath;
      pTraj->Update();
    } else {  // We'll finish in this volume
      decrS = remPath / (hadron ? pCur->iLenPerCm : pCur->radLenPerCm);
      if (helix) {
	status = pTraj->HelixPath(pCur->field, decrS);
      } else {
	status = pTraj->LinePath(decrS);
      }

      nVolume++;
      pTraj->Update();
      return 1;
    }
  }
}

//------------------------------------------------------------------------
 void LCDFMCCalRecon::DoMerge() {
  // do Cluster merge taking into account detector granuarity...

  if (ClsList == 0) return;

  Int_t nCluster = ClsList->GetEntries();
  if (nCluster < 1) return;

  //----------> Sort order of energy...
  ClsList->Sort();

  //----------> Cluster Mearge
  Int_t       ikal;
  LCDCluster* pcali;
  TVector3    xcali;
  Int_t       emhad_i;
  Int_t       barend_i;
  Int_t       jkal;
  LCDCluster* pcalj;
  Int_t       barend_j;
  Int_t       emhad_j;
  TVector3    xcalj;
  Double_t    merge_siz;
  Double_t    merge_angle;
  Double_t    diff_c;
  Double_t    diff_l;
  for (ikal=0 ; ikal < nCluster-1 ; ikal++) {
    pcali=(LCDCluster*)(ClsList->At(ikal));

    emhad_i  = pcali->GetSystem(); 
    barend_i = pcali->GetBarEnd();
    xcali.SetTheta(pcali->GetEnergyTheta());
    xcali.SetPhi(pcali->GetEnergyTheta());
    xcali.SetXYZ(pcali->GetEnergyR()*TMath::Sin(pcali->GetEnergyTheta())
		 *TMath::Cos(pcali->GetEnergyPhi()),
		 pcali->GetEnergyR()*TMath::Sin(pcali->GetEnergyTheta())
		 *TMath::Sin(pcali->GetEnergyPhi()),
		 pcali->GetEnergyR()*TMath::Cos(pcali->GetEnergyTheta()));

    for (jkal=ikal+1 ; jkal < nCluster ; jkal++) {
      pcalj=(LCDCluster*)ClsList->At(jkal);

      barend_j = pcalj->GetBarEnd();
      emhad_j  = pcalj->GetSystem(); 

      if (barend_i != barend_j || emhad_i != emhad_j) continue;

      xcalj.SetXYZ(pcalj->GetEnergyR()*TMath::Sin(pcalj->GetEnergyTheta())
		   *TMath::Cos(pcalj->GetEnergyPhi()),
		   pcalj->GetEnergyR()*TMath::Sin(pcalj->GetEnergyTheta())
		   *TMath::Sin(pcalj->GetEnergyPhi()),
		   pcalj->GetEnergyR()*TMath::Cos(pcalj->GetEnergyTheta()));

      merge_siz   = merge_siz_EM; //EM
      merge_angle = merge_angle_EM; //EM
      if (emhad_i != 0) { // HAD
	merge_siz   = merge_siz_HAD;
	merge_angle = merge_angle_HAD;
      }

      diff_c=xcalj.Angle(xcali);
      diff_l=xcalj.Mag()*diff_c;
      if ( diff_l < merge_siz || diff_c < merge_angle ){//Marge cluster
	AddCluster(pcali,pcalj);
	ClsList->Remove(pcalj);
	//delete pcalj;
      }

    }

    ClsList->Compress();
    nCluster=ClsList->GetEntries();

  }

}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
 void LCDFMCCalRecon::AddCluster(LCDCluster* kal1, LCDCluster* kal2) {
  
  //Energy, Theta, Phi, R
  Double_t energy = kal1->GetEnergy();
  Double_t eTheta = kal1->GetEnergyTheta();
  Double_t ePhi   = kal1->GetEnergyPhi();
  Double_t eR     = kal1->GetEnergyR();
  Double_t ekal   = kal2->GetEnergy();
  Double_t tkal   = kal2->GetEnergyTheta();
  Double_t fkal   = kal2->GetEnergyPhi();
  Double_t rkal   = kal2->GetEnergyR();

  Double_t esum     =0;
  Double_t ephisum  =0;
  Double_t ethetasum=0;
  Double_t ersum    =0;

  esum   = energy + ekal;
  Double_t phiDiff= fkal - ePhi;
  if      (phiDiff >  TMath::Pi()) { phiDiff -= 2.*TMath::Pi();}
  else if (phiDiff < -TMath::Pi()) { phiDiff += 2.*TMath::Pi();}
  ephisum = ekal*phiDiff;
  Double_t theDiff= tkal - eTheta;
  if      (theDiff >  TMath::Pi()) { theDiff -= 2.*TMath::Pi();}
  else if (theDiff < -TMath::Pi()) { theDiff += 2.*TMath::Pi();}
  ethetasum = ekal*theDiff;
  ersum     = energy*eR + ekal*rkal;

  energy   = esum;
  ePhi    += ephisum/esum;
  eTheta  += ethetasum/esum;
  eR       = ersum/esum;

  Double_t phiRMS   = kal1->GetEnergyPhiRMS();
  Double_t thetaRMS = kal1->GetEnergyThetaRMS();
  Double_t rRMS     = kal1->GetEnergyRRMS();;
  //Temporary...T.Abe 4/20/2001
  kal1->SetUpCluster(energy,
		     ePhi,    phiRMS,
		     eTheta,  thetaRMS,
		     eR,      rRMS,
		     kal1->GetBarEnd(),kal1->GetSystem());

  //McList
  Int_t        nMC2=kal2->GetMcListEntries();
  Int_t*    mcList2=kal2->GetMcList();
  Float_t* emcList2=kal2->GetMcListE();
  for (Int_t imc=0 ; imc < nMC2 ; imc++) {
    kal1->AddParticle(mcList2[imc],emcList2[imc]);
  }

}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
 void LCDFMCCalRecon::Smear(Double_t epart, LCDSwimTraj* ptraj,
			Bool_t hadron, LCDCluster* kal) {
  // Inputs are
  //    -- "perfect" particle position, after propagation
  //    -- hadron/lepton boolean
  //    -- smearing parameters
  //    -- random number generator instance

  Double_t a, b;

  if (hadron) {
    a = gp->GetCalEnergyHADA(); b = gp->GetCalEnergyHADB();
  } else {
    a = gp->GetCalEnergyEMA() ; b = gp->GetCalEnergyEMB() ;
  }

  Double_t energySigma = (TMath::Sqrt((a*a)/epart + b*b))*epart;
  Double_t energy = m_rand.Gaus(epart, energySigma);
  if (energy < 0.0) energy = 0.0;

  // Now for position
  // First step should be to calculate shower max position from what
  // I've got (shower start position).  This will add something in
  // the direction of the momentum vector (or do I need to propagate
  // along a helix for charged particles in a field??) but I don't
  // know how to do this yet.

  // Then smear position in plane perp. to momentum vector.
  if (hadron) {
    a = gp->GetCalTransHADA(); b = gp->GetCalTransHADB();
  } else {
    a = gp->GetCalTransEMA();  b = gp->GetCalTransEMB();
  }

  // Find sigma. Use MC truth -- not smeared -- energy
  //  sigma = sqrt((a*a)/E + (b*b))
  // and use it to throw a distance in transverse plane.
  // Also get random phi (in same plane).  Since distance can be
  // negative, just let phi range from 0 to pi.
  Double_t transSigma = TMath::Sqrt((a*a)/epart + b*b);
  Double_t transDist  = m_rand.Gaus(0.0, transSigma);
  Double_t phi        = TMath::Pi()*m_rand.Rndm(0);

  TVector3 xop(transDist*TMath::Cos(phi),transDist*TMath::Sin(phi),0);
  TVector3 nx(1,0,0),ny(0,1,0),nz(0,0,1);
  TVector3 nzp(*(ptraj->GetNewMomPtr())); nzp.SetMag(1.0);
  //TVector3 nzp(*(ptraj->GetNewPosPtr())); nzp.SetMag(1.0);
  TVector3 nxp(nz.Cross(nzp))           ; nxp.SetMag(1.0);
  TVector3 nyp(nzp.Cross(nxp))          ; nyp.SetMag(1.0);
  TVector3 xpos   ((nx*nxp)*xop.X(),(ny*nxp)*xop.X(),(nz*nxp)*xop.X());
  xpos += TVector3((nx*nyp)*xop.Y(),(ny*nyp)*xop.Y(),(nz*nyp)*xop.Y());
  xpos += TVector3((nx*nzp)*xop.Z(),(ny*nzp)*xop.Z(),(nz*nzp)*xop.Z());
  xpos += *(ptraj->GetNewPosPtr());

  //
  Double_t eTheta = xpos.Theta();
  Double_t ePhi   = xpos.Phi();
  Double_t eR     = xpos.Mag();
  //Double_t eR     = ptraj->GetNewPosPtr()->Mag();

  // errors in phi and theta added in quad. should give total trans. err
  Double_t eThetaRMS = transSigma/TMath::Sqrt(2.0);
  Double_t ePhiRMS   = eThetaRMS;
  Double_t eRRMS     = 0.0;

  //
  kal->SetUpCluster(energy,
		    ePhi,   ePhiRMS,
		    eTheta, eThetaRMS,
		    eR,     eRRMS,
		    0, 0);
}