source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMicroManager.cc @ 1337

//
// ********************************************************************
// * 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-04-beta-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;
}