source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/include/G4QNucleus.hh @ 836

//
// ********************************************************************
// * 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: G4QNucleus.hh,v 1.33 2007/10/31 13:23:07 mkossov Exp $
// GEANT4 tag $Name:  $
//
//      ---------------- G4QNucleus ----------------
//             by Mikhail Kossov, Sept 1999.
//  class header for Quasmon initiated Candidates used by the CHIPS Model
// ----------------------------------------------------------------------

#ifndef G4QNucleus_h
#define G4QNucleus_h 1

#include "G4QCandidateVector.hh"
#include "G4QHadronVector.hh"
#include "G4QChipolino.hh"
#include <utility>
#include <vector>
#include "globals.hh"

class G4QNucleus : public G4QHadron
{
public:
  G4QNucleus();                                            // Default Constructor
  G4QNucleus(G4int nucPDG);                                // At Rest PDG-Constructor
  G4QNucleus(G4LorentzVector p, G4int nucPDG);             // Full PDG-Constructor
  G4QNucleus(G4QContent nucQC);                            // At Rest QuarkCont-Constructor
  G4QNucleus(G4QContent nucQC, G4LorentzVector p);         // Full QuarkCont-Constructor
  G4QNucleus(G4int z, G4int n, G4int s=0);                 // At Rest ZNS-Constructor
  G4QNucleus(G4int z, G4int n, G4int s, G4LorentzVector p);// Full ZNS-Constructor
  //G4QNucleus(const G4QNucleus& right);                     // Copy Constructor by value
  G4QNucleus(G4QNucleus* right);                           // Copy Constructor by pointer
  ~G4QNucleus();                                           // Public Destructor
  // Overloaded Operators
  const G4QNucleus& operator=(const G4QNucleus& right);
  G4bool operator==(const G4QNucleus &right) const {return this==&right;}
  G4bool operator!=(const G4QNucleus &right) const {return this!=&right;}
  // Specific Selectors
  G4int      GetPDG()      const {return 90000000+1000*(1000*S+Z)+N;}// PDG Code of Nucleus
  G4int      GetZ()        const {return Z;}                                // Get a#of protons
  G4int      GetN()        const {return N;}                                // Get a#of neutrons
  G4int      GetS()        const {return S;}                                // Get a#of lambdas
  G4int      GetA()        const {return Z+N+S;}    // Get A of the nucleus
  G4int      GetDZ()       const {return dZ;}                                            // Get a#of protons in dense region
  G4int      GetDN()       const {return dN;}                                            // Get a#of neutrons in dense region
  G4int      GetDS()       const {return dS;}                                            // Get a#of lambdas in dense region
  G4int      GetDA()       const {return dZ+dN+dS;} // Get A of the dense part of nucleus
  G4int      GetMaxClust() const {return maxClust;} // Get Max BarNum of Clusters
  G4double   GetProbability(G4int bn=0) const {return probVect[bn];} // clust(BarN)probabil
  G4double   GetMZNS()     const {return GetQPDG().GetNuclMass(Z,N,S);} // not H or Q
  G4double   GetGSMass()   const {return GetQPDG().GetMass();}//Nucleus GSMass (not Hadron)
  G4QContent GetQCZNS()    const                   // Get ZNS quark content of Nucleus
  {
    if(S>=0) return G4QContent(Z+N+N+S,Z+Z+N+S,S,0,0,0);
    else     return G4QContent(Z+N+N+S,Z+Z+N+S,0,0,0,-S);
  }
  G4int      GetNDefMesonC() const{return nDefMesonC;}; // max#of predefed mesonCandidates
  G4int      GetNDefBaryonC()const{return nDefBaryonC;};// max#of predefed baryonCandidates
  std::pair<G4double, G4double> RefetchImpactXandY() const {return theImpactParameter;}
  G4double GetDensity(const G4ThreeVector&aPos) {return rho0*GetRelativeDensity(aPos);}
  G4double GetRelativeDensity(const G4ThreeVector& aPosition); // Densyty/rho0
  G4double GetRadius(const G4double maxRelativeDenisty=0.5); // Radius of %ofDensity
  G4double GetOuterRadius();                        // Get radius of the most far nucleon
  G4double GetDeriv(const G4ThreeVector& point);    // Derivitive of density
  G4double GetFermiMomentum(G4double density);      // Returns modul of FermyMomentum(dens)
  G4ThreeVector Get3DFermiMomentum(G4double density, G4double maxMom=-1.)
  {
    if(maxMom<0) maxMom=GetFermiMomentum(density);
    return maxMom*RandomUnitSphere();
  }
  G4QHadron* GetNextNucleon()
    {return (currentNucleon>=0&&currentNucleon<GetA()) ? theNucleons[currentNucleon++] :0;}
  //std::vector<G4double>* GetBThickness() const {return Tb;} // T(b) function, step .1 fm
  std::vector<G4double> const* GetBThickness() {return &Tb;} // T(b) function, step .1 fm

