//
//    File: DFCALShower_factory.cc
// Created: Tue May 17 11:57:50 EST 2005
// Creator: remitche (on Linux mantrid00 2.4.20-18.8smp i686)
// Edited: B. Schaefer 3/23/2012 removed radiation hard insert functionality

#include <math.h>
#include <DVector3.h>
//#include <DLorentzVector.h>
using namespace std;

#include "DFCALShower_factory.h"
#include "DFCALGeometry.h"
//#include "DFCALCluster.h"
#include "DFCALHit.h"
#include <JANA/JEvent.h>
#include <JANA/JApplication.h>
using namespace jana;

//----------------
// Constructor
//----------------
DFCALShower_factory::DFCALShower_factory()
{

// Set of coefficients for non-linear energy corrections 

  SHOWER_ENERGY_THRESHOLD = 50*k_MeV;

  /*
    NON_LIN_COEF_A[0] = 0.439287;
  NON_LIN_COEF_B[0] = 0.503378; 
  NON_LIN_COEF_C[0] = 1.80842;
  NON_LIN_COEF_alfa[0] = 1+0.0724789;
  */
  NON_LIN_COEF_A[0] = 0.7218;
  NON_LIN_COEF_B[0] = -11.96;
  NON_LIN_COEF_C[0] = -1.461;
  NON_LIN_COEF_alfa[0] = -0.513;
 
  NON_LIN_COEF_A[1] = 0.5879;
  NON_LIN_COEF_B[1] = 0.30737; 
  NON_LIN_COEF_C[1] = 2.1461;
  NON_LIN_COEF_alfa[1] = 1.0;
  

  //Loose timing cut
  SHOWER_TIMING_WINDOW = 5.0;

// Parameters to make shower-depth correction taken from Radphi, 
// slightly modifed to match photon-polar angle
        FCAL_RADIATION_LENGTH[0] = 3.1;
        FCAL_CRITICAL_ENERGY[0] = 0.035;
        FCAL_SHOWER_OFFSET[0] = 1.0;

	FCAL_RADIATION_LENGTH[1] = 0.89;
        FCAL_CRITICAL_ENERGY[1] = 0.00964;
        FCAL_SHOWER_OFFSET[1] = 1.0;


	/*
	gPARMS->SetDefaultParameter("FCAL:SHOWER_ENERGY_THRESHOLD", SHOWER_ENERGY_THRESHOLD);

	gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_A", NON_LIN_COEF_A);
	if (NON_LIN_COEF_A<=0.) {
	  cout  << "Warning: DFCALShower : parameter A=" <<  NON_LIN_COEF_A 
		<< " is not valid!" << endl; 
	}


	gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_B", NON_LIN_COEF_B);
	gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_C", NON_LIN_COEF_C);
	gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_alfa", NON_LIN_COEF_alfa);

	gPARMS->SetDefaultParameter("FCAL:SHOWER_TIMING_WINDOW", SHOWER_TIMING_WINDOW);
	
	gPARMS->SetDefaultParameter("FCAL:FCAL_RADIATION_LENGTH", FCAL_RADIATION_LENGTH);
	gPARMS->SetDefaultParameter("FCAL:FCAL_CRITICAL_ENERGY", FCAL_CRITICAL_ENERGY);
	gPARMS->SetDefaultParameter("FCAL:FCAL_SHOWER_OFFSET", FCAL_SHOWER_OFFSET);

	*/
}

//------------------
// brun
//------------------
// take merging out
jerror_t DFCALShower_factory::brun(JEventLoop *loop, int runnumber)
{
  /*// Get calibration constants
	map<string, double> cluster_merging;
	loop->GetCalib("FCAL/cluster_merging", cluster_merging);
	if(cluster_merging.find("MIN_CLUSTER_SEPARATION")!=cluster_merging.end()){
		MIN_CLUSTER_SEPARATION = cluster_merging["MIN_CLUSTER_SEPARATION"];
		if(debug_level>0)jout<<"MIN_CLUSTER_SEPARATION = "<<MIN_CLUSTER_SEPARATION<<endl;
	}else{
		jerr<<"Unable to get from MIN_CLUSTER_SEPARATION FCAL/cluster_merging in Calib database!"<<endl;
		loop->GetJApplication()->Quit();
        }*/

  // Get calibration constants
  FCAL_C_EFFECTIVE[0]=FCAL_C_EFFECTIVE[1]=15.;
  map<string, double> fcal_parms;
  loop->GetCalib("FCAL/fcal_parms", fcal_parms);
  if (fcal_parms.find("FCAL_C_EFFECTIVE")!=fcal_parms.end()){
    FCAL_C_EFFECTIVE[0] = fcal_parms["FCAL_C_EFFECTIVE"];
    FCAL_C_EFFECTIVE[1] =  FCAL_C_EFFECTIVE[0];
    if(debug_level>0)jout<<"FCAL_C_EFFECTIVE = "<<FCAL_C_EFFECTIVE<<endl;
  } else {
    jerr<<"Unable to get FCAL_C_EFFECTIVE from FCAL/fcal_parms in Calib database!"<<endl;
  }
  
  m_zTarget=65.0*k_cm;
  m_FCALfront[0]= DFCALGeometry::fcalFaceZ();
  m_FCALback[0]= m_FCALfront[0]+DFCALGeometry::blockLength(0);
  m_FCALfront[1]=638.8;
  m_FCALback[1]=656.8;
  DApplication *dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
  if (dapp) {
    const DGeometry *geom = dapp->GetDGeometry(runnumber);
    geom->GetTargetZ(m_zTarget);
    geom->GetFCALZ(m_FCALfront[0]);

    vector<double>block;
    geom->Get("//box[@name='LGBL']/@X_Y_Z",block);
    m_FCALback[0]=m_FCALfront[0]+block[2];

    geom->Get("//box[@name='LTBK']/@X_Y_Z",block);
    m_FCALfront[1]=0.5*(m_FCALback[0]+m_FCALfront[0]-block[2]);
    m_FCALback[1]=m_FCALfront[1]+block[2];


    printf("back %f front %f\n",m_FCALback[0],m_FCALfront[0]);
    printf("back %f front %f\n",m_FCALback[1],m_FCALfront[1]);
  }

  return NOERROR;
}


