diff --git a/geant4/LEMuSR/G4Modified/G4VDecayChannel.cc b/geant4/LEMuSR/G4Modified/G4VDecayChannel.cc new file mode 100644 index 0000000..8475c50 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4VDecayChannel.cc @@ -0,0 +1,460 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * 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. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4VDecayChannel.cc,v 1.16 2004/12/02 08:09:00 kurasige Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-03 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 27 July 1996 H.Kurashige +// 30 May 1997 H.Kurashige +// 23 Mar. 2000 H.Weber : add GetAngularMomentum +// ------------------------------------------------------------ + +#include "G4ParticleDefinition.hh" +#include "G4ParticleTable.hh" +#include "G4DecayTable.hh" +#include "G4DecayProducts.hh" +#include "G4VDecayChannel.hh" + +const G4String G4VDecayChannel::noName = " "; + +G4VDecayChannel::G4VDecayChannel(const G4String &aName, G4int Verbose) + :kinematics_name(aName), + rbranch(0.0), + numberOfDaughters(0), + parent_name(0), daughters_name(0), + particletable(0), + parent(0), daughters(0), + parent_mass(0.0), daughters_mass(0), + verboseLevel(Verbose) +{ + // set pointer to G4ParticleTable (static and singleton object) + particletable = G4ParticleTable::GetParticleTable(); +} + +G4VDecayChannel::G4VDecayChannel(const G4String &aName, + const G4String& theParentName, + G4double theBR, + G4int theNumberOfDaughters, + const G4String& theDaughterName1, + const G4String& theDaughterName2, + const G4String& theDaughterName3, + const G4String& theDaughterName4 ) + :kinematics_name(aName), + rbranch(theBR), + numberOfDaughters(theNumberOfDaughters), + parent_name(0), daughters_name(0), + particletable(0), + parent(0), daughters(0), + parent_mass(0.0), daughters_mass(0), + verboseLevel(1) +{ + // set pointer to G4ParticleTable (static and singleton object) + particletable = G4ParticleTable::GetParticleTable(); + + // parent name + parent_name = new G4String(theParentName); + + // cleate array + daughters_name = new G4String*[numberOfDaughters]; + for (G4int index=0;index0) daughters_name[0] = new G4String(theDaughterName1); + if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2); + if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3); + if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4); +} + + + +G4VDecayChannel::G4VDecayChannel(const G4VDecayChannel &right) +{ + kinematics_name = right.kinematics_name; + verboseLevel = right.verboseLevel; + rbranch = right.rbranch; + + // copy parent name + parent_name = new G4String(*right.parent_name); + parent = 0; + parent_mass = 0.0; + + //create array + numberOfDaughters = right.numberOfDaughters; + + if ( numberOfDaughters >0 ) { + daughters_name = new G4String*[numberOfDaughters]; + //copy daughters name + for (G4int index=0; index < numberOfDaughters; index++) + { + daughters_name[index] = new G4String(*right.daughters_name[index]); + } + } + + // + daughters_mass = 0; + daughters = 0; + + // particle table + particletable = G4ParticleTable::GetParticleTable(); +} + +G4VDecayChannel & G4VDecayChannel::operator=(const G4VDecayChannel &right) +{ + if (this != &right) { + kinematics_name = right.kinematics_name; + verboseLevel = right.verboseLevel; + rbranch = right.rbranch; + + // copy parent name + parent_name = new G4String(*right.parent_name); + + // clear daughters_name array + ClearDaughtersName(); + + // recreate array + numberOfDaughters = right.numberOfDaughters; + if ( numberOfDaughters >0 ) { + daughters_name = new G4String*[numberOfDaughters]; + //copy daughters name + for (G4int index=0; index < numberOfDaughters; index++) { + daughters_name[index] = new G4String(*right.daughters_name[index]); + } + } + } + + // + parent = 0; + daughters = 0; + parent_mass = 0.0; + daughters_mass = 0; + + // particle table + particletable = G4ParticleTable::GetParticleTable(); + + return *this; +} + + +G4VDecayChannel::~G4VDecayChannel() +{ + if (parent_name != 0) delete parent_name; + ClearDaughtersName(); + if (daughters_mass != 0) delete [] daughters_mass; +} + +void G4VDecayChannel::ClearDaughtersName() +{ + if ( daughters_name != 0) { + if (numberOfDaughters>0) { +#ifdef G4VERBOSE + if (verboseLevel>1) { + G4cout << "G4VDecayChannel::ClearDaughtersName "; + G4cout << "clear all daughters " << G4endl; + } +#endif + for (G4int index=0; index < numberOfDaughters; index++) { + if (daughters_name[index] != 0) delete daughters_name[index]; + } + } + delete [] daughters_name; + daughters_name = 0; + } + // + if (daughters != 0) delete [] daughters; + if (daughters_mass != 0) delete [] daughters_mass; + daughters = 0; + daughters_mass = 0; + + numberOfDaughters = 0; +} + +void G4VDecayChannel::SetNumberOfDaughters(G4int size) +{ + if (size >0) { + // remove old contents + ClearDaughtersName(); + // cleate array + daughters_name = new G4String*[size]; + for (G4int index=0;index0) { + G4cout << "G4VDecayChannel::SetDaughter: "; + G4cout << "Number of daughters is not defined" << G4endl; + } +#endif + return; + } + + // check existence of daughters_name array + if (daughters_name == 0) { + // cleate array + daughters_name = new G4String*[numberOfDaughters]; + for (G4int index=0;index=numberOfDaughters) ) { +#ifdef G4VERBOSE + if (verboseLevel>0) { + G4cout << "G4VDecayChannel::SetDaughter"; + G4cout << "index out of range " << anIndex << G4endl; + } +#endif + } else { + // delete the old name if it exists + if (daughters_name[anIndex]!=0) delete daughters_name[anIndex]; + // fill the name + daughters_name[anIndex] = new G4String(particle_name); + // refill the array of daughters[] if it exists + if (daughters != 0) FillDaughters(); +#ifdef G4VERBOSE + if (verboseLevel>1) { + G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :"; + G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<GetParticleName()); +} + +void G4VDecayChannel::FillDaughters() +{ + G4int index; + +#ifdef G4VERBOSE + if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <GetPDGMass(); + + // + G4double sumofdaughtermass = 0.0; + if ((numberOfDaughters <=0) || (daughters_name == 0) ){ +#ifdef G4VERBOSE + if (verboseLevel>0) { + G4cout << "G4VDecayChannel::FillDaughters "; + G4cout << "[ " << parent->GetParticleName() << " ]"; + G4cout << "numberOfDaughters is not defined yet"; + } +#endif + daughters = 0; + G4Exception("G4VDecayChannel::FillDaughters"); + } + + //create and set the array of pointers to daughter particles + daughters = new G4ParticleDefinition*[numberOfDaughters]; + if (daughters_mass != 0) delete [] daughters_mass; + daughters_mass = new G4double[numberOfDaughters]; + // loop over all daughters + for (index=0; index < numberOfDaughters; index++) { + if (daughters_name[index] == 0) { + // daughter name is not defined +#ifdef G4VERBOSE + if (verboseLevel>0) { + G4cout << "G4VDecayChannel::FillDaughters "; + G4cout << "[ " << parent->GetParticleName() << " ]"; + G4cout << index << "-th daughter is not defined yet" << G4endl; + } +#endif + daughters[index] = 0; + G4Exception("G4VDecayChannel::FillDaughters"); + } + //search daughter particles in the particle table + daughters[index] = particletable->FindParticle(*daughters_name[index]); + if (daughters[index] == 0) { + // can not find the daughter particle +#ifdef G4VERBOSE + if (verboseLevel>0) { + G4cout << "G4VDecayChannel::FillDaughters "; + G4cout << "[ " << parent->GetParticleName() << " ]"; + G4cout << index << ":" << *daughters_name[index]; + G4cout << " is not defined !!" << G4endl; + G4cout << " The BR of this decay mode is set to zero " << G4endl; + } +#endif + SetBR(0.0); + return; + } +#ifdef G4VERBOSE + if (verboseLevel>1) { + G4cout << index << ":" << *daughters_name[index]; + G4cout << ":" << daughters[index] << G4endl; + } +#endif + daughters_mass[index] = daughters[index]->GetPDGMass(); + sumofdaughtermass += daughters[index]->GetPDGMass(); + } // end loop over all daughters + + // check sum of daghter mass + G4double widthMass = parent->GetPDGWidth(); + if ( (parent->GetParticleType() != "nucleus") && + (sumofdaughtermass > parentmass + 5*widthMass) ){ + // !!! illegal mass !!! +#ifdef G4VERBOSE + if (GetVerboseLevel()>0) { + G4cout << "G4VDecayChannel::FillDaughters "; + G4cout << "[ " << parent->GetParticleName() << " ]"; + G4cout << " Energy/Momentum conserevation breaks " <1) { + G4cout << " parent:" << *parent_name; + G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <GetPDGMass()/GeV; + G4cout << "[GeV/c/c]" <0) { + G4cout << "G4VDecayChannel::FillParent "; + G4cout << ": parent name is not defined !!" << G4endl; + } +#endif + parent = 0; + G4Exception("G4VDecayChannel::FillParent"); + } + // search parent particle in the particle table + parent = particletable->FindParticle(*parent_name); + if (parent == 0) { + // parent particle does not exist +#ifdef G4VERBOSE + if (verboseLevel>0) { + G4cout << "G4VDecayChannel::FillParent "; + G4cout << *parent_name << " does not exist !!" << G4endl; + } +#endif + G4Exception("G4VDecayChannel::FillParent"); + } + parent_mass = parent->GetPDGMass(); +} + +void G4VDecayChannel::SetParent(const G4ParticleDefinition * parent_type) +{ + if (parent_type != 0) SetParent(parent_type->GetParticleName()); +} + +G4int G4VDecayChannel::GetAngularMomentum() +{ + // determine angular momentum + + // fill pointers to daughter particles if not yet set + if (daughters == 0) FillDaughters(); + + const G4int PiSpin = parent->GetPDGiSpin(); + const G4int PParity = parent->GetPDGiParity(); + if (2==numberOfDaughters) { // up to now we can only handle two particle decays + const G4int D1iSpin = daughters[0]->GetPDGiSpin(); + const G4int D1Parity = daughters[0]->GetPDGiParity(); + const G4int D2iSpin = daughters[1]->GetPDGiSpin(); + const G4int D2Parity = daughters[1]->GetPDGiParity(); + const G4int MiniSpin = std::abs (D1iSpin - D2iSpin); + const G4int MaxiSpin = D1iSpin + D2iSpin; + const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int + G4int lMin; +#ifdef G4VERBOSE + if (verboseLevel>1) { + G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl; + G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl; + } +#endif + for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings + lMin = std::abs(PiSpin-j)/2; +#ifdef G4VERBOSE + if (verboseLevel>1) + G4cout << "-> checking 2*j=" << j << G4endl; +#endif + for (G4int l=lMin; l<=lMax; l++) { +#ifdef G4VERBOSE + if (verboseLevel>1) + G4cout << " checking l=" << l << G4endl; +#endif + if (l%2==0) { + if (PParity == D1Parity*D2Parity) { // check parity for this l + return l; + } + } else { + if (PParity == -1*D1Parity*D2Parity) { // check parity for this l + return l; + } + } + } + } + } else { + G4Exception ("G4VDecayChannel::GetAngularMomentum: Sorry, can't handle 3 particle decays (up to now)"); + } + G4Exception ("G4VDecayChannel::GetAngularMomentum: Can't find angular momentum for this decay!"); + return 0; +} + +void G4VDecayChannel::DumpInfo() +{ + G4cout << " BR: " << rbranch << " [" << kinematics_name << "]"; + G4cout << " : " ; + for (G4int index=0; index < numberOfDaughters; index++) + { + if(daughters_name[index] != 0) { + G4cout << " " << *(daughters_name[index]); + } else { + G4cout << " not defined "; + } + } + G4cout << G4endl; +} + +const G4String& G4VDecayChannel::GetNoName() const +{ + return noName; +} diff --git a/geant4/LEMuSR/G4Modified/G4VDecayChannel.hh b/geant4/LEMuSR/G4Modified/G4VDecayChannel.hh new file mode 100644 index 0000000..01e18c2 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4VDecayChannel.hh @@ -0,0 +1,288 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * 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. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4VDecayChannel.hh,v 1.11 2004/12/02 08:08:58 kurasige Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-03 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 27 July 1996 H.Kurashige +// 30 May 1997 H.Kurashige +// 23 Mar. 2000 H.Weber : add GetAngularMomentum() +// ------------------------------------------------------------ +#ifndef G4VDecayChannel_h +#define G4VDecayChannel_h 1 + +#include "G4ios.hh" +#include "globals.hh" +#include +#include "G4ThreeVector.hh" + +class G4ParticleDefinition; +class G4DecayProducts; +class G4ParticleTable; + +class G4VDecayChannel +{ + // Class Description + // This class is a abstract class to describe decay kinematics + // + + public: + //Constructors + G4VDecayChannel(const G4String &aName, G4int Verbose = 1); + G4VDecayChannel(const G4String &aName, + const G4String& theParentName, + G4double theBR, + G4int theNumberOfDaughters, + const G4String& theDaughterName1, + const G4String& theDaughterName2 = "", + const G4String& theDaughterName3 = "", + const G4String& theDaughterName4 = "" ); + + // Destructor + virtual ~G4VDecayChannel(); + + private: + // copy constructor and assignment operatotr + G4VDecayChannel(const G4VDecayChannel &); + G4VDecayChannel & operator=(const G4VDecayChannel &); + + public: + // equality operators + G4int operator==(const G4VDecayChannel &right) const {return (this == &right);} + G4int operator!=(const G4VDecayChannel &right) const {return (this != &right);} + + // less-than operator is defined for G4DecayTable + G4int operator<(const G4VDecayChannel &right) const; + + public: // With Description + virtual G4DecayProducts* DecayIt(G4double parentMass = -1.0) = 0; + + public: // With Description + //get kinematics name + G4String GetKinematicsName() const; + //get branching ratio + G4double GetBR() const; + //get number of daughter particles + G4int GetNumberOfDaughters() const; + + //get the pointer to the parent particle + G4ParticleDefinition * GetParent(); + //get the pointer to a daughter particle + G4ParticleDefinition * GetDaughter(G4int anIndex); + + //get the angular momentum of the decay + G4int GetAngularMomentum(); + //get the name of the parent particle + const G4String& GetParentName() const; + //get the name of a daughter particle + const G4String& GetDaughterName(G4int anIndex) const; + + // get mass of parent + G4double GetParentMass() const; + G4double GetDaughterMass(G4int anIndex) const; + G4ThreeVector GetParentPolarization(); + //set the parent particle (by name or by pointer) + void SetParent(const G4ParticleDefinition * particle_type); + void SetParent(const G4String &particle_name); + //set branching ratio + void SetBR(G4double value); + //set number of daughter particles + void SetNumberOfDaughters(G4int value); + //set a daughter particle (by name or by pointer) + void SetDaughter(G4int anIndex, + const G4ParticleDefinition * particle_type); + void SetDaughter(G4int anIndex, + const G4String &particle_name); + void SetParentPolarization(G4ThreeVector polar); + protected: + // kinematics name + G4String kinematics_name; + // branching ratio [0.0 - 1.0] + G4double rbranch; + // number of daughters + G4int numberOfDaughters; + // parent particle + G4String* parent_name; + //daughter particles + G4String** daughters_name; + + protected: // With Description + // celar daughters array + void ClearDaughtersName(); + + protected: + // pointer to particle table + G4ParticleTable* particletable; + + // temporary buffers of pointers to G4ParticleDefinition + G4ParticleDefinition* parent; + G4ParticleDefinition** daughters; + + // parent mass + G4double parent_mass; + G4double* daughters_mass; + + //parent theParentPolarization + G4ThreeVector theParentPolarization; + + // fill daughters array + void FillDaughters(); + // fill parent + void FillParent(); + + public: // With Description + void SetVerboseLevel(G4int value); + G4int GetVerboseLevel() const; + void DumpInfo(); + + private: + const G4String& GetNoName() const; + + private: + // controle flag for output message + G4int verboseLevel; + // 0: Silent + // 1: Warning message + // 2: More + + static const G4String noName; +}; + +inline + G4int G4VDecayChannel::operator<(const G4VDecayChannel &right) const +{ + return (this->rbranch < right.rbranch); +} + +inline + G4ParticleDefinition* G4VDecayChannel::GetDaughter(G4int anIndex) + { + //pointers to daughter particles are filled, if they are not set yet + if (daughters == 0) FillDaughters(); + + //get the pointer to a daughter particle + if ( (anIndex>=0) && (anIndex0) + G4cout << "G4VDecayChannel::GetDaughter index out of range "<=0) && (anIndex0){ + G4cout << "G4VDecayChannel::GetDaughterName "; + G4cout << "index out of range " << anIndex << G4endl; + } + return GetNoName(); + } +} + +inline + G4double G4VDecayChannel::GetDaughterMass(G4int anIndex) const +{ + if ( (anIndex>=0) && (anIndex0){ + G4cout << "G4VDecayChannel::GetDaughterMass "; + G4cout << "index out of range " << anIndex << G4endl; + } + return 0.0; + } +} + +inline + G4ParticleDefinition* G4VDecayChannel::GetParent() +{ + //the pointer to the parent particle is filled, if it is not set yet + if (parent == 0) FillParent(); + //get the pointer to the parent particle + return parent; +} + +inline + const G4String& G4VDecayChannel::GetParentName() const +{ + return *parent_name; +} + +inline + G4double G4VDecayChannel::GetParentMass() const +{ + return parent_mass; +} + + +inline + void G4VDecayChannel::SetParent(const G4String &particle_name) +{ + if (parent_name != 0) delete parent_name; + parent_name = new G4String(particle_name); + parent = 0; +} + +inline + G4int G4VDecayChannel::GetNumberOfDaughters() const +{ + return numberOfDaughters; +} + +inline + G4String G4VDecayChannel::GetKinematicsName() const { return kinematics_name; } + +inline + void G4VDecayChannel::SetBR(G4double value){ rbranch = value; } + +inline + G4double G4VDecayChannel::GetBR() const { return rbranch; } + +inline + void G4VDecayChannel::SetVerboseLevel(G4int value){ verboseLevel = value; } + +inline + G4int G4VDecayChannel::GetVerboseLevel() const { return verboseLevel; } + + +inline G4ThreeVector G4VDecayChannel::GetParentPolarization(){return theParentPolarization;} +inline void G4VDecayChannel::SetParentPolarization(G4ThreeVector polar){theParentPolarization=polar;} + + +#endif + + + + + +