  // Specific Modifiers
  G4bool     EvaporateBaryon(G4QHadron* h1,G4QHadron* h2); // Evaporate Baryon from Nucleus
  void       EvaporateNucleus(G4QHadron* hA, G4QHadronVector* oHV);// Evaporate Nucleus
  //void DecayBaryon(G4QHadron* dB, G4QHadronVector* oHV); // gamma+N or Delt->N+Pi @@later
  void       DecayDibaryon(G4QHadron* dB, G4QHadronVector* oHV);   // deuteron is kept
  void       DecayIsonucleus(G4QHadron* dB, G4QHadronVector* oHV); // nP+(Pi+) or nN+(Pi-)
  void       DecayMultyBaryon(G4QHadron* dB, G4QHadronVector* oHV);// A*p, A*n or A*L
  void       DecayAntiStrange(G4QHadron* dB, G4QHadronVector* oHV);// nuclei with K+/K0
  void       DecayAlphaBar(G4QHadron* dB, G4QHadronVector* oHV);   // alpha+p or alpha+n
  void       DecayAlphaDiN(G4QHadron* dB, G4QHadronVector* oHV);   // alpha+p+p
  void       DecayAlphaAlpha(G4QHadron* dB, G4QHadronVector* oHV); // alpha+alpha
  G4int      SplitBaryon();                         // Is it possible to split baryon/alpha
  G4int      HadrToNucPDG(G4int hPDG);              // Converts hadronic PDGCode to nuclear
  G4int      NucToHadrPDG(G4int nPDG);              // Converts nuclear PDGCode to hadronic
  G4bool     Split2Baryons();                       // Is it possible to split two baryons?
  void       ActivateBThickness();                  // Calculate T(b) for nucleus (db=.1fm)
  void       InitByPDG(G4int newPDG);               // Init existing nucleus by new PDG
  void       InitByQC(G4QContent newQC)             // Init existing nucleus by new QCont
                                {G4int PDG=G4QPDGCode(newQC).GetPDGCode(); InitByPDG(PDG);}
  void       IncProbability(G4int bn);              // Add one cluster to probability
  void       Increase(G4int PDG, G4LorentzVector LV = G4LorentzVector(0.,0.,0.,0.));
  void       Increase(G4QContent QC, G4LorentzVector LV = G4LorentzVector(0.,0.,0.,0.));
  void       Reduce(G4int PDG);                     // Reduce Nucleus by PDG fragment
  void       CalculateMass() {Set4Momentum(G4LorentzVector(0.,0.,0.,GetGSMass()));}
  void       SetMaxClust(G4int maxC){maxClust=maxC;}// Set Max BarNum of Clusters
  void       InitCandidateVector(G4QCandidateVector& theQCandidates,
                                 G4int nM=45, G4int nB=72, G4int nC=117);
  void       PrepareCandidates(G4QCandidateVector& theQCandidates, G4bool piF=false, G4bool
                               gaF=false, G4LorentzVector LV=G4LorentzVector(0.,0.,0.,0.));
  G4int      UpdateClusters(G4bool din);            // Return a#of clusters & calc.probab's
  G4QNucleus operator+=(const G4QNucleus& rhs);     // Add a cluster to the  nucleus
  G4QNucleus operator-=(const G4QNucleus& rhs);     // Subtract a cluster from a nucleus
  G4QNucleus operator*=(const G4int& rhs);          // Multiplication of the Nucleus
  G4bool StartLoop();                               // returns size of theNucleons (cN=0)
  G4bool ReduceSum(G4ThreeVector* momentum, G4double*); // Reduce momentum nonconservation
  void DoLorentzBoost(const G4LorentzVector& theBoost); // Boost nucleons by 4-vector
  void DoLorentzBoost(const G4ThreeVector& theBeta);// Boost nucleons by v/c
  void DoLorentzContraction(const G4LorentzVector&B){DoLorentzContraction(B.vect()/B.e());}
  void DoLorentzContraction(const G4ThreeVector& theBeta); // Lorentz Contraction by v/c
  void DoTranslation(const G4ThreeVector& theShift); // Used only in GHAD-TFT