//------------------
// evnt
//------------------
jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, int eventnumber)
{
  vector<const DFCALCluster*> fcalClusters;
  eventLoop->Get(fcalClusters);
  if(fcalClusters.size()<1)return NOERROR;
 
  // Use the center of the target as an approximation for the vertex position
  DVector3 target(0.0, 0.0, m_zTarget);

  // Loop over list of DFCALCluster objects and calculate the "Non-linear" corrected
  // energy and position for each. We'll use a logarithmic energy-weighting to 
  // find the final position and error. 
  for( vector< const DFCALCluster* >::const_iterator clItr = fcalClusters.begin();
       clItr != fcalClusters.end();  ++clItr ){
    const DFCALCluster* cluster=*clItr;

    double cTime = cluster->getTime();
 		
    double errX = cluster->getRMS_x();
    double errY = cluster->getRMS_y();
    double errZ;  // will be filled by call to GetCorrectedEnergyAndPosition()
		
    // Get corrected energy, position, and errZ
    double Ecorrected;
    DVector3 pos_corrected;
    GetCorrectedEnergyAndPosition( cluster , Ecorrected, pos_corrected, errZ, &target);

    if (Ecorrected>0.){		
      double x=pos_corrected.X();
      double y=pos_corrected.Y();
      unsigned int index=0;
      if (fabs(x)<60.5 && fabs(y)<60.5) index=1;

      //up to this point, all times have been times at which light reaches
      //the back of the detector. Here we correct for the time that it 
      //takes the Cherenkov light to reach the back of the detector
      //so that the t reported is roughly the time of the shower at the
      //position pos_corrected	
      cTime -= ( m_FCALback[index] - pos_corrected.Z() )/FCAL_C_EFFECTIVE[index];

      // Make the DFCALShower object
      DFCALShower* shower = new DFCALShower;
      
      shower->setEnergy( Ecorrected );
      shower->setPosition( pos_corrected );   
      shower->setPosError( errX, errY, errZ );
      shower->setTime ( cTime );
      
      shower->AddAssociatedObject(cluster);

      _data.push_back(shower);
    }
  }

  return NOERROR;
}

//--------------------------------
// GetCorrectedEnergyAndPosition
//
// Non-linear and depth corrections should be fixed within DFCALShower member functions
//--------------------------------
void DFCALShower_factory::GetCorrectedEnergyAndPosition(const DFCALCluster* cluster, double &Ecorrected, DVector3 &pos_corrected, double &errZ, const DVector3 *vertex)
{
  // Non-linear energy correction are done here
  int MAXITER = 1000;

  DVector3  posInCal = cluster->getCentroid();
  float x0 = posInCal.Px();
  float y0 = posInCal.Py();

  double Eclust = cluster->getEnergy();
  
  unsigned int index=0;
  if (fabs(x0)<60.5 && fabs(y0)<60.5) index=1;

  double A  = NON_LIN_COEF_A[index];
  double B  = NON_LIN_COEF_B[index];
  double C  = NON_LIN_COEF_C[index];
  double alfa  = NON_LIN_COEF_alfa[index];
	 
  double Egamma = 0.;

  if (A > 0 ) { 
    
    Egamma = Eclust/A;

    for ( int niter=0; 1; niter++) {

      double energy = Egamma;
      double non_lin_part = pow(Egamma,1+alfa)/(B+C*Egamma);
      Egamma = Eclust/A - non_lin_part;
      if ( fabs( (Egamma-energy)/energy ) < 0.001 ) {
	break;
	
      }
      else if ( niter > MAXITER ) {
	
	cout << " Iteration failed for cluster energy " << Eclust << endl;
	Egamma  = 0;
        
	break;
	
      }
      
    }
    
  }
  else {
    //cout  << "Warning: DFCALShower : parameter A=" <<  NON_LIN_COEF_A 
    //<< " is not valid!" << endl; 
    Egamma=Eclust;
  }

  // then depth corrections 
  if ( Egamma > 0 ) { 
    float dxV = x0-vertex->X();
    float dyV = y0-vertex->Y();
    float zV = vertex->Z();
   
    double z0 = m_FCALfront[index] - zV;
    double zMax = FCAL_RADIATION_LENGTH[index]*(FCAL_SHOWER_OFFSET[index] 
					 + log(Egamma/FCAL_CRITICAL_ENERGY[index]));
    double zed = z0;
    double zed1 = z0 + zMax;

    double r0 = sqrt(dxV*dxV + dyV*dyV );

    int niter;
    for ( niter=0; niter<100; niter++) {
      double tt = r0/zed1;
      zed = z0 + zMax/sqrt( 1 + tt*tt );
      if ( fabs( (zed-zed1) ) < 0.001) {
	break;
      }
      zed1 = zed;
    }
    
    posInCal.SetZ( zed + zV );
    errZ = zed - zed1;
  }
  
  Ecorrected = Egamma;
  //Ecorrected=pow(Eclust/0.8348,0.989);
  //Ecorrected=Eclust;
  //posInCal.SetZ(m_FCALfront);
  pos_corrected = posInCal;
}