//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: G4StatMFMicroManager.cc,v 1.6 2008/07/25 11:20:47 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-03-cand-01 $
//
// Hadronic Process: Nuclear De-excitations
// by V. Lara


#include "G4StatMFMicroManager.hh"
#include "G4HadronicException.hh"


// Copy constructor
G4StatMFMicroManager::G4StatMFMicroManager(const G4StatMFMicroManager & )
{
    throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroManager::copy_constructor meant to not be accessable");
}

// Operators

G4StatMFMicroManager & G4StatMFMicroManager::
operator=(const G4StatMFMicroManager & )
{
    throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroManager::operator= meant to not be accessable");
    return *this;
}


G4bool G4StatMFMicroManager::operator==(const G4StatMFMicroManager & ) const
{
    return false;
}
 

G4bool G4StatMFMicroManager::operator!=(const G4StatMFMicroManager & ) const
{
    return true;
}



// constructor
G4StatMFMicroManager::G4StatMFMicroManager(const G4Fragment & theFragment, const G4int multiplicity,
					   const G4double FreeIntE, const G4double SCompNuc) : 
    _Normalization(0.0)
{
    // Perform class initialization
    Initialize(theFragment,multiplicity,FreeIntE,SCompNuc);
}


// destructor
G4StatMFMicroManager::~G4StatMFMicroManager() 
{
  if (!_Partition.empty()) 
    {
      std::for_each(_Partition.begin(),_Partition.end(),
		      DeleteFragment());
    }
}



// Initialization method

void G4StatMFMicroManager::Initialize(const G4Fragment & theFragment, const G4int m, 
				      const G4double FreeIntE, const G4double SCompNuc) 
{
    G4int i;

    G4double U = theFragment.GetExcitationEnergy();

    G4double A = theFragment.GetA();
    G4double Z = theFragment.GetZ();
	
    // Statistical weights
    _WW = 0.0;

    // Mean breakup multiplicity
    _MeanMultiplicity = 0.0;

    // Mean channel temperature
    _MeanTemperature = 0.0;

    // Mean channel entropy
    _MeanEntropy = 0.0;	
	
    // Keep fragment atomic numbers
// 	G4int * FragmentAtomicNumbers = new G4int(static_cast<G4int>(A+0.5));
//	G4int * FragmentAtomicNumbers = new G4int(m);
    G4int FragmentAtomicNumbers[4];
	
    // We distribute A nucleons between m fragments mantaining the order
    // FragmentAtomicNumbers[m-1]>FragmentAtomicNumbers[m-2]>...>FragmentAtomicNumbers[0]
    // Our initial distribution is 
    // FragmentAtomicNumbers[m-1]=A, FragmentAtomicNumbers[m-2]=0, ..., FragmentAtomicNumbers[0]=0
    FragmentAtomicNumbers[m-1] = static_cast<G4int>(A);
    for (i = 0; i <  (m - 1); i++) FragmentAtomicNumbers[i] = 0;

    // We try to distribute A nucleons in partitions of m fragments
    // MakePartition return true if it is possible 
    // and false if it is not	
    while (MakePartition(m,FragmentAtomicNumbers)) {
	// Allowed partitions are stored and its probability calculated
			
	G4StatMFMicroPartition * aPartition = new G4StatMFMicroPartition(static_cast<G4int>(A),
									 static_cast<G4int>(Z));
	G4double PartitionProbability = 0.0;
			
	for (i = m-1; i >= 0; i--) aPartition->SetPartitionFragment(FragmentAtomicNumbers[i]);
	PartitionProbability = aPartition->CalcPartitionProbability(U,FreeIntE,SCompNuc);
	_Partition.push_back(aPartition);
			
	_WW += PartitionProbability;
	_MeanMultiplicity += m*PartitionProbability;
	_MeanTemperature += aPartition->GetTemperature() * PartitionProbability;
	if (PartitionProbability > 0.0) 
	    _MeanEntropy += PartitionProbability * aPartition->GetEntropy();
			
    }
		
	
    // garbage collection
// 	delete [] FragmentAtomicNumbers;
	
}


G4bool G4StatMFMicroManager::MakePartition(const G4int k, G4int * ANumbers)
    // Distributes A nucleons between k fragments
    // mantaining the order ANumbers[k-1] > ANumbers[k-2] > ... > ANumbers[0]
    // If it is possible returns true. In other case returns false
{
    G4int l = 1;
    while (l < k) {
	G4int tmp = ANumbers[l-1] + ANumbers[k-1];
	ANumbers[l-1] += 1;
	ANumbers[k-1] -= 1;
	if (ANumbers[l-1] > ANumbers[l] || ANumbers[k-2] > ANumbers[k-1]) {
	    ANumbers[l-1] = 1;
	    ANumbers[k-1] = tmp - 1;
	    l++;
	} else return true;
    }
    return false;
}



void G4StatMFMicroManager::Normalize(const G4double Norm)
{
    _Normalization = Norm;
    _WW /= Norm;
    _MeanMultiplicity /= Norm;
    _MeanTemperature /= Norm;
    _MeanEntropy /= Norm; 
	
    return;
}

G4StatMFChannel * G4StatMFMicroManager::ChooseChannel(const G4double A0, const G4double Z0, 
						      const G4double MeanT)
{
    G4double RandNumber = _Normalization * _WW * G4UniformRand();
    G4double AccumWeight = 0.0;
	
    for (std::vector<G4StatMFMicroPartition*>::iterator i = _Partition.begin();
	 i != _Partition.end(); ++i)
    {
	AccumWeight += (*i)->GetProbability();
	if (RandNumber < AccumWeight)
	    return (*i)->ChooseZ(A0,Z0,MeanT);
    }

    throw G4HadronicException(__FILE__, __LINE__, 
	"G4StatMFMicroCanonical::ChooseChannel: Couldn't find a channel.");
    return 0;
}