  // Static functions
  static void SetParameters(G4double fN=.1,G4double fD=.05, G4double cP=4., G4double mR=1.,
                            G4double nD=.8*fermi);

  // Specific General Functions
  G4ThreeVector RandomUnitSphere();                 // Randomize position inside UnitSphere
  G4int RandomizeBinom(G4double p,G4int N);         // Randomize according to Binomial Law
  G4double CoulombBarrier(const G4double& cZ=1, const G4double& cA=1, G4double dZ=0.,
                          G4double dA=0.);          // CoulombBarrier in MeV
  G4double FissionCoulombBarrier(const G4double& cZ, const G4double& cA, G4double dZ=0.,
                          G4double dA=0.);          // Fission CoulombBarrier in MeV
  G4double BindingEnergy(const G4double& cZ=0, const G4double& cA=0, G4double dZ=0.,
                         G4double dA=0.);
  G4double CoulBarPenProb(const G4double& CB, const G4double& E, const G4int& C,
                          const G4int& B);
  std::pair<G4double, G4double> ChooseImpactXandY(G4double maxImpact); // Randomize bbar
  void ChooseNucleons();                            // Initializes 3D Nucleons
  void ChoosePositions();                           // Initializes positions of 3D nucleons
  void ChooseFermiMomenta();                        // Initializes FermyMoms of 3D nucleons
  void InitDensity();                               // Initializes density distribution
  void Init3D();                                    // automatically starts the LOOP
private:  
  // Specific Encapsulated Functions
  void       SetZNSQC(G4int z, G4int n, G4int s);   // Set QC, using Z,N,S
  G4QNucleus GetThis() const {return G4QNucleus(Z,N,S);} // @@ Check for memory leak

// Body
private:
  // Static Parameters
  static const G4int nDefMesonC =45;
  static const G4int nDefBaryonC=72;
  //
  static G4double freeNuc;        // probability of the quasi-free baryon on surface
  static G4double freeDib;        // probability of the quasi-free dibaryon on surface
  static G4double clustProb;      // clusterization probability in dense region
  static G4double mediRatio;      // relative vacuum hadronization probability
  static G4double nucleonDistance;// Distance between nucleons (0.8 fm)
  // The basic  
  G4int Z;                        // Z of the Nucleus
  G4int N;                        // N of the Nucleus
  G4int S;                        // S of the Nucleus
  // The secondaries
  G4int dZ;                       // Z of the dense region of the nucleus
  G4int dN;                       // N of the dense region of the nucleus
  G4int dS;                       // S of the dense region of the nucleus
  G4int maxClust;                 // Baryon Number of the last calculated cluster
  G4double probVect[256];         // Cluster probability ("a#of issues" can be real) Vector
  // 3D
  std::pair<G4double, G4double> theImpactParameter; // 2D impact parameter vector bbar
  G4QHadronVector theNucleons;  // Vector of nucleons of which Nucleus consists of
  G4int currentNucleon;         // Current nucleon for the NextNucleon (? M.K.)
  G4double rho0;                // Normalazation density
  G4double radius;              // Nuclear radius
  //std::vector<G4double>* Tb;    // T(b) function with step .1 fm (@@ make .1 a parameter)
  std::vector<G4double> Tb;    // T(b) function with step .1 fm (@@ make .1 a parameter)
};

std::ostream& operator<<(std::ostream& lhs, G4QNucleus& rhs);
std::ostream& operator<<(std::ostream& lhs, const G4QNucleus& rhs);


#endif