//////////////////////////////////////////////////////////////////////////
//
// Function: getScintillatorEnergy.cc
// Purpose: Implementation of getScintillatorEnergy non-member function.
//          
// Created: 25/05/99  Pierre Savard  (copied code from Anwar Bhatti)
// History:  20/02/01 P. Savard: use 3 different codes for spike killings
//                                             0 = run1 default
//                                             1 = run2 simple
//                                             2 = run1 fancy
//           28/09/01 B. Heinemann: implement versioning to account for changeing
//                                  attentuation length and sigmas
//           29/04/03 A. Wyatt: if a pmt is dead, take the energy from the non-dead pmt
//////////////////////////////////////////////////////////////////////////

#include <cassert>
#include <iostream>
#include "CalorGeometry/CalParameters.hh"
#include "CalorObjects/getScintillatorEnergy.hh"
#include "CalorObjects/CalThresholds.hh"

// namespace calor{

void getScintillatorEnergy(int detector, float ecal0, float ecal1, 
			   float& ecal, float& phical, int& flag, int& killSpike, bool badpmt0, bool badpmt1,
			   int version){
  // flag flags the towers which are spikes
    flag =0 ;
    ecal =0.;

    float rcal,sglcal,nsigcal,ksigcal,rcali;

    // version determines which constants are used for which year
    float  rcem=RCEM_92;;
    float  rcha=RCHA_92;
    float  rwha=RWHA_92;
    float  rcemi=RCEMI_92;
    float  rchai=RCHAI_92;
    float  rwhai=RWHAI_92;
    float  atcemm = exp(-RCEM_92*0.5);
    float  atcham = exp(-RCHA_92*0.5);
    float  atwham = exp(-RWHA_92*0.5);

    float  sglcem  = SGLCEM;
    float  signcem = SIGNCEM;
    float  sigkcem = SIGKCEM;
    float  sglcha  = SGLCHA;
    float  signcha = SIGNCHA;
    float  sigkcha = SIGKCHA;
    float  sglwha  = SGLWHA;
    float  signwha = SIGNWHA;
    float  sigkwha = SIGKWHA;
    if (version==1) {
      rcem=RCEM_01;
      rcha=RCHA;
      rwha=RWHA;
      rcemi=RCEMI_01;
      rchai=RCHAI;
      rwhai=RWHAI;
      atcemm = exp(-RCEM*0.5);
      atcham = exp(-RCHA*0.5);
      atwham = exp(-RWHA*0.5);
    } else {
      rcem=RCEM;
      rcha=RCHA;
      rwha=RWHA;
      rcemi=RCEMI;
      rchai=RCHAI;
      rwhai=RWHAI;
    }
    atcemm = exp(-RCEM*0.5);
    atcham = exp(-RCHA*0.5);
    atwham = exp(-RWHA*0.5);

    switch (detector) {
    case 1:
      sglcal =sglcem;
      nsigcal=signcem;       //number of sigma 
      ksigcal=sigkcem;       //width of sigma 
      rcal   =rcem;
      rcali  =rcemi;
      break;
    case 2:
      sglcal =sglcha;
      nsigcal=signcha;
      ksigcal=sigkcha;
      rcal   =rcha;
      rcali  =rchai;
      break;
    case 3:
      sglcal =sglwha;
      nsigcal=signwha;
      ksigcal=sigkwha;
      rcal   =rwha;
      rcali  =rwhai;
      break;
    }


    // if a tower has a dead/hot PMT (defined in DB tables ***BadChannels)
    // then ignore the dead/hot PMT (should have the lowest/highest energy) by setting its energy to the 
    // energy of the other. Dead/hot towers are stored in validcode
    if (badpmt0) ecal0 = ecal1;
    if (badpmt1) ecal1 = ecal0;
    if (badpmt0 && badpmt1) {ecal0 = 0.; ecal1 = 0.;}



    if( ecal0<=0. || ecal1 <= 0.) {
      float ex = ecal0;
      if (ecal1>ecal0) {
	ex = ecal1;
      }
      
      phical = 0.5;
      if(ex < sglcal) {
	ecal = .5*ex;
      }
      else {
	if(killSpike) {
	  flag=1;
	}
      }
    } 
    
    else { 
      
      double lrat = log((double)ecal0/ecal1);
      // simple run2 spike killer
      if(killSpike == 1 ) {

	//the cut here roughly corresponds to 1000 counts
	if(ecal0+ecal1 > 3.0){
                float asym = fabs(ecal0-ecal1)/(ecal0+ecal1);
                if(asym > 0.95) flag = 1;
	}
	phical = 0.5-0.5*rcali*lrat; 
	ecal   = sqrt(ecal0*ecal1);
      }
      // fancy run1 spike killer
      if(killSpike == 2 ) {
	float cut = rcal + nsigcal*ksigcal/sqrt(.5*(ecal0+ecal1));
	if(fabs(lrat) > cut) {
	  if (ecal0<ecal1) {
	    ecal = ecal0;
	  }
	  else {   
	    ecal = ecal1;
	  }
	  
	  phical = .5;
	  flag=1;
	}
	else {
	  phical = 0.5-0.5*rcali*lrat; 
	  ecal   = sqrt(ecal0*ecal1);
	}
      }
      if(killSpike < 1){
	phical = 0.5 - 0.5*rcali*lrat; 
	ecal   = sqrt(ecal0*ecal1);
      }
      
      if (1.<phical) {
	phical = 1.;
      }
    if (0.>phical){
      phical = 0.;
    }
    
    }
}
// } // namespace calor