// // ******************************************************************** // * 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: G4DiffractiveExcitation.cc,v 1.22 2010/09/20 15:50:46 vuzhinsk Exp $ // ------------------------------------------------------------ // GEANT 4 class implemetation file // // ---------------- G4DiffractiveExcitation -------------- // by Gunter Folger, October 1998. // diffractive Excitation used by strings models // Take a projectile and a target // excite the projectile and target // Essential changed by V. Uzhinsky in November - December 2006 // in order to put it in a correspondence with original FRITIOF // model. Variant of FRITIOF with nucleon de-excitation is implemented. // Other changes by V.Uzhinsky in May 2007 were introduced to fit // meson-nucleon interactions. Additional changes by V. Uzhinsky // were introduced in December 2006. They treat diffraction dissociation // processes more exactly. // --------------------------------------------------------------------- #include "globals.hh" #include "Randomize.hh" #include "G4DiffractiveExcitation.hh" #include "G4FTFParameters.hh" #include "G4ElasticHNScattering.hh" #include "G4LorentzRotation.hh" #include "G4RotationMatrix.hh" #include "G4ThreeVector.hh" #include "G4ParticleDefinition.hh" #include "G4VSplitableHadron.hh" #include "G4ExcitedString.hh" #include "G4ParticleTable.hh" #include "G4Neutron.hh" #include "G4ParticleDefinition.hh" //#include "G4ios.hh" //#include "UZHI_diffraction.hh" G4DiffractiveExcitation::G4DiffractiveExcitation() { } // --------------------------------------------------------------------- G4bool G4DiffractiveExcitation:: ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const { // -------------------- Projectile parameters ----------------------- G4LorentzVector Pprojectile=projectile->Get4Momentum(); if(Pprojectile.z() < 0.) { target->SetStatus(2); return false; } G4double ProjectileRapidity = Pprojectile.rapidity(); G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding(); G4int absProjectilePDGcode=std::abs(ProjectilePDGcode); G4bool PutOnMassShell(false); // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation G4double M0projectile = Pprojectile.mag(); // Without de-excitation if(M0projectile < projectile->GetDefinition()->GetPDGMass()) { PutOnMassShell=true; M0projectile=projectile->GetDefinition()->GetPDGMass(); } G4double M0projectile2 = M0projectile * M0projectile; G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass(); G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass(); G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff(); // -------------------- Target parameters ------------------------- G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding(); G4int absTargetPDGcode=std::abs(TargetPDGcode); //G4cout<<"Excit "<Get4Momentum(); G4double M0target = Ptarget.mag(); // G4double TargetRapidity = Ptarget.rapidity(); if(M0target < target->GetDefinition()->GetPDGMass()) { PutOnMassShell=true; M0target=target->GetDefinition()->GetPDGMass(); } G4double M0target2 = M0target * M0target; G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass(); G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass(); G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff(); G4double AveragePt2=theParameters->GetAveragePt2(); // G4double ProbOfDiffraction=ProbProjectileDiffraction + // ProbTargetDiffraction; G4double SumMasses=M0projectile+M0target+200.*MeV; // Kinematical properties of the interactions -------------- G4LorentzVector Psum; // 4-momentum in CMS Psum=Pprojectile+Ptarget; G4double S=Psum.mag2(); // Transform momenta to cms and then rotate parallel to z axis; G4LorentzRotation toCms(-1*Psum.boostVector()); G4LorentzVector Ptmp=toCms*Pprojectile; if ( Ptmp.pz() <= 0. ) { target->SetStatus(2); // "String" moving backwards in CMS, abort collision !! return false; } toCms.rotateZ(-1*Ptmp.phi()); toCms.rotateY(-1*Ptmp.theta()); G4LorentzRotation toLab(toCms.inverse()); Pprojectile.transform(toCms); Ptarget.transform(toCms); G4double PZcms2, PZcms; G4double SqrtS=std::sqrt(S); if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses)) {target->SetStatus(2); return false;} // The model cannot work for // p+p-interactions // at Plab < 1.62 GeV/c. if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) {target->SetStatus(2); return false;} // The model cannot work for // Pi+p-interactions // at Plab < 1. GeV/c. if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) || (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) || (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) {target->SetStatus(2); return false;} // The model cannot work for // K+p-interactions // at Plab < ??? GeV/c. ??? PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) /4./S; if(PZcms2 < 0) {target->SetStatus(2); return false;} // It can be in an interaction with // off-shell nuclear nucleon PZcms = std::sqrt(PZcms2); if(PutOnMassShell) { if(Pprojectile.z() > 0.) { Pprojectile.setPz( PZcms); Ptarget.setPz( -PZcms); } else { Pprojectile.setPz(-PZcms); Ptarget.setPz( PZcms); }; Pprojectile.setE(std::sqrt(M0projectile2 + Pprojectile.x()*Pprojectile.x()+ Pprojectile.y()*Pprojectile.y()+ PZcms2)); Ptarget.setE(std::sqrt(M0target2 + Ptarget.x()*Ptarget.x()+ Ptarget.y()*Ptarget.y()+ PZcms2)); } G4double maxPtSquare; // = PZcms2; /* G4cout<<"Start --------------------"<GetMagQuarkExchange(); G4double SlopeQuarkExchange =theParameters->GetSlopeQuarkExchange(); G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange(); //G4cout<<"Q exc "<FindParticle(2112))->GetPDGMass(); G4double DeltaMass= (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass(); //G4cout<>Uzhi; // Check for possible quark exchange ----------------------------------- if(G4UniformRand() < MagQuarkExchange* std::exp(-SlopeQuarkExchange*ProjectileRapidity)) //TargetRapidity))) 1.45 { // std::exp(-SlopeQuarkExchange*(ProjectileRapidity - 1.36))) //TargetRapidity))) 1.45 //G4cout<<"Q exchange"< 0 ) // ProjQ1 is quark { ProjExchangeQ = ProjQ1; if(ProjExchangeQ != TargQ1) { TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ; } else if(ProjExchangeQ != TargQ2) { TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ; } else { TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ; } } else // ProjQ2 is quark { ProjExchangeQ = ProjQ2; if(ProjExchangeQ != TargQ1) { TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ; } else if(ProjExchangeQ != TargQ2) { TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ; } else { TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ; } } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark G4int aProjQ1=std::abs(ProjQ1); G4int aProjQ2=std::abs(ProjQ2); if(aProjQ1 == aProjQ2) {NewProjCode = 111;} // Pi0-meson else // |ProjQ1| # |ProjQ2| { if(aProjQ1 > aProjQ2) {NewProjCode = aProjQ1*100+aProjQ2*10+1;} else {NewProjCode = aProjQ2*100+aProjQ1*10+1;} } NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) && (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar else if( target->GetDefinition()->GetPDGiIsospin() == 3 ) //Delta was the target { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;} //Save Delta isobar else {} // De-excite initial Delta isobar } else if((G4UniformRand() < DeltaProbAtQuarkExchange) && //Nucleon was the target (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar else {} //Save initial nucleon // target->SetDefinition( // Fix 15.12.09 // G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09 } else { // projectile is baryon ---------------- G4double Same=0.; //0.3; //0.5; G4bool ProjDeltaHasCreated(false); G4bool TargDeltaHasCreated(false); G4double Ksi=G4UniformRand(); if(G4UniformRand() < 0.5) // Sampling exchange quark from proj. or targ. { // Sampling exchanged quark from the projectile --- if( Ksi < 0.333333 ) {ProjExchangeQ = ProjQ1;} else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) {ProjExchangeQ = ProjQ2;} else {ProjExchangeQ = ProjQ3;} //G4cout<<"ProjExchangeQ "<FindParticle(NewProjCode); if(TestParticle) TestParticleMass=TestParticle->GetPDGMass(); if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;} else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta { if(G4UniformRand() > DeltaProbAtQuarkExchange) {NewProjCode +=2; ProjDeltaHasCreated=true;} else {NewProjCode +=0; ProjDeltaHasCreated=false;} } else // Projectile was Nucleon { if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) {NewProjCode +=2; ProjDeltaHasCreated=true;} else {NewProjCode +=0; ProjDeltaHasCreated=false;} } G4ParticleDefinition* NewTestParticle= G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode); //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<GetPDGMass()<GetPDGMass()) {NewProjCode=TestParticleID;} //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // ***************************** //G4cout<<"TargQ1, TargQ2, TargQ3 "<FindParticle(NewTargCode); if(TestParticle) TestParticleMass=TestParticle->GetPDGMass(); if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta { if(G4UniformRand() > DeltaProbAtQuarkExchange) {NewTargCode +=2; TargDeltaHasCreated=true;} else {NewTargCode +=0; TargDeltaHasCreated=false;} } else // Target was Nucleon { if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2; TargDeltaHasCreated=true;} else {NewTargCode +=0; TargDeltaHasCreated=false;} } NewTestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode); //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<GetPDGMass()<GetPDGMass()) {NewTargCode=TestParticleID;} //G4cout<<"NewProjCode NewTargCode "<>Uzhi; if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode)) { // Nothing was changed! It is not right!? } // Forming baryons -------------------------------------------------- if(ProjDeltaHasCreated) {ProbProjectileDiffraction=1.; ProbTargetDiffraction=0.;} if(TargDeltaHasCreated) {ProbProjectileDiffraction=0.; ProbTargetDiffraction=1.;} if(ProjDeltaHasCreated) { M0projectile= (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass(); M0projectile2 = M0projectile * M0projectile; ProjectileDiffStateMinMass =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV } // if(M0target < // (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass()) if(TargDeltaHasCreated) { M0target= (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass(); M0target2 = M0target * M0target; TargetDiffStateMinMass =M0target+210.*MeV; //210 MeV=m_pi+70 MeV; TargetNonDiffStateMinMass=M0target+210.*MeV; //210 MeV=m_pi+70 MeV; } } // End of if projectile is baryon --------------------------- // If we assume that final state hadrons after the charge exchange will be // in the ground states, we have to put ---------------------------------- /* // if(M0projectile < // (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass()) if(ProjDeltaHasCreated) { M0projectile= (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass(); M0projectile2 = M0projectile * M0projectile; ProjectileDiffStateMinMass =M0projectile+160.*MeV; //160 MeV=m_pi+20 MeV ProjectileNonDiffStateMinMass=M0projectile+160.*MeV; //160 MeV=m_pi+20 MeV } // if(M0target < // (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass()) if(TargDeltaHasCreated) { M0target= (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass(); M0target2 = M0target * M0target; TargetDiffStateMinMass =M0target+160.*MeV; //160 MeV=m_pi+20 MeV; TargetNonDiffStateMinMass=M0target+160.*MeV; //160 MeV=m_pi+20 MeV; } */ PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) /4./S; //G4cout<<"PZcms2 1 "<SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); //---------------------------------------------------------- PZcms = std::sqrt(PZcms2); Pprojectile.setPz( PZcms); Pprojectile.setE(std::sqrt(M0projectile2+PZcms2)); Ptarget.setPz( -PZcms); Ptarget.setE(std::sqrt(M0target2+PZcms2)); // ---------------------------------------------------------- // G4double Wexcit=1.-1.97*std::exp(-0.5*ProjectileRapidity); G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity); //G4cout<>Uzhi; //Wexcit=0.; if(G4UniformRand() > Wexcit) { // Make elastic scattering //G4cout<<"Make elastic scattering"<SetTimeOfCreation(target->GetTimeOfCreation()); projectile->SetPosition(target->GetPosition()); projectile->Set4Momentum(Pprojectile); target->Set4Momentum(Ptarget); G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters); return Result; } // end of if(G4UniformRand() > Wexcit) } // End of charge exchange part ------------------------------ // ------------------------------------------------------------------ G4double ProbOfDiffraction=ProbProjectileDiffraction + ProbTargetDiffraction; /* G4cout<<"Excite --------------------"<>Uzhi; */ /* if(ProjectileNonDiffStateMinMass + TargetNonDiffStateMinMass > SqrtS) // 24.07.10 { if(ProbOfDiffraction!=0.) { ProbProjectileDiffraction/=ProbOfDiffraction; ProbOfDiffraction=1.; } else {return false;} } */ if(ProbOfDiffraction!=0.) { ProbProjectileDiffraction/=ProbOfDiffraction; } else { ProbProjectileDiffraction=0.; } //G4cout<<"Prob ProjDiff TargDiff "<= 500 && (whilecount%100)==0) // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" // << ", loop count/ maxPtSquare : " // << whilecount << " / " << maxPtSquare << G4endl; // whilecount++; if (whilecount > 1000 ) { Qmomentum=G4LorentzVector(0.,0.,0.,0.); target->SetStatus(2); return false; // Ignore this interaction }; // --------------- Check that the interaction is possible ----------- ProjMassT2=ProjectileDiffStateMinMass2; ProjMassT =ProjectileDiffStateMinMass; TargMassT2=M0target2; TargMassT =M0target; //G4cout<<"Masses "<SetStatus(2); return false; } maxPtSquare=PZcms2; Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); ProjMassT2=ProjectileDiffStateMinMass2+Pt2; ProjMassT =std::sqrt(ProjMassT2); TargMassT2=M0target2+Pt2; TargMassT =std::sqrt(TargMassT2); PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) /4./S; if(PZcms2 < 0 ) continue; PZcms =std::sqrt(PZcms2); PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; PMinusMax=SqrtS-TargMassT; PMinusNew=ChooseP(PMinusMin, PMinusMax); // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax)); //PMinusNew=1./sqr(1./std::sqrt(PMinusMin)-G4UniformRand()*(1./std::sqrt(PMinusMin)-1./std::sqrt(PMinusMax))); TMinusNew=SqrtS-PMinusNew; Qminus=Ptarget.minus()-TMinusNew; TPlusNew=TargMassT2/TMinusNew; Qplus=Ptarget.plus()-TPlusNew; Qmomentum.setPz( (Qplus-Qminus)/2 ); Qmomentum.setE( (Qplus+Qminus)/2 ); } while ((Pprojectile+Qmomentum).mag2() < ProjectileDiffStateMinMass2); //|| //Repeat the sampling because there was not any excitation //((Ptarget -Qmomentum).mag2() < M0target2 )) ); } else { // -------------- Target diffraction ---------------- //G4cout<<"Target diffraction"<= 500 && (whilecount%100)==0) // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" // << ", loop count/ maxPtSquare : " // << whilecount << " / " << maxPtSquare << G4endl; // whilecount++; if (whilecount > 1000 ) { Qmomentum=G4LorentzVector(0.,0.,0.,0.); target->SetStatus(2); return false; // Ignore this interaction }; //G4cout<<"Qm while "<SetStatus(2); return false; } maxPtSquare=PZcms2; Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); //G4cout<<"Qm while "<= 500 && (whilecount%100)==0) // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" // << ", loop count/ maxPtSquare : " // << whilecount << " / " << maxPtSquare << G4endl; // whilecount++; if (whilecount > 1000 ) { Qmomentum=G4LorentzVector(0.,0.,0.,0.); target->SetStatus(2); return false; // Ignore this interaction }; // --------------- Check that the interaction is possible ----------- ProjMassT2=ProjectileNonDiffStateMinMass2; ProjMassT =ProjectileNonDiffStateMinMass; TargMassT2=TargetNonDiffStateMinMass2; TargMassT =TargetNonDiffStateMinMass; PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) /4./S; if(PZcms2 < 0 ) { target->SetStatus(2); return false; } maxPtSquare=PZcms2; Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2; ProjMassT =std::sqrt(ProjMassT2); TargMassT2=TargetNonDiffStateMinMass2+Pt2; TargMassT =std::sqrt(TargMassT2); PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) /4./S; //G4cout<<"PZcms2 ND"<>Uzhi; */ } while ( ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction ((Ptarget -Qmomentum).mag2() < TargetNonDiffStateMinMass2 )); } Pprojectile += Qmomentum; Ptarget -= Qmomentum; //G4cout<<"Pr Y "<SetTimeOfCreation(target->GetTimeOfCreation()); projectile->SetPosition(target->GetPosition()); // Creation time and position of target nucleon were determined at // ReggeonCascade() of G4FTFModel // ------------------------------------------------------ //G4cout<<"Mproj "<Set4Momentum(Pprojectile); target->Set4Momentum(Ptarget); projectile->IncrementCollisionCount(1); target->IncrementCollisionCount(1); return true; } // --------------------------------------------------------------------- void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, G4bool isProjectile, G4ExcitedString * &FirstString, G4ExcitedString * &SecondString, G4FTFParameters *theParameters) const { hadron->SplitUp(); G4Parton *start= hadron->GetNextParton(); if ( start==NULL) { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; FirstString=0; SecondString=0; return; } G4Parton *end = hadron->GetNextParton(); if ( end==NULL) { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; FirstString=0; SecondString=0; return; } G4LorentzVector Phadron=hadron->Get4Momentum(); G4LorentzVector Pstart(0.,0.,0.,0.); G4LorentzVector Pend(0.,0.,0.,0.); G4LorentzVector Pkink(0.,0.,0.,0.); G4LorentzVector PkinkQ1(0.,0.,0.,0.); G4LorentzVector PkinkQ2(0.,0.,0.,0.); G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding()); G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding()); //-------------------------------------------------------------------------------- G4double Wmin(0.); if(isProjectile) { Wmin=theParameters->GetProjMinDiffMass(); } else { Wmin=theParameters->GetTarMinDiffMass(); } // end of if(isProjectile) G4double W = hadron->Get4Momentum().mag(); G4double W2=W*W; G4double Pt(0.), x1(0.), x2(0.), x3(0.); G4bool Kink=false; if(W > Wmin) { // Kink is possible G4double Pt2kink=theParameters->GetPt2Kink(); Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.)); if(Pt > 500.*MeV) { G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.)); G4double Y=Ymax*(1.- 2.*G4UniformRand()); x1=1.-Pt/W*std::exp( Y); x3=1.-Pt/W*std::exp(-Y); x2=2.-x1-x3; G4double Mass_startQ = 650.*MeV; if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV; if(PDGcode_startQ == 3) Mass_startQ = 500.*MeV; if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV; G4double Mass_endQ = 650.*MeV; if(PDGcode_endQ < 3) Mass_endQ = 325.*MeV; if(PDGcode_endQ == 3) Mass_endQ = 500.*MeV; if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV; G4double P2_1=W2*x1*x1/4.-Mass_endQ *Mass_endQ; G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ; G4double P2_2=sqr((2.-x1-x3)*W/2.); if((P2_1 <= 0.) || (P2_3 <= 0.)) { Kink=false;} else { G4double P_1=std::sqrt(P2_1); G4double P_2=std::sqrt(P2_2); G4double P_3=std::sqrt(P2_3); G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2); G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3); // Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) { Kink=false;} else { Kink=true; Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1); Pkink.setPx( Pt); Pkink.setPy( 0.); Pkink.setPz( P_2*CosT12); Pstart.setE(x3*W/2.); Pkink.setE(Pkink.vect().mag()); Pend.setE(x1*W/2.); G4double XkQ=GetQuarkFractionOfKink(0.,1.); if(Pkink.getZ() > 0.) { if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;} } else { if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;} } PkinkQ2=Pkink - PkinkQ1; //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------ G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/ std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13)); G4double Psi=std::acos(Cos2Psi); G4LorentzRotation Rotate; if(isProjectile) {Rotate.rotateY(Psi);} else {Rotate.rotateY(pi-Psi);} Rotate.rotateZ(twopi*G4UniformRand()); Pstart*=Rotate; Pkink*=Rotate; PkinkQ1*=Rotate; PkinkQ2*=Rotate; Pend*=Rotate; } } // end of if((P2_1 < 0.) || (P2_3 <0.)) } // end of if(Pt > 500.*MeV) } // end of if(W > Wmin) Check for a kink //-------------------------------------------------------------------------------- if(Kink) { // Kink is possible std::vector QuarkProbabilitiesAtGluonSplitUp = theParameters->GetQuarkProbabilitiesAtGluonSplitUp(); G4int QuarkInGluon(1); G4double Ksi=G4UniformRand(); for(unsigned int Iq=0; Iq <3; Iq++) { if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;} G4Parton * Gquark = new G4Parton(QuarkInGluon); G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon); //------------------------------------------------------------------------------- G4LorentzRotation toCMS(-1*Phadron.boostVector()); G4LorentzRotation toLab(toCMS.inverse()); Pstart.transform(toLab); start->Set4Momentum(Pstart); PkinkQ1.transform(toLab); PkinkQ2.transform(toLab); Pend.transform(toLab); end->Set4Momentum(Pend); G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding()); if(absPDGcode < 1000) { // meson if ( isProjectile ) { // Projectile if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end { // Quark on the end FirstString = new G4ExcitedString(end ,Ganti_quark, +1); SecondString= new G4ExcitedString(Gquark,start ,+1); Ganti_quark->Set4Momentum(PkinkQ1); Gquark->Set4Momentum(PkinkQ2); } else { // Anti_Quark on the end FirstString = new G4ExcitedString(end ,Gquark, +1); SecondString= new G4ExcitedString(Ganti_quark,start ,+1); Gquark->Set4Momentum(PkinkQ1); Ganti_quark->Set4Momentum(PkinkQ2); } // end of if(end->GetPDGcode() > 0) } else { // Target if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end { // Quark on the end FirstString = new G4ExcitedString(Ganti_quark,end ,-1); SecondString= new G4ExcitedString(start ,Gquark,-1); Ganti_quark->Set4Momentum(PkinkQ2); Gquark->Set4Momentum(PkinkQ1); } else { // Anti_Quark on the end FirstString = new G4ExcitedString(Gquark,end ,-1); SecondString= new G4ExcitedString(start ,Ganti_quark,-1); Gquark->Set4Momentum(PkinkQ2); Ganti_quark->Set4Momentum(PkinkQ1); } // end of if(end->GetPDGcode() > 0) } // end of if ( isProjectile ) } else // if(absPDGCode < 1000) { // Baryon/AntiBaryon if ( isProjectile ) { // Projectile if((end->GetDefinition()->GetParticleType() == "diquarks") && (end->GetDefinition()->GetPDGEncoding() > 0 ) ) { // DiQuark on the end FirstString = new G4ExcitedString(end ,Gquark, +1); SecondString= new G4ExcitedString(Ganti_quark,start ,+1); Gquark->Set4Momentum(PkinkQ1); Ganti_quark->Set4Momentum(PkinkQ2); } else { // Anti_DiQuark on the end or quark FirstString = new G4ExcitedString(end ,Ganti_quark, +1); SecondString= new G4ExcitedString(Gquark,start ,+1); Ganti_quark->Set4Momentum(PkinkQ1); Gquark->Set4Momentum(PkinkQ2); } // end of if(end->GetPDGcode() > 0) } else { // Target if((end->GetDefinition()->GetParticleType() == "diquarks") && (end->GetDefinition()->GetPDGEncoding() > 0 ) ) { // DiQuark on the end FirstString = new G4ExcitedString(Gquark,end ,-1); SecondString= new G4ExcitedString(start ,Ganti_quark,-1); Gquark->Set4Momentum(PkinkQ1); Ganti_quark->Set4Momentum(PkinkQ2); } else { // Anti_DiQuark on the end or Q FirstString = new G4ExcitedString(Ganti_quark,end ,-1); SecondString= new G4ExcitedString(start ,Gquark,-1); Gquark->Set4Momentum(PkinkQ2); Ganti_quark->Set4Momentum(PkinkQ1); } // end of if(end->GetPDGcode() > 0) } // end of if ( isProjectile ) } // end of if(absPDGcode < 1000) FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); FirstString->SetPosition(hadron->GetPosition()); SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation()); SecondString->SetPosition(hadron->GetPosition()); // ------------------------------------------------------------------------- } else // End of kink is possible { // Kink is impossible if ( isProjectile ) { FirstString= new G4ExcitedString(end,start, +1); } else { FirstString= new G4ExcitedString(start,end, -1); } SecondString=0; FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); FirstString->SetPosition(hadron->GetPosition()); // momenta of string ends // G4double Momentum=hadron->Get4Momentum().vect().mag(); G4double Plus=hadron->Get4Momentum().e() + Momentum; G4double Minus=hadron->Get4Momentum().e() - Momentum; G4ThreeVector tmp; if(Momentum > 0.) { tmp.set(hadron->Get4Momentum().px(), hadron->Get4Momentum().py(), hadron->Get4Momentum().pz()); tmp/=Momentum; } else { tmp.set(0.,0.,1.); } G4LorentzVector Pstart(tmp,0.); G4LorentzVector Pend(tmp,0.); if(isProjectile) { Pstart*=(-1.)*Minus/2.; Pend *=(+1.)*Plus /2.; } else { Pstart*=(+1.)*Plus/2.; Pend *=(-1.)*Minus/2.; } Momentum=-Pstart.mag(); Pstart.setT(Momentum); // It is assumed that quark has m=0. Momentum=-Pend.mag(); Pend.setT(Momentum); // It is assumed that di-quark has m=0. start->Set4Momentum(Pstart); end->Set4Momentum(Pend); SecondString=0; } // End of kink is impossible #ifdef G4_FTFDEBUG G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl; G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <Get4Momentum().mag() << G4endl; G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <Get4Momentum().mag()<< G4endl; G4cout << " sum of ends " << Pstart+Pend << G4endl; G4cout << " Original " << hadron->Get4Momentum() << G4endl; #endif return; } // --------- private methods ---------------------- // --------------------------------------------------------------------- G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const { // choose an x between Xmin and Xmax with P(x) ~ 1/x // to be improved... G4double range=Pmax-Pmin; if ( Pmin <= 0. || range <=0. ) { G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl; throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments "); } G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); return P; } // --------------------------------------------------------------------- G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const { // @@ this method is used in FTFModel as well. Should go somewhere common! G4double Pt2(0.); if(AveragePt2 <= 0.) {Pt2=0.;} else { Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.)); } G4double Pt=std::sqrt(Pt2); G4double phi=G4UniformRand() * twopi; return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); } // --------------------------------------------------------------------- G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const { G4double z, yf; do { z = zmin + G4UniformRand()*(zmax-zmin); yf = z*z +sqr(1 - z); } while (G4UniformRand() > yf); return z; } // --------------------------------------------------------------------- void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09 { G4int absIdPDG = std::abs(IdPDG); Q1= absIdPDG/ 100; Q2= (absIdPDG %100)/10; G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 ); if (IdPDG < 0 ) anti *=-1; Q1 *= anti; Q2 *= -1 * anti; return; } // --------------------------------------------------------------------- void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09 { Q1 = IdPDG / 1000; Q2 = (IdPDG % 1000) / 100; Q3 = (IdPDG % 100) / 10; return; } // --------------------------------------------------------------------- G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09 { G4int TmpQ(0); if( Q3 > Q2 ) { TmpQ = Q2; Q2 = Q3; Q3 = TmpQ; } else if( Q3 > Q1 ) { TmpQ = Q1; Q1 = Q3; Q3 = TmpQ; } if( Q2 > Q1 ) { TmpQ = Q1; Q1 = Q2; Q2 = TmpQ; } G4int NewCode = Q1*1000 + Q2* 100 + Q3* 10 + 2; return NewCode; } // --------------------------------------------------------------------- G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &) { throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called"); } G4DiffractiveExcitation::~G4DiffractiveExcitation() { } const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &) { throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called"); return *this; } int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const { throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called"); return false; } int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const { throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called"); return true; }