//	BLGlobalField.cc
/*
This source file is part of G4beamline, http://g4beamline.muonsinc.com
Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved.

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

http://www.gnu.org/copyleft/gpl.html
*/

#include <time.h>

#include "Randomize.hh"
#include "G4EqMagElectricField.hh"
#include "G4TransportationManager.hh"
#include "G4FieldManager.hh"
#include "G4PropagatorInField.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4ChordFinder.hh"
#include "G4ClassicalRK4.hh"
#include "G4SimpleRunge.hh"
#include "G4RunManager.hh"
#include "G4Event.hh"

#include "BLGlobalField.hh"
#include "BLParam.hh"

// define parameters:
BLSetParam minStep("minStep","0.01","Minimum step size (mm)");
BLSetParam deltaChord("deltaChord","3.0","Geant4 tracking parameter");
BLSetParam deltaOneStep("deltaOneStep","0.01","Geant4 tracking parameter");
BLSetParam deltaIntersection("deltaIntersection","0.1","Geant4 tracking parameter");
BLSetParam epsMin("epsMin","2.5e-7","Geant4 tracking parameter");
BLSetParam epsMax("epsMax","0.05","Geant4 tracking parameter");

bool BLGlobalField::spinTracking = false;
G4EqEMFieldWithSpin *BLGlobalField::spinEqRhs = 0;
BLGlobalField *BLGlobalField::object = 0;

BLGlobalField::BLGlobalField() : G4ElectroMagneticField(), fields()
{
	first = true;
	fp = 0;
	nfp = 0;
	forceEfieldZero = false;

	clear();

	// Get transportation, field, and propagator  managers
	G4TransportationManager  *pTransportMgr= 
               G4TransportationManager::GetTransportationManager();
	G4FieldManager* fieldMgr = pTransportMgr->GetFieldManager();
	G4PropagatorInField *pFieldPropagator =
					pTransportMgr->GetPropagatorInField(); 

	// Need to SetFieldChangesEnergy to account for a time varying electric
	// field (r.f. fields)
	fieldMgr->SetFieldChangesEnergy(true);

	fieldMgr->SetDetectorField(this);

	// Construct equation of motion of particles through e.m. fields
	G4MagIntegratorStepper *pStepper = 0;
	if(spinTracking) {
		spinEqRhs = new G4EqEMFieldWithSpin(this);
		G4EquationOfMotion *equation = spinEqRhs;
		pStepper = new G4ClassicalRK4(equation,12);
	} else {
		spinEqRhs = 0;
		G4EquationOfMotion *equation = new G4EqMagElectricField(this);
		pStepper = new G4ClassicalRK4(equation,8);
	}
	G4double min_step = Param.getDouble("minStep");
	G4MagInt_Driver *fIntgrDriver = new G4MagInt_Driver(min_step, 
				pStepper, pStepper->GetNumberOfVariables() ); 
	G4ChordFinder *pChordFinder = new G4ChordFinder(fIntgrDriver); 
	fieldMgr->SetChordFinder(pChordFinder); 

	// Set accuracy parameters
	G4double delta_chord = Param.getDouble("deltaChord");
	pChordFinder->SetDeltaChord( delta_chord );
	G4double delta_onestep = Param.getDouble("deltaOneStep");
	fieldMgr->SetAccuraciesWithDeltaOneStep(delta_onestep);
	G4double delta_intersection = Param.getDouble("deltaIntersection");
	fieldMgr->SetDeltaIntersection(delta_intersection); 
	G4double eps_min = Param.getDouble("epsMin");
	G4double eps_max = Param.getDouble("epsMax");
	pFieldPropagator->SetMinimumEpsilonStep(eps_min);
	pFieldPropagator->SetMaximumEpsilonStep(eps_max);
//	G4cout << "Accuracy Parameters:" <<
//    		" MinStep=" << min_step <<
//   		" DeltaChord=" << delta_chord <<
//  		" DeltaOneStep=" << delta_onestep << G4endl;
//	G4cout << "                    " <<
//		" DeltaIntersection=" << delta_intersection <<
//    		" EpsMin=" << eps_min <<
//   		" EpsMax=" << eps_max <<  G4endl;
	fieldMgr->SetChordFinder(pChordFinder);

	// set object
	object = this;
}

BLGlobalField::~BLGlobalField()
{
	clear();
}

BLGlobalField *BLGlobalField::getObject()
{
	if(!object) new BLGlobalField();
	return object;
}

void BLGlobalField::GetFieldValue(const G4double *point, G4double *field) const
{
	// NOTE: this routine dominates the CPU time for tracking.
	// The test case is the MICE beamline + detector, which has 38 
	// entries in fields[].
	// Using bounding boxes sped it up by a factor of 8 (13 ev/sec => 
	// 100 ev/sec).
	// Using the simple array fp[] instead of fields[] directly sped
	// it up by an additional factor of two (200 ev/sec).
	// Making this a friend of BLElementField and putting the code for
	// isInBoundingBox() inline here didn't help at all.

	field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
	// protect against Geant4 bug that calls us with point[] NaN.
	if(point[0] != point[0]) return;

	if(first)		// (can't use nfp or fp, as they may change)
		((BLGlobalField*)this)->setupArray();	// (cast away const)

	for(int i=0; i<nfp; ++i) {
		const BLElementField *p = fp[i];
		if(p->isInBoundingBox(point))
			p->addFieldValue(point,field);
	}

	if(forceEfieldZero) {
		field[3] = field[4] = field[5] = 0.0;
	}
}

void BLGlobalField::clear()
{
	std::vector<const BLElementField *>::iterator i;
	for(i=fields.begin(); i!=fields.end(); ++i) {
		delete *i;
	}
	fields.clear();

	if(fp) delete fp;
	first = true;
	fp = 0;
	nfp = 0;
}

void BLGlobalField::setupArray()
{
	first = false;
	nfp = fields.size();
	fp = new const BLElementField *[nfp+1]; // add 1 so it's never 0
	for(int i=0; i<nfp; ++i)
		fp[i] = fields[i];
}

void BLGlobalField::setSpinTracking(bool v)
{
	if(object) G4Exception("BLGlobalField","Invalid setSpinTracking() call",
				FatalException,"");
	spinTracking = v;
}