G4StatMFMicroManager.cc

Go to the documentation of this file.

//
// ********************************************************************
// * 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 67983 2013-03-13 10:42:03Z gcosmo $
//
// 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, G4int multiplicity,
                       G4double FreeIntE, 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, G4int im, 
                      G4double FreeIntE, G4double SCompNuc) 
{
    G4int i;

    G4double U = theFragment.GetExcitationEnergy();

    G4int A = theFragment.GetA_asInt();
    G4int Z = theFragment.GetZ_asInt();
    
    // 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[im-1] = A;
    for (i = 0; i <  (im - 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(im,FragmentAtomicNumbers)) {
    // Allowed partitions are stored and its probability calculated
            
    G4StatMFMicroPartition * aPartition = new G4StatMFMicroPartition(A,Z);
    G4double PartitionProbability = 0.0;
            
    for (i = im-1; i >= 0; i--) aPartition->SetPartitionFragment(FragmentAtomicNumbers[i]);
    PartitionProbability = aPartition->CalcPartitionProbability(U,FreeIntE,SCompNuc);
    _Partition.push_back(aPartition);
            
    _WW += PartitionProbability;
    _MeanMultiplicity += im*PartitionProbability;
    _MeanTemperature += aPartition->GetTemperature() * PartitionProbability;
    if (PartitionProbability > 0.0) 
        _MeanEntropy += PartitionProbability * aPartition->GetEntropy();
            
    }
}

G4bool G4StatMFMicroManager::MakePartition(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(G4double Norm)
{
    _Normalization = Norm;
    _WW /= Norm;
    _MeanMultiplicity /= Norm;
    _MeanTemperature /= Norm;
    _MeanEntropy /= Norm; 
    
    return;
}

G4StatMFChannel* 
G4StatMFMicroManager::ChooseChannel(G4int A0, G4int Z0, 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;
}