Geant4: G4RPGTwoCluster Class Reference

G4bool G4RPGTwoCluster::ReactionStage	(	const G4HadProjectile *	,
		G4ReactionProduct &	,
		G4bool &	,
		const G4DynamicParticle *	,
		G4ReactionProduct &	,
		G4bool &	,
		const G4Nucleus &	,
		G4ReactionProduct &	,
		G4FastVector< G4ReactionProduct, 256 > &	,
		G4int &	,
		G4bool	,
		G4ReactionProduct &
	)
Reimplemented from G4RPGReaction.
Definition at line 44 of file G4RPGTwoCluster.cc.
Referenced by G4RPGInelastic::CalculateMomenta().
00056 {
00057   // Derived from H. Fesefeldt's FORTRAN code TWOCLU
00058   //
00059   // A simple two cluster model is used to generate x- and pt- values for 
00060   // incident, target, and all secondary particles. 
00061   // This should be sufficient for low energy interactions.
00062 
00063   G4int i;
00064   G4ParticleDefinition* aProton = G4Proton::Proton();
00065   G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
00066   G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus();
00067   G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus();
00068   G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
00069   G4bool veryForward = false;
00070 
00071   const G4double protonMass = aProton->GetPDGMass()/MeV;
00072   const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
00073   const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
00074   const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
00075   const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
00076   G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
00077   G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
00078                                           targetMass*targetMass +
00079                                           2.0*targetMass*etOriginal);  // GeV
00080   G4double currentMass = currentParticle.GetMass()/GeV;
00081   targetMass = targetParticle.GetMass()/GeV;
00082 
00083   if (currentMass == 0.0 && targetMass == 0.0) {
00084     G4double ek = currentParticle.GetKineticEnergy();
00085     G4ThreeVector mom = currentParticle.GetMomentum();
00086     currentParticle = *vec[0];
00087     targetParticle = *vec[1];
00088     for (i = 0; i < (vecLen-2); ++i) *vec[i] = *vec[i+2];
00089     if (vecLen < 2) {
00090       for (G4int j = 0; j < vecLen; j++) delete vec[j];
00091       vecLen = 0;
00092       throw G4HadReentrentException(__FILE__, __LINE__,
00093       "G4RPGTwoCluster::ReactionStage : Negative number of particles");
00094     }
00095     delete vec[vecLen-1];
00096     delete vec[vecLen-2];
00097     vecLen -= 2;
00098     currentMass = currentParticle.GetMass()/GeV;
00099     targetMass = targetParticle.GetMass()/GeV;
00100     incidentHasChanged = true;
00101     targetHasChanged = true;
00102     currentParticle.SetKineticEnergy(ek);
00103     currentParticle.SetMomentum(mom);
00104     veryForward = true;
00105   }
00106 
00107   const G4double atomicWeight = targetNucleus.GetA_asInt();
00108   const G4double atomicNumber = targetNucleus.GetZ_asInt();
00109 
00110   // particles have been distributed in forward and backward hemispheres
00111   // in center of mass system of the hadron nucleon interaction
00112 
00113   // Incident particle always in forward hemisphere
00114 
00115   G4int forwardCount = 1;        // number of particles in forward hemisphere
00116   currentParticle.SetSide( 1 );
00117   G4double forwardMass = currentParticle.GetMass()/GeV;
00118   G4double cMass = forwardMass;
00119     
00120   // Target particle always in backward hemisphere
00121   G4int backwardCount = 1;       // number of particles in backward hemisphere
00122   targetParticle.SetSide( -1 );
00123   G4double backwardMass = targetParticle.GetMass()/GeV;
00124   G4double bMass = backwardMass;
00125 
00126   //  G4int backwardNucleonCount = 1;  // number of nucleons in backward hemisphere
00127   for (i = 0; i < vecLen; ++i) {
00128     if (vec[i]->GetSide() < 0) vec[i]->SetSide(-1);   // to take care of 
00129     // case where vec has been preprocessed by GenerateXandPt
00130     // and some of them have been set to -2 or -3
00131     if (vec[i]->GetSide() == -1) {
00132       ++backwardCount;
00133       backwardMass += vec[i]->GetMass()/GeV;
00134     } else {
00135       ++forwardCount;
00136       forwardMass += vec[i]->GetMass()/GeV;
00137     }
00138   }
00139 
00140   // Add nucleons and some pions from intra-nuclear cascade
00141   G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
00142   if(term1 < 0) term1 = 0.0001; // making sure xtarg<0;
00143   const G4double afc = 0.312 + 0.2 * std::log(term1);
00144   G4double xtarg;
00145   if( centerofmassEnergy < 2.0+G4UniformRand() )        // added +2 below, JLC 4Jul97
00146     xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
00147   else
00148     xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
00149   if( xtarg <= 0.0 )xtarg = 0.01;
00150   G4int nuclearExcitationCount = G4Poisson( xtarg );
00151 
00152   if(atomicWeight<1.0001) nuclearExcitationCount = 0;
00153   //  G4int extraNucleonCount = 0;
00154   //  G4double extraMass = 0.0;
00155   //  G4double extraNucleonMass = 0.0;
00156   if( nuclearExcitationCount > 0 )
00157   {
00158     G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );     
00159     const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
00160     //
00161     //  NOTE: in TWOCLU, these new particles were given negative codes
00162     //        here we use  NewlyAdded = true  instead
00163     //
00164     for( i=0; i<nuclearExcitationCount; ++i )
00165     {
00166       G4ReactionProduct* pVec = new G4ReactionProduct();
00167       if( G4UniformRand() < nucsup[momentumBin] )  // add proton or neutron
00168       {
00169         if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
00170           pVec->SetDefinition( aProton );
00171         else
00172           pVec->SetDefinition( aNeutron );
00173         // Not used        ++backwardNucleonCount;
00174         // Not used        ++extraNucleonCount;
00175         // Not used        extraNucleonMass += pVec->GetMass()/GeV;
00176       }
00177       else
00178       {                                       // add a pion
00179         G4double ran = G4UniformRand();
00180         if( ran < 0.3181 )
00181           pVec->SetDefinition( aPiPlus );
00182         else if( ran < 0.6819 )
00183           pVec->SetDefinition( aPiZero );
00184         else
00185           pVec->SetDefinition( aPiMinus );
00186 
00187         // DHW: add following two lines to correct energy balance
00188         //        ++backwardCount;
00189         //        backwardMass += pVec->GetMass()/GeV;
00190       }
00191       pVec->SetSide( -2 );    // backside particle
00192       // Not used     extraMass += pVec->GetMass()/GeV;
00193       pVec->SetNewlyAdded( true );
00194       vec.SetElement( vecLen++, pVec );
00195     }
00196   }
00197 
00198   // Masses of particles added from cascade not included in energy balance.
00199   // That's correct for nucleons from the intra-nuclear cascade but not for 
00200   // pions from the cascade.
00201  
00202   G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
00203   G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
00204   G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
00205   G4bool secondaryDeleted;
00206   G4double pMass;
00207 
00208   while( eAvailable <= 0.0 )   // must eliminate a particle
00209   {
00210     secondaryDeleted = false;
00211     for( i=(vecLen-1); i>=0; --i )
00212     {
00213       if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
00214       {
00215         pMass = vec[i]->GetMass()/GeV;
00216         for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];     // shift up
00217         --forwardCount;
00218         forwardEnergy += pMass;
00219         forwardMass -= pMass;
00220         secondaryDeleted = true;
00221         break;
00222       }
00223       else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
00224       {
00225         pMass = vec[i]->GetMass()/GeV;
00226         for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
00227         --backwardCount;
00228         backwardEnergy += pMass;
00229         backwardMass -= pMass;
00230         secondaryDeleted = true;
00231         break;
00232       }
00233     }  // breaks go down to here
00234 
00235     if( secondaryDeleted )
00236     {
00237       delete vec[vecLen-1];
00238       --vecLen;
00239         // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00240     }
00241     else
00242     {
00243       if( vecLen == 0 ) return false;  // all secondaries have been eliminated
00244       if( targetParticle.GetSide() == -1 )
00245       {
00246         pMass = targetParticle.GetMass()/GeV;
00247         targetParticle = *vec[0];
00248         for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
00249         --backwardCount;
00250         backwardEnergy += pMass;
00251         backwardMass -= pMass;
00252         secondaryDeleted = true;
00253       }
00254       else if( targetParticle.GetSide() == 1 )
00255       {
00256         pMass = targetParticle.GetMass()/GeV;
00257         targetParticle = *vec[0];
00258         for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
00259         --forwardCount;
00260         forwardEnergy += pMass;
00261         forwardMass -= pMass;
00262         secondaryDeleted = true;
00263       }
00264 
00265       if( secondaryDeleted )
00266       {
00267         delete vec[vecLen-1];
00268         --vecLen;
00269       }
00270       else
00271       {
00272         if( currentParticle.GetSide() == -1 )
00273         {
00274           pMass = currentParticle.GetMass()/GeV;
00275           currentParticle = *vec[0];
00276           for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
00277           --backwardCount;
00278           backwardEnergy += pMass;
00279           backwardMass -= pMass;
00280           secondaryDeleted = true;
00281         }
00282         else if( currentParticle.GetSide() == 1 )
00283         {
00284           pMass = currentParticle.GetMass()/GeV;
00285           currentParticle = *vec[0];
00286           for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
00287           --forwardCount;
00288           forwardEnergy += pMass;
00289           forwardMass -= pMass;
00290           secondaryDeleted = true;
00291         }
00292         if( secondaryDeleted )
00293         {
00294           delete vec[vecLen-1];
00295           --vecLen;
00296         }
00297         else break;
00298 
00299       }  // secondary not deleted 
00300     }  // secondary not deleted
00301 
00302     eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
00303   } // while
00304 
00305   //
00306   // This is the start of the TwoCluster function
00307   // Choose multi-particle resonance masses by sampling 
00308   //    P(M) = gc[g(M-M0)]**(c-1) *exp[-(g(M-M0))**c] 
00309   // for M > M0
00310   //
00311   // Use for the forward and backward clusters, but not
00312   // the cascade cluster
00313 
00314   const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
00315   const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
00316   G4int ntc = 0;
00317 
00318   if (forwardCount < 1 || backwardCount < 1) return false;  // array bounds protection
00319 
00320   G4double rmc = forwardMass;
00321   if (forwardCount > 1) {
00322     ntc = std::min(3,forwardCount-2);
00323     rmc += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
00324   }
00325   G4double rmd = backwardMass;
00326   if( backwardCount > 1 ) {
00327     ntc = std::min(3,backwardCount-2);
00328     rmd += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
00329   }
00330 
00331   while( rmc+rmd > centerofmassEnergy )
00332   {
00333     if( (rmc <= forwardMass) && (rmd <= backwardMass) )
00334     {
00335       G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
00336       rmc *= temp;
00337       rmd *= temp;
00338     }
00339     else
00340     {
00341       rmc = 0.1*forwardMass + 0.9*rmc;
00342       rmd = 0.1*backwardMass + 0.9*rmd;
00343     }
00344   }
00345 
00346   G4ReactionProduct pseudoParticle[8];
00347   for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
00348     
00349   pseudoParticle[1].SetMass( mOriginal*GeV );
00350   pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
00351   pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00352     
00353   pseudoParticle[2].SetMass( protonMass*MeV );
00354   pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
00355   pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
00356   //
00357   //  transform into center of mass system
00358   //
00359   pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
00360   pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
00361   pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
00362 
00363   // Calculate cm momentum for forward and backward masses
00364   // W = sqrt(pf*pf + rmc*rmc) + sqrt(pf*pf + rmd*rmd)
00365   // Solve for pf
00366 
00367   const G4double pfMin = 0.0001;
00368   G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
00369   pf *= pf;
00370   pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
00371   pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
00372   //
00373   //  set final state masses and energies in centre of mass system
00374   //
00375   pseudoParticle[3].SetMass( rmc*GeV );
00376   pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
00377     
00378   pseudoParticle[4].SetMass( rmd*GeV );
00379   pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
00380 
00381   //
00382   // Get cm scattering angle by sampling t from tmin to tmax
00383   //
00384   const G4double bMin = 0.01;
00385   const G4double b1 = 4.0;
00386   const G4double b2 = 1.6;
00387   G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
00388   G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
00389   G4double factor = 1.0 - std::exp(-dtb);
00390   G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
00391 
00392   costheta = std::max(-1.0, std::min(1.0, costheta));
00393   G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00394   G4double phi = G4UniformRand() * twopi;
00395   //
00396   // calculate final state momenta in centre of mass system
00397   //
00398   pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
00399                                  pf*sintheta*std::sin(phi)*GeV,
00400                                  pf*costheta*GeV );
00401   pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());
00402 
00403   // Backward cluster of nucleons and pions from intra-nuclear cascade
00404   // Set up in lab system and transform to cms
00405 
00406   G4double pp, pp1;
00407   if( nuclearExcitationCount > 0 )
00408   {
00409     const G4double ga = 1.2;
00410     G4double ekit1 = 0.04;
00411     G4double ekit2 = 0.6;   // Max KE of cascade particle
00412     if( ekOriginal <= 5.0 )
00413     {
00414       ekit1 *= ekOriginal*ekOriginal/25.0;
00415       ekit2 *= ekOriginal*ekOriginal/25.0;
00416     }
00417     G4double scale = std::pow(ekit2/ekit1, 1.0-ga) - 1.0;
00418     for( i=0; i<vecLen; ++i )
00419     {
00420       if( vec[i]->GetSide() == -2 )
00421       {
00422         G4double kineticE = ekit1*std::pow((1.0 + G4UniformRand()*scale), 1.0/(1.0-ga) );
00423         vec[i]->SetKineticEnergy( kineticE*GeV );
00424         G4double vMass = vec[i]->GetMass()/MeV;
00425         G4double totalE = kineticE*GeV + vMass;
00426         pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
00427         G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
00428         G4double sint = std::sqrt(1.0-cost*cost);
00429         phi = twopi*G4UniformRand();
00430         vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
00431                              pp*sint*std::sin(phi)*MeV,
00432                              pp*cost*MeV );
00433         vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
00434       }
00435     }
00436   }
00437 
00438   //
00439   // Fragmentation of forward and backward clusters
00440   //
00441 
00442   currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
00443   currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
00444     
00445   targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
00446   targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
00447     
00448   pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
00449   pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
00450   pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
00451     
00452   pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
00453   pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
00454   pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
00455   
00456   G4double wgt;
00457       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00458   if( forwardCount > 1 )     // tempV will contain the forward particles
00459   {
00460     G4FastVector<G4ReactionProduct,256> tempV;
00461     tempV.Initialize( forwardCount );
00462     G4bool constantCrossSection = true;
00463     G4int tempLen = 0;
00464     if( currentParticle.GetSide() == 1 )
00465       tempV.SetElement( tempLen++, &currentParticle );
00466     if( targetParticle.GetSide() == 1 )
00467       tempV.SetElement( tempLen++, &targetParticle );
00468     for( i=0; i<vecLen; ++i )
00469     {
00470       if( vec[i]->GetSide() == 1 )
00471       {
00472         if( tempLen < 18 )
00473           tempV.SetElement( tempLen++, vec[i] );
00474         else
00475         {
00476           vec[i]->SetSide( -1 );
00477           continue;
00478         }
00479       }
00480     }
00481     if( tempLen >= 2 )
00482     {
00483       wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
00484                                 constantCrossSection, tempV, tempLen );
00485       if( currentParticle.GetSide() == 1 )
00486         currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
00487       if( targetParticle.GetSide() == 1 )
00488         targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
00489       for( i=0; i<vecLen; ++i )
00490       {
00491         if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
00492       }
00493     }
00494   }
00495       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00496   if( backwardCount > 1 )   //  tempV will contain the backward particles,
00497   {                         //  but not those created from the intranuclear cascade
00498     G4FastVector<G4ReactionProduct,256> tempV;
00499     tempV.Initialize( backwardCount );
00500     G4bool constantCrossSection = true;
00501     G4int tempLen = 0;
00502     if( currentParticle.GetSide() == -1 )
00503       tempV.SetElement( tempLen++, &currentParticle );
00504     if( targetParticle.GetSide() == -1 )
00505       tempV.SetElement( tempLen++, &targetParticle );
00506     for( i=0; i<vecLen; ++i )
00507     {
00508       if( vec[i]->GetSide() == -1 )
00509       {
00510         if( tempLen < 18 )
00511           tempV.SetElement( tempLen++, vec[i] );
00512         else
00513         {
00514           vec[i]->SetSide( -2 );
00515           vec[i]->SetKineticEnergy( 0.0 );
00516           vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
00517           continue;
00518         }
00519       }
00520     }
00521     if( tempLen >= 2 )
00522     {
00523       wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
00524                                 constantCrossSection, tempV, tempLen );
00525       if( currentParticle.GetSide() == -1 )
00526         currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
00527       if( targetParticle.GetSide() == -1 )
00528         targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
00529       for( i=0; i<vecLen; ++i )
00530       {
00531         if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
00532       }
00533     }
00534   }
00535 
00536       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00537   //
00538   // Lorentz transformation in lab system
00539   //
00540   currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
00541   targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
00542   for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
00543 
00544       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00545   //
00546   // sometimes the leading strange particle is lost, set it back
00547   //
00548   G4bool dum = true;
00549   if( leadFlag )
00550   {
00551     // leadFlag will be true
00552     //  iff original particle is strange AND if incident particle is strange
00553     //  leadFlag is set to the incident particle
00554     //  or
00555     //  target particle is strange leadFlag is set to the target particle
00556 
00557     if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
00558       dum = false;
00559     else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
00560       dum = false;
00561     else
00562     {
00563       for( i=0; i<vecLen; ++i )
00564       {
00565         if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
00566         {
00567           dum = false;
00568           break;
00569         }
00570       }
00571     }
00572     if( dum )
00573     {
00574       G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
00575       G4double ekin;
00576       if( ((leadMass <  protonMass) && (targetParticle.GetMass()/MeV <  protonMass)) ||
00577           ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
00578       {
00579         ekin = targetParticle.GetKineticEnergy()/GeV;
00580         pp1 = targetParticle.GetMomentum().mag()/MeV; // old momentum
00581         targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
00582         targetParticle.SetKineticEnergy( ekin*GeV );
00583         pp = targetParticle.GetTotalMomentum()/MeV;   // new momentum
00584         if( pp1 < 1.0e-3 ) {
00585           G4ThreeVector iso = Isotropic(pp);
00586           targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00587         } else {
00588           targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
00589         }
00590         targetHasChanged = true;
00591       }
00592       else
00593       {
00594         ekin = currentParticle.GetKineticEnergy()/GeV;
00595         pp1 = currentParticle.GetMomentum().mag()/MeV;
00596         currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
00597         currentParticle.SetKineticEnergy( ekin*GeV );
00598         pp = currentParticle.GetTotalMomentum()/MeV;
00599         if( pp1 < 1.0e-3 ) {
00600           G4ThreeVector iso = Isotropic(pp);
00601           currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00602         } else {
00603           currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00604         }
00605         incidentHasChanged = true;
00606       }
00607     }
00608   }    // end of if( leadFlag )
00609 
00610   // Get number of final state nucleons and nucleons remaining in
00611   // target nucleus
00612     
00613   std::pair<G4int, G4int> finalStateNucleons = 
00614     GetFinalStateNucleons(originalTarget, vec, vecLen);
00615 
00616   G4int protonsInFinalState = finalStateNucleons.first;
00617   G4int neutronsInFinalState = finalStateNucleons.second;
00618 
00619   G4int numberofFinalStateNucleons = 
00620     protonsInFinalState + neutronsInFinalState;
00621 
00622   if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00623       targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00624       originalIncident->GetDefinition()->GetPDGMass() < 
00625                                  G4Lambda::Lambda()->GetPDGMass())
00626     numberofFinalStateNucleons++;
00627 
00628   numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
00629 
00630   G4int PinNucleus = std::max(0, 
00631     G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
00632   G4int NinNucleus = std::max(0,
00633     G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
00634   //
00635   //  for various reasons, the energy balance is not sufficient,
00636   //  check that,  energy balance, angle of final system, etc.
00637   //
00638   pseudoParticle[4].SetMass( mOriginal*GeV );
00639   pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
00640   pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00641     
00642   G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
00643   G4int diff = 0;
00644   if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() )  diff = 1;
00645   if(numberofFinalStateNucleons == 1) diff = 0;
00646   pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
00647   pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
00648   pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
00649     
00650   G4double theoreticalKinetic =
00651     pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
00652     
00653   pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
00654   pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
00655   pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
00656 
00657   if( vecLen < 16 )
00658   {
00659     G4ReactionProduct tempR[130];
00660     tempR[0] = currentParticle;
00661     tempR[1] = targetParticle;
00662     for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
00663 
00664     G4FastVector<G4ReactionProduct,256> tempV;
00665     tempV.Initialize( vecLen+2 );
00666     G4bool constantCrossSection = true;
00667     G4int tempLen = 0;
00668     for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
00669 
00670     if( tempLen >= 2 )
00671     {
00672       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00673       wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
00674                                 pseudoParticle[5].GetTotalEnergy()/MeV,
00675                                 constantCrossSection, tempV, tempLen );
00676       if (wgt == -1) {
00677         G4double Qvalue = 0;
00678         for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
00679         wgt = GenerateNBodyEvent( Qvalue/MeV,
00680                                   constantCrossSection, tempV, tempLen );
00681       }
00682       theoreticalKinetic = 0.0;
00683       for( i=0; i<vecLen+2; ++i )
00684       {
00685         pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
00686         pseudoParticle[7].SetMass( tempV[i]->GetMass() );
00687         pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
00688         pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
00689         theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
00690       }
00691     }
00692       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00693   }
00694   else
00695   {
00696     theoreticalKinetic -=
00697       ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
00698     for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
00699   }
00700   G4double simulatedKinetic =
00701     currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
00702   for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
00703 
00704   // make sure that kinetic energies are correct
00705   // the backward nucleon cluster is not produced within proper kinematics!!!
00706     
00707   if( simulatedKinetic != 0.0 )
00708   {
00709     wgt = (theoreticalKinetic)/simulatedKinetic;
00710     currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
00711     pp = currentParticle.GetTotalMomentum()/MeV;
00712     pp1 = currentParticle.GetMomentum().mag()/MeV;
00713     if( pp1 < 0.001*MeV ) {
00714       G4ThreeVector iso = Isotropic(pp);
00715       currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00716     } else {
00717       currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00718     }
00719 
00720     targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
00721     pp = targetParticle.GetTotalMomentum()/MeV;
00722     pp1 = targetParticle.GetMomentum().mag()/MeV;
00723     if( pp1 < 0.001*MeV ) {
00724       G4ThreeVector iso = Isotropic(pp);
00725       targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00726     } else {
00727       targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
00728     }
00729 
00730     for( i=0; i<vecLen; ++i )
00731     {
00732       vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
00733       pp = vec[i]->GetTotalMomentum()/MeV;
00734       pp1 = vec[i]->GetMomentum().mag()/MeV;
00735       if( pp1 < 0.001 ) {
00736         G4ThreeVector iso = Isotropic(pp);
00737         vec[i]->SetMomentum( iso.x(), iso.y(), iso.z() );
00738       } else {
00739         vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
00740       }
00741     }
00742   }
00743       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00744 
00745   Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
00746           modifiedOriginal, originalIncident, targetNucleus,
00747           currentParticle, targetParticle, vec, vecLen );
00748 
00749   //  Add black track particles
00750   //  the total number of particles produced is restricted to 198
00751   //  this may have influence on very high energies
00752 
00753   if( atomicWeight >= 1.5 )
00754   {
00755     // npnb is number of proton/neutron black track particles
00756     // ndta is the number of deuterons, tritons, and alphas produced
00757     // epnb is the kinetic energy available for proton/neutron black track 
00758     //   particles
00759     // edta is the kinetic energy available for deuteron/triton/alpha 
00760     //   particles
00761 
00762     G4int npnb = 0;
00763     G4int ndta = 0;
00764 
00765     G4double epnb, edta;
00766     if (veryForward) {
00767       epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
00768       edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
00769     } else {
00770       epnb = targetNucleus.GetPNBlackTrackEnergy();
00771       edta = targetNucleus.GetDTABlackTrackEnergy();
00772     }
00773 
00774     const G4double pnCutOff = 0.001;     // GeV
00775     const G4double dtaCutOff = 0.001;    // GeV
00776     //    const G4double kineticMinimum = 1.e-6;
00777     //    const G4double kineticFactor = -0.005;
00778       
00779     //    G4double sprob = 0.0;   // sprob = probability of self-absorption in 
00780                             // heavy molecules
00781     // Not currently used (DHW 9 June 2008)  const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
00782     //    if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
00783       
00784     if( epnb >= pnCutOff )
00785     {
00786       npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
00787       if( numberofFinalStateNucleons + npnb > atomicWeight )
00788         npnb = G4int(atomicWeight - numberofFinalStateNucleons);
00789       npnb = std::min( npnb, 127-vecLen );
00790     }
00791     if( edta >= dtaCutOff )
00792     {
00793       ndta = G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
00794       ndta = std::min( ndta, 127-vecLen );
00795     }
00796     if (npnb == 0 && ndta == 0) npnb = 1;
00797 
00798     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00799 
00800     AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal, 
00801                            PinNucleus, NinNucleus, targetNucleus,
00802                            vec, vecLen );
00803     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00804   }
00805 
00806   //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
00807   //  MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
00808   //
00809   //  calculate time delay for nuclear reactions
00810   //
00811   if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
00812     currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
00813   else
00814     currentParticle.SetTOF( 1.0 );
00815 
00816   return true;
00817 }

Public Member Functions
	G4RPGTwoCluster ()
G4bool	ReactionStage (const G4HadProjectile , G4ReactionProduct &, G4bool &, const G4DynamicParticle , G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
G4RPGTwoCluster Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation
