Geant4: G4ChiralInvariantPhaseSpace Class Reference

G4HadFinalState * G4ChiralInvariantPhaseSpace::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	aTargetNucleus,
		G4HadFinalState *	aChange = `0`
	)

Definition at line 52 of file G4ChiralInvariantPhaseSpace.cc.

References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4ParticleTable::FindIon(), G4ParticleTable::FindParticle(), G4QEnvironment::Fragment(), G4cerr, G4cout, G4endl, G4Gamma::Gamma(), G4QCHIPSWorld::Get(), G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4IonTable::GetIonMass(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4QCHIPSWorld::GetParticles(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGEncoding(), G4Nucleus::GetZ_asInt(), isAlive, G4Neutron::Neutron(), G4DynamicParticle::SetDefinition(), G4HadFinalState::SetEnergyChange(), G4DynamicParticle::SetMomentum(), G4HadFinalState::SetMomentumChange(), G4Quasmon::SetParameters(), G4QNucleus::SetParameters(), G4HadFinalState::SetStatusChange(), and stopAndKill.

Referenced by G4StringChipsParticleLevelInterface::ApplyYourself(), G4StringChipsInterface::ApplyYourself(), G4QStringChipsParticleLevelInterface::ApplyYourself(), G4GammaNuclearReaction::ApplyYourself(), G4ElectroNuclearReaction::ApplyYourself(), G4ProtonAntiProtonAtRestChips::AtRestDoIt(), and G4PionMinusNuclearAtRestChips::AtRestDoIt().

00056 {
00057   G4HadFinalState * aResult;
00058   if(aChange != 0)
00059   {
00060     aResult = aChange;
00061   }
00062   else
00063   {
00064     aResult = & theResult;
00065     aResult->Clear();
00066     aResult->SetStatusChange(stopAndKill);
00067   }
00068   //projectile properties needed in constructor of quasmon
00069   G4LorentzVector proj4Mom;
00070   proj4Mom = aTrack.Get4Momentum();
00071   G4int projectilePDGCode = aTrack.GetDefinition()
00072                                   ->GetPDGEncoding();
00073   
00074   //target properties needed in constructor of quasmon
00075   G4int targetZ = aTargetNucleus.GetZ_asInt();
00076   G4int targetA = aTargetNucleus.GetA_asInt();
00077   G4int targetPDGCode = 90000000 + 1000*targetZ + (targetA-targetZ);
00078   // NOT NECESSARY ______________
00079   G4double targetMass = G4ParticleTable::GetParticleTable()->GetIonTable()
00080                                                            ->GetIonMass(targetZ, targetA);
00081   G4LorentzVector targ4Mom(0.,0.,0.,targetMass);  
00082   // END OF NOT NECESSARY^^^^^^^^
00083 
00084   // V.Ivanchenko set the same value as for other CHIPS models
00085   G4int nop = 152; 
00086   //G4int nop = 164; // nuclear clusters up to A=21
00087   G4double fractionOfSingleQuasiFreeNucleons = 0.4;
00088   G4double fractionOfPairedQuasiFreeNucleons = 0.0;
00089   if(targetA>27) fractionOfPairedQuasiFreeNucleons = 0.04;
00090   G4double clusteringCoefficient = 4.;
00091   G4double temperature = 180.;
00092   G4double halfTheStrangenessOfSee = 0.1; // = s/d = s/u
00093   G4double etaToEtaPrime = 0.3;
00094   
00095   // construct and fragment the quasmon
00096   //G4QCHIPSWorld aWorld(nop);              // Create CHIPS World of nop particles
00097   G4QCHIPSWorld::Get()->GetParticles(nop);  // Create CHIPS World of nop particles
00098   G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
00099                             fractionOfPairedQuasiFreeNucleons,
00100                             clusteringCoefficient);
00101   G4Quasmon::SetParameters(temperature,
00102                            halfTheStrangenessOfSee,
00103                            etaToEtaPrime);
00104   //  G4QEnvironment::SetParameters(solidAngle);
00105 //  G4cout << "Input info "<< projectilePDGCode << " " 
00106 //         << targetPDGCode <<" "
00107 //       << 1./MeV*proj4Mom<<" "
00108 //       << 1./MeV*targ4Mom << " "
00109 //       << nop << G4endl;
00110   G4QHadronVector projHV;
00111   G4QHadron* iH = new G4QHadron(projectilePDGCode, 1./MeV*proj4Mom);
00112   projHV.push_back(iH);
00113 
00114   //AND
00115   //A. Dotti 24 Aug. : Trying to handle situation when G4QEnvironment::Fragment throws an exception
00116   //                   Seen by ATLAS for gamma+Nuclear interaction for Gamma@2.4GeV on Al
00117   //                   The poor-man solution here is to re-try interaction if a G4QException is catched
00118   //                   Warning: G4QExcpetion does NOT inherit from base class G4HadException
00119   G4QHadronVector* output=0;
00120 
00121   bool retry=true;
00122   int retryes=0;
00123   int maxretries=10;
00124 
00125   while ( retry && retryes < maxretries ) 
00126     {
00127       G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
00128 
00129       try
00130         {
00131           ++retryes;
00132           output = pan->Fragment();
00133           retry=false;//If here, Fragment did not throw! (AND)
00134         }
00135       catch( ... )
00136         {
00137           G4cerr << "***WARNING*** Exception thrown passing through G4ChiralInvariantPhaseSpace "<<G4endl;
00138           G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
00139           G4cerr << " Dumping the information in the pojectile list"<<G4endl;
00140           for(size_t i=0; i< projHV.size(); i++)
00141             {
00142               G4cerr <<"  Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
00143                      <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
00144             }
00145           G4cerr << "Retrying interaction "<<G4endl; //AND
00146           //throw; //AND
00147         }
00148       delete pan;
00149     } //AND
00150   if ( retryes >= maxretries ) //AND
00151     {
00152       G4cerr << "***ERROR*** Maximum number of retries ("<<maxretries<<") reached for G4QEnvironment::Fragment(), exception is being thrown" << G4endl;
00153       throw G4HadronicException(__FILE__,__LINE__,"G4ChiralInvariantPhaseSpace::ApplyYourself(...) - Maximum number of re-tries reached, abandon interaction");
00154     }
00155   
00156   // Fill the particle change.
00157   G4DynamicParticle * theSec;
00158 #ifdef CHIPSdebug
00159   G4cout << "G4ChiralInvariantPhaseSpace: NEW EVENT #ofHadrons="
00160          <<output->size()<<G4endl;
00161 #endif
00162 
00163 
00164   unsigned int particle;
00165   for( particle = 0; particle < output->size(); particle++)
00166   {
00167     if(output->operator[](particle)->GetNFragments() != 0) 
00168     {
00169       delete output->operator[](particle);
00170       continue;
00171     }
00172     theSec = new G4DynamicParticle;  
00173     G4int pdgCode = output->operator[](particle)->GetPDGCode();
00174 #ifdef CHIPSdebug
00175     G4cout << "G4ChiralInvariantPhaseSpace: h#"<<particle
00176            <<", PDG="<<pdgCode<<G4endl;
00177 #endif
00178     G4ParticleDefinition * theDefinition;
00179     // Note that I still have to take care of strange nuclei
00180     // For this I need the mass calculation, and a changed interface
00181     // for ion-tablel ==> work for Hisaya @@@@@@@
00182     // Then I can sort out the pdgCode. I also need a decau process 
00183     // for strange nuclei; may be another chips interface
00184     if(pdgCode>90000000) 
00185     {
00186       G4int aZ = (pdgCode-90000000)/1000;
00187       G4int anN = pdgCode-90000000-1000*aZ;
00188       theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
00189       if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
00190     }
00191     //AND
00192     else if ( pdgCode == 90000000 && output->operator[](particle)->Get4Momentum().m()<1*MeV )
00193       {
00194         //A. Dotti: 
00195         //We cure the case the model returns a (A,Z)=0,0 G4QHadron with a very small mass
00196         //We convert this to a gamma. According to the author of the model this is the 
00197         //correct thing to do and it is done also in other parts of the CHIPS package
00198         theDefinition = G4Gamma::Gamma();
00199       }
00200     //AND
00201     else
00202     {
00203       theDefinition = G4ParticleTable::GetParticleTable()
00204         ->FindParticle(output->operator[](particle)->GetPDGCode());
00205     }
00206 
00207     //AND
00208     //A. Dotti: Handle problematic cases in which one of the products has not been recognized
00209     //          This should never happen but we want to add an extra-protection
00210     if ( theDefinition == NULL )
00211       {
00212         //If we arrived here something bad really happened. We do not know how to handle the products of the interaction, we give up, resetting the 
00213         //result and keeping the primary alive.
00214         G4cerr<<"**WARNING*** G4ChiralInvariantPhaseSpace::ApplyYourself(...) : G4QEnvironment::Fragment() returns an invalid fragment\n with fourMom(MeV)=";
00215         G4cerr<<output->operator[](particle)->Get4Momentum()<<" and mass(MeV)="<<output->operator[](particle)->Get4Momentum().m();
00216         G4cerr<<". Offending PDG is:"<<pdgCode<<" abandon interaction. \n Taget PDG was:"<<targetPDGCode<<" \n Dumping the information in the projectile list:\n";
00217         for(size_t i=0; i< projHV.size(); i++)
00218           {
00219             G4cerr <<" Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
00220                  <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<"\n";
00221           }
00222         G4cerr<<"\n Please report as bug \n***END OF MESSAGE***"<<G4endl;
00223 
00224         for ( unsigned int cparticle=0 ; cparticle<output->size();++cparticle)
00225           delete output->operator[](cparticle);
00226         delete output;
00227         std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00228         projHV.clear();
00229 
00230         aResult->Clear();
00231         aResult->SetStatusChange(isAlive);
00232         aResult->SetEnergyChange(aTrack.GetKineticEnergy());
00233         aResult->SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00234         return aResult;
00235       }
00236     //AND
00237 
00238 
00239     theSec->SetDefinition(theDefinition);
00240     theSec->SetMomentum(output->operator[](particle)->Get4Momentum().vect());
00241     aResult->AddSecondary(theSec); 
00242     delete output->operator[](particle);
00243   }
00244   delete output;
00245   std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00246   projHV.clear();
00247   return aResult;
00248 }


Public Member Functions
	G4ChiralInvariantPhaseSpace ()
	~G4ChiralInvariantPhaseSpace ()
G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus, G4HadFinalState *aChange=0)

G4ChiralInvariantPhaseSpace Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation