// $Header: CalRecon.cxx $

#include "CalRecon.h"
#include "CalSmeared.h"

#define  false  0
#define  true   1

#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 PIZERO_TYPE 111     // not included among the smearables

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

//______________________________________________________________________
//                                                                       
// CalRecon
//                                                                
// CalRecon 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(CalRecon)   

 CalRecon::CalRecon(const TString& name)  { 
  // Constructor.   Primary job is to assemble various parameters
  // describing detector geometry and resolution.

  Int_t ok;
  ok = GetDetector(name);
  if (ok) {
    //  InitSmearPar();
    smears = new TObjArray(10);
  }
  else { // make it clear we have no detector
    nVolumes = 0;
    pInnerVolume = 0;
    pVolumes = 0;
  }
}
      
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 void CalRecon::cleanup() {
  // Get rid of allocated memory, obsolete pointers.

  // Get rid of reference to old MCparticles
  MCparticles = 0;

  // Get rid of old smears if any
  smears->Delete(0);
  
};

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

  SmearTraj  traj;
  Bool_t     ok;

  MCparticles = event->MCparticles();
  nPart = MCparticles->GetLast() + 1;

  if (!nVolumes) return;  // no detector
  
  for (Int_t ix = 0; ix < nPart; ix++) { // Process particle with index ix.
    
    McPart*       pPart = (McPart*) (*MCparticles)[ix];
    CalSmeared*   pSmeared;
    Int_t         canSmear;

    if (canSmear = Smearable(pPart)) {
      Bool_t hadron = (canSmear == SMEARABLE_HAD);
      SmearVolume *pCur;
      traj.SetParticle(pPart, ix);
      if (pCur = Swim(&traj, hadron)) {
	pSmeared = new CalSmeared(traj, hadron, pCur->pFuzz, pRand); 
    
        smears->Add(pSmeared);
      }
    }

    // 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
  }
};


//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 void CalRecon::spew(FILE* ofile) const
{
  // Output standard cluster quantities to ascii file.

  if (!nVolumes) return;  // no detector

  Int_t  nSmear = smears->GetLast() + 1; 
  fprintf(ofile, "Clusters n %i n", nSmear);
  //  fprintf(ofile, "%d MC particles produced %d smeared clustersn", 
  //	  nPart, nSmear);
  for (Int_t iSmear = 0; iSmear < nSmear; iSmear++) {
    CalSmeared* pSmear = (CalSmeared *) (*smears)[iSmear];
    pSmear->spew(ofile);
  }
  fprintf(ofile, "endn");
}

#define VERY_SMALL 0.0000001
#define MCSTAT_FINALSTATE  1
 Int_t CalRecon::Smearable(McPart* 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,
    //    PIZERO_TYPE,
    K0SHORT_TYPE,
    K0LONG_TYPE,
    KPLUS_TYPE,
    PROTON_TYPE,
    NEUTRON_TYPE,
    LAMBDA_TYPE};

  Bool_t ok = false;

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

  Int_t partType = pPart->GetType();
  if (partType < 0) {
    partType = - partType;
  }

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

  // Energy is very small -- forget it
  Float_t* iniMom = pPart->GetMomentum();
  if (*(iniMom + 3)  < VERY_SMALL) return SMEARABLE_NOT;

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

 Int_t CalRecon::GetDetector(const TString& name) {
  // Somehow get access to needed detector parameters.  For now, use the
  // name to identify a set of parameters.  Could also use it to identify
  // a file, or to just ask (from some parameter server class) for information
  // about the "ambient" detector, if any.
  //
  // Result of this call is to initialize the following fields in CalRecon:
  //       nBarrel   pBarrels
  //       nEndcap   pEndcaps
  // so that at the end have a list of pointers to barrels, arranged in
  // order of increasing R, and a list of endcap in order of increasing z
  
  pRand = new TRandom();

  if (name.CompareTo("Small") == 0)   GetSmall();
  else if (name.CompareTo("Large") == 0) GetLarge();
  else {
    printf("Unsupported detectorn");
    return 0;
  }

  return 1;
}

 SmearVolume*  CalRecon::Swim(SmearTraj* 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)

  SmearTraj::PropagateRet status;
  Int_t                   charge = pTraj->GetCharge();
  SmearVolume*            pCur = (SmearVolume *) pInnerVolume;
  Double_t                decrS;              // measured in cm
  Double_t      totalPath, remPath; // measured in interaction/rad len
  Double_t      decrPath;
  Double_t*     pPos;
  Double_t*     pMom;
  SmearVolume*  pNextVolume;

  int           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 == SmearTraj::INTERSECT_NONE) ||
      (status == SmearTraj::INTERSECT_ERROR)   ) return 0;

  // Check that we didn't exit to the air volume next to the beampipe.
  pPos = pTraj->GetNewPos();
  if ((status == SmearTraj::INTERSECT_ZPLUS) || 
      (status == SmearTraj::INTERSECT_ZMINUS)) {
    Double_t x = *pPos;
    Double_t y = *(pPos + 1);

    if (pCur->pNextZ == 0) return 0;
    pNextVolume = (SmearVolume *) pCur->pNextZ;
    if (TMath::Sqrt(x*x + y*y) < pNextVolume->outerR) return 0;
  }

  pTraj->Update();

  // We can update mcPart, but first we have to change doubles to floats
  pMom = pTraj->GetNewMom();
  {
    Float_t pos[3];
    Float_t mom[3];

    for (int i = 0; i < 3; i++) {
      pos[i] = *(pPos + i);
      mom[i] = *(pMom + i);
    }
    pTraj->GetMc()->SetParticleAtCal(&mom[0], &pos[0]);
  }

  // Throw an interaction length or rad. length depending on particle type
  totalPath = remPath = RanExpDecay(pRand);

  while (true) {
    Bool_t helix;
    nVolume++;
    switch(status) {   // find new volume
    case SmearTraj::INTERSECT_ZMINUS:
    case SmearTraj::INTERSECT_ZPLUS:
      pCur = (SmearVolume *) 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
      {
	Double_t *newPos = pTraj->GetNewPos();
	Double_t x = *newPos;
	Double_t y = *(newPos + 1);
	Double_t newR = TMath::Sqrt(x*x + y*y);

	while (newR > pCur->outerR) {
	  pCur = (SmearVolume *) pCur->pNextR;
	  if (!pCur) return 0;
	}
      }
      break;
    case SmearTraj::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.
      
      {
	Double_t *newPos = pTraj->GetNewPos();
	Double_t newZ = TMath::Abs(*(newPos + 2));
	pCur = (SmearVolume *) pCur->pPrevR;
	if (!pCur) return 0;

	while (newZ > pCur->outerZ) {
	  pCur = (SmearVolume *) pCur->pNextZ;
	  if (!pCur) return 0;
	}
      }
      break;
    case SmearTraj::INTERSECT_CYLOUT:
      pCur = (SmearVolume *) 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 == SmearTraj::INTERSECT_NONE) ||
	(status == SmearTraj::INTERSECT_ERROR)  ) return 0;

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

    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 pCur;
    }
  }
}