#include <G4DiffractiveExcitation.hh>
Public Member Functions | |
G4DiffractiveExcitation () | |
virtual | ~G4DiffractiveExcitation () |
virtual G4bool | ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const |
virtual void | CreateStrings (G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const |
Definition at line 48 of file G4DiffractiveExcitation.hh.
G4DiffractiveExcitation::G4DiffractiveExcitation | ( | ) |
G4DiffractiveExcitation::~G4DiffractiveExcitation | ( | ) | [virtual] |
void G4DiffractiveExcitation::CreateStrings | ( | G4VSplitableHadron * | aHadron, | |
G4bool | isProjectile, | |||
G4ExcitedString *& | FirstString, | |||
G4ExcitedString *& | SecondString, | |||
G4FTFParameters * | theParameters | |||
) | const [virtual] |
Definition at line 1076 of file G4DiffractiveExcitation.cc.
References G4cout, G4endl, G4UniformRand, G4Parton::Get4Momentum(), G4VSplitableHadron::Get4Momentum(), G4Parton::GetDefinition(), G4VSplitableHadron::GetNextParton(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetParticleType(), G4Parton::GetPDGcode(), G4ParticleDefinition::GetPDGEncoding(), G4VSplitableHadron::GetPosition(), G4FTFParameters::GetProjMinDiffMass(), G4FTFParameters::GetPt2Kink(), G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp(), G4VSplitableHadron::GetStatus(), G4FTFParameters::GetTarMinDiffMass(), G4VSplitableHadron::GetTimeOfCreation(), G4INCL::Math::pi, G4Parton::Set4Momentum(), G4ExcitedString::SetPosition(), G4ExcitedString::SetTimeOfCreation(), G4VSplitableHadron::SplitUp(), and sqr().
01081 { 01082 /* 01083 G4cout<<"Create Strings SplitUp "<<hadron<<G4endl; 01084 G4cout<<"Defin "<<hadron->GetDefinition()<<G4endl; 01085 G4cout<<"Defin "<<hadron->GetDefinition()->GetPDGEncoding()<<G4endl; 01086 */ 01087 hadron->SplitUp(); 01088 G4Parton *start= hadron->GetNextParton(); 01089 01090 if ( start==NULL) 01091 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; 01092 FirstString=0; SecondString=0; 01093 return; 01094 } 01095 G4Parton *end = hadron->GetNextParton(); 01096 if ( end==NULL) 01097 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; 01098 FirstString=0; SecondString=0; 01099 return; 01100 } 01101 //G4cout<<start<<" "<<start->GetPDGcode()<<" "<<end<<" "<<end->GetPDGcode()<<G4endl; 01102 //G4cout<<"Create string "<<start->GetPDGcode()<<" "<<end->GetPDGcode()<<G4endl; 01103 G4LorentzVector Phadron=hadron->Get4Momentum(); 01104 //G4cout<<"String mom "<<Phadron<<G4endl; 01105 G4LorentzVector Pstart(0.,0.,0.,0.); 01106 G4LorentzVector Pend(0.,0.,0.,0.); 01107 G4LorentzVector Pkink(0.,0.,0.,0.); 01108 G4LorentzVector PkinkQ1(0.,0.,0.,0.); 01109 G4LorentzVector PkinkQ2(0.,0.,0.,0.); 01110 01111 G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding()); 01112 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding()); 01113 //G4cout<<"PDGcode_startQ "<<PDGcode_startQ<<" PDGcode_endQ "<<PDGcode_endQ <<G4endl; 01114 01115 //-------------------------------------------------------------------------------- 01116 G4double Wmin(0.); 01117 if(isProjectile) 01118 { 01119 Wmin=theParameters->GetProjMinDiffMass(); 01120 } else 01121 { 01122 Wmin=theParameters->GetTarMinDiffMass(); 01123 } // end of if(isProjectile) 01124 01125 G4double W = hadron->Get4Momentum().mag(); 01126 //G4cout<<"Wmin W "<<Wmin<<" "<<W<<G4endl; 01127 G4double W2=W*W; 01128 01129 G4double Pt(0.), x1(0.), x3(0.); // x2(0.), 01130 G4bool Kink=false; 01131 01132 if(!(((start->GetDefinition()->GetParticleSubType() == "di_quark") && 01133 ( end->GetDefinition()->GetParticleSubType() == "di_quark") ) || 01134 ((start->GetDefinition()->GetParticleSubType() == "quark") && 01135 ( end->GetDefinition()->GetParticleSubType() == "quark") ))) 01136 { // Kinky strings are allowed only for qq-q strings 01137 // Kinky strings are impossible for other systems (qq-qqbar, q-qbar) 01138 // according to the analysis of Pbar P interactions 01139 //G4cout<<G4endl<<"Check for Kink!##############"<<G4endl<<G4endl; 01140 if(W > Wmin) 01141 { // Kink is possible 01142 if(hadron->GetStatus() == 0) // VU 10.04.2012 01143 { 01144 G4double Pt2kink=theParameters->GetPt2Kink(); // For non-diffractive 01145 Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.)); 01146 } 01147 else 01148 {Pt=0.;} // For diffractive 01149 01150 if(Pt > 500.*MeV) 01151 { 01152 G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.)); 01153 G4double Y=Ymax*(1.- 2.*G4UniformRand()); 01154 01155 x1=1.-Pt/W*std::exp( Y); 01156 x3=1.-Pt/W*std::exp(-Y); 01157 // x2=2.-x1-x3; 01158 01159 G4double Mass_startQ = 650.*MeV; 01160 if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV; 01161 if(PDGcode_startQ == 3) Mass_startQ = 500.*MeV; 01162 if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV; 01163 01164 G4double Mass_endQ = 650.*MeV; 01165 if(PDGcode_endQ < 3) Mass_endQ = 325.*MeV; 01166 if(PDGcode_endQ == 3) Mass_endQ = 500.*MeV; 01167 if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV; 01168 01169 G4double P2_1=W2*x1*x1/4.-Mass_endQ *Mass_endQ; 01170 G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ; 01171 01172 G4double P2_2=sqr((2.-x1-x3)*W/2.); 01173 01174 if((P2_1 <= 0.) || (P2_3 <= 0.)) 01175 { Kink=false;} 01176 else 01177 { 01178 G4double P_1=std::sqrt(P2_1); 01179 G4double P_2=std::sqrt(P2_2); 01180 G4double P_3=std::sqrt(P2_3); 01181 01182 G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2); 01183 G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3); 01184 // Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 01185 01186 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) 01187 { Kink=false;} 01188 else 01189 { 01190 Kink=true; 01191 Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 01192 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); 01193 Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1); 01194 Pkink.setPx( Pt); Pkink.setPy( 0.); Pkink.setPz( P_2*CosT12); 01195 Pstart.setE(x3*W/2.); 01196 Pkink.setE(Pkink.vect().mag()); 01197 Pend.setE(x1*W/2.); 01198 01199 G4double XkQ=GetQuarkFractionOfKink(0.,1.); 01200 if(Pkink.getZ() > 0.) 01201 { 01202 if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;} 01203 } else { 01204 if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;} 01205 } 01206 01207 PkinkQ2=Pkink - PkinkQ1; 01208 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------ 01209 01210 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/ 01211 std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13)); 01212 G4double Psi=std::acos(Cos2Psi); 01213 01214 G4LorentzRotation Rotate; 01215 if(isProjectile) {Rotate.rotateY(Psi);} 01216 else {Rotate.rotateY(pi-Psi);} 01217 Rotate.rotateZ(twopi*G4UniformRand()); 01218 01219 Pstart*=Rotate; 01220 Pkink*=Rotate; 01221 PkinkQ1*=Rotate; 01222 PkinkQ2*=Rotate; 01223 Pend*=Rotate; 01224 01225 } 01226 } // end of if((P2_1 < 0.) || (P2_3 <0.)) 01227 } // end of if(Pt > 500.*MeV) 01228 } // end of if(W > Wmin) Check for a kink 01229 } // end of qq-q string selection 01230 //-------------------------------------------------------------------------------- 01231 /* 01232 G4cout<<"Kink "<<Kink<<" " 01233 <<start->GetDefinition()->GetParticleSubType()<<" " 01234 << end->GetDefinition()->GetParticleSubType() <<G4endl; 01235 01236 G4cout<<"Kink "<<Kink<<" " 01237 <<start->GetDefinition()->GetPDGEncoding()<<" " 01238 << end->GetDefinition()->GetPDGEncoding() <<G4endl; 01239 G4int Uzhi; G4cin>>Uzhi; 01240 */ 01241 01242 if(Kink) 01243 { // Kink is possible 01244 //G4cout<<"Kink is sampled!"<<G4endl; 01245 std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp = 01246 theParameters->GetQuarkProbabilitiesAtGluonSplitUp(); 01247 01248 G4int QuarkInGluon(1); G4double Ksi=G4UniformRand(); 01249 for(unsigned int Iq=0; Iq <3; Iq++) 01250 { 01251 //G4cout<<"Iq "<<Iq<<G4endl; 01252 01253 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;} 01254 01255 //G4cout<<"Last Iq "<<QuarkInGluon<<G4endl; 01256 G4Parton * Gquark = new G4Parton(QuarkInGluon); 01257 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon); 01258 //G4cout<<"Lorentz "<<G4endl; 01259 01260 //------------------------------------------------------------------------------- 01261 G4LorentzRotation toCMS(-1*Phadron.boostVector()); 01262 01263 G4LorentzRotation toLab(toCMS.inverse()); 01264 //G4cout<<"Pstart "<<Pstart<<G4endl; 01265 //G4cout<<"Pend "<<Pend<<G4endl; 01266 Pstart.transform(toLab); start->Set4Momentum(Pstart); 01267 PkinkQ1.transform(toLab); 01268 PkinkQ2.transform(toLab); 01269 Pend.transform(toLab); end->Set4Momentum(Pend); 01270 //G4cout<<"Pstart "<<Pstart<<G4endl; 01271 //G4cout<<"Pend "<<Pend<<G4endl; 01272 //G4cout<<"Defin "<<hadron->GetDefinition()<<G4endl; 01273 //G4cout<<"Defin "<<hadron->GetDefinition()->GetPDGEncoding()<<G4endl; 01274 01275 // G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding()); 01276 G4int absPDGcode=1500; // 23 Dec 01277 if((start->GetDefinition()->GetParticleSubType() == "quark") && 01278 ( end->GetDefinition()->GetParticleSubType() == "quark") ) 01279 absPDGcode=110; 01280 01281 //G4cout<<"absPDGcode "<<absPDGcode<<G4endl; 01282 01283 if(absPDGcode < 1000) 01284 { // meson 01285 if ( isProjectile ) 01286 { // Projectile 01287 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end 01288 { // Quark on the end 01289 FirstString = new G4ExcitedString(end ,Ganti_quark, +1); 01290 SecondString= new G4ExcitedString(Gquark,start ,+1); 01291 Ganti_quark->Set4Momentum(PkinkQ1); 01292 Gquark->Set4Momentum(PkinkQ2); 01293 01294 } else 01295 { // Anti_Quark on the end 01296 FirstString = new G4ExcitedString(end ,Gquark, +1); 01297 SecondString= new G4ExcitedString(Ganti_quark,start ,+1); 01298 Gquark->Set4Momentum(PkinkQ1); 01299 Ganti_quark->Set4Momentum(PkinkQ2); 01300 01301 } // end of if(end->GetPDGcode() > 0) 01302 } else { // Target 01303 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end 01304 { // Quark on the end 01305 FirstString = new G4ExcitedString(Ganti_quark,end ,-1); 01306 SecondString= new G4ExcitedString(start ,Gquark,-1); 01307 Ganti_quark->Set4Momentum(PkinkQ2); 01308 Gquark->Set4Momentum(PkinkQ1); 01309 01310 } else 01311 { // Anti_Quark on the end 01312 FirstString = new G4ExcitedString(Gquark,end ,-1); 01313 SecondString= new G4ExcitedString(start ,Ganti_quark,-1); 01314 Gquark->Set4Momentum(PkinkQ2); 01315 Ganti_quark->Set4Momentum(PkinkQ1); 01316 01317 } // end of if(end->GetPDGcode() > 0) 01318 } // end of if ( isProjectile ) 01319 } else // if(absPDGCode < 1000) 01320 { // Baryon/AntiBaryon 01321 if ( isProjectile ) 01322 { // Projectile 01323 if((end->GetDefinition()->GetParticleType() == "diquarks") && 01324 (end->GetDefinition()->GetPDGEncoding() > 0 ) ) 01325 { // DiQuark on the end 01326 FirstString = new G4ExcitedString(end ,Gquark, +1); 01327 SecondString= new G4ExcitedString(Ganti_quark,start ,+1); 01328 Gquark->Set4Momentum(PkinkQ1); 01329 Ganti_quark->Set4Momentum(PkinkQ2); 01330 01331 } else 01332 { // Anti_DiQuark on the end or quark 01333 FirstString = new G4ExcitedString(end ,Ganti_quark, +1); 01334 SecondString= new G4ExcitedString(Gquark,start ,+1); 01335 Ganti_quark->Set4Momentum(PkinkQ1); 01336 Gquark->Set4Momentum(PkinkQ2); 01337 01338 } // end of if(end->GetPDGcode() > 0) 01339 } else { // Target 01340 01341 if((end->GetDefinition()->GetParticleType() == "diquarks") && 01342 (end->GetDefinition()->GetPDGEncoding() > 0 ) ) 01343 { // DiQuark on the end 01344 FirstString = new G4ExcitedString(Gquark,end ,-1); 01345 01346 SecondString= new G4ExcitedString(start ,Ganti_quark,-1); 01347 Gquark->Set4Momentum(PkinkQ1); 01348 Ganti_quark->Set4Momentum(PkinkQ2); 01349 01350 } else 01351 { // Anti_DiQuark on the end or Q 01352 FirstString = new G4ExcitedString(Ganti_quark,end ,-1); 01353 SecondString= new G4ExcitedString(start ,Gquark,-1); 01354 Gquark->Set4Momentum(PkinkQ2); 01355 Ganti_quark->Set4Momentum(PkinkQ1); 01356 01357 } // end of if(end->GetPDGcode() > 0) 01358 } // end of if ( isProjectile ) 01359 } // end of if(absPDGcode < 1000) 01360 01361 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); 01362 FirstString->SetPosition(hadron->GetPosition()); 01363 01364 SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation()); 01365 SecondString->SetPosition(hadron->GetPosition()); 01366 01367 // ------------------------------------------------------------------------- 01368 } else // End of kink is possible 01369 { // Kink is impossible 01370 //G4cout<<start<<" "<<start->GetPDGcode()<<" "<<end<<" "<<end->GetPDGcode()<<G4endl; 01371 if ( isProjectile ) 01372 { 01373 FirstString= new G4ExcitedString(end,start, +1); 01374 } else { 01375 FirstString= new G4ExcitedString(start,end, -1); 01376 } 01377 01378 SecondString=0; 01379 01380 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); 01381 FirstString->SetPosition(hadron->GetPosition()); 01382 01383 // momenta of string ends 01384 // 01385 G4double Momentum=hadron->Get4Momentum().vect().mag(); 01386 G4double Plus=hadron->Get4Momentum().e() + Momentum; 01387 G4double Minus=hadron->Get4Momentum().e() - Momentum; 01388 01389 G4ThreeVector tmp; 01390 if(Momentum > 0.) 01391 { 01392 tmp.set(hadron->Get4Momentum().px(), 01393 hadron->Get4Momentum().py(), 01394 hadron->Get4Momentum().pz()); 01395 tmp/=Momentum; 01396 } 01397 else 01398 { 01399 tmp.set(0.,0.,1.); 01400 } 01401 01402 G4LorentzVector Pstart1(tmp,0.); 01403 G4LorentzVector Pend1(tmp,0.); 01404 01405 if(isProjectile) 01406 { 01407 Pstart1*=(-1.)*Minus/2.; 01408 Pend1 *=(+1.)*Plus /2.; 01409 } 01410 else 01411 { 01412 Pstart1*=(+1.)*Plus/2.; 01413 Pend1 *=(-1.)*Minus/2.; 01414 } 01415 01416 Momentum=-Pstart1.mag(); 01417 Pstart1.setT(Momentum); // It is assumed that quark has m=0. 01418 01419 Momentum=-Pend1.mag(); 01420 Pend1.setT(Momentum); // It is assumed that di-quark has m=0. 01421 01422 start->Set4Momentum(Pstart1); 01423 end->Set4Momentum(Pend1); 01424 SecondString=0; 01425 } // End of kink is impossible 01426 01427 //G4cout<<"Quarks in the string at creation"<<FirstString->GetRightParton()->GetPDGcode()<<" "<<FirstString->GetLeftParton()->GetPDGcode()<<G4endl; 01428 //G4cout<<FirstString<<" "<<SecondString<<G4endl; 01429 01430 #ifdef G4_FTFDEBUG 01431 G4cout << " generated string flavors " 01432 << start->GetPDGcode() << " / " 01433 << end->GetPDGcode() << G4endl; 01434 G4cout << " generated string momenta: quark " 01435 << start->Get4Momentum() << "mass : " 01436 <<start->Get4Momentum().mag() << G4endl; 01437 G4cout << " generated string momenta: Diquark " 01438 << end ->Get4Momentum() 01439 << "mass : " <<end->Get4Momentum().mag()<< G4endl; 01440 G4cout << " sum of ends " << Pstart+Pend << G4endl; 01441 G4cout << " Original " << hadron->Get4Momentum() << G4endl; 01442 #endif 01443 01444 return; 01445 01446 }
G4bool G4DiffractiveExcitation::ExciteParticipants | ( | G4VSplitableHadron * | aPartner, | |
G4VSplitableHadron * | bPartner, | |||
G4FTFParameters * | theParameters, | |||
G4ElasticHNScattering * | theElastic | |||
) | const [virtual] |
Definition at line 74 of file G4DiffractiveExcitation.cc.
References G4ElasticHNScattering::ElasticScattering(), G4ParticleTable::FindParticle(), G4UniformRand, G4VSplitableHadron::Get4Momentum(), G4FTFParameters::GetAveragePt2(), G4VSplitableHadron::GetDefinition(), G4FTFParameters::GetDeltaProbAtQuarkExchange(), G4FTFParameters::GetMagQuarkExchange(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGiIsospin(), G4ParticleDefinition::GetPDGMass(), G4VSplitableHadron::GetPosition(), G4FTFParameters::GetProbabilityOfProjDiff(), G4FTFParameters::GetProbabilityOfTarDiff(), G4FTFParameters::GetProbLogDistr(), G4FTFParameters::GetProbOfSameQuarkExchange(), G4FTFParameters::GetProjMinDiffMass(), G4FTFParameters::GetProjMinNonDiffMass(), G4FTFParameters::GetSlopeQuarkExchange(), G4VSplitableHadron::GetStatus(), G4FTFParameters::GetTarMinDiffMass(), G4FTFParameters::GetTarMinNonDiffMass(), G4VSplitableHadron::GetTimeOfCreation(), G4VSplitableHadron::IncrementCollisionCount(), G4VSplitableHadron::Set4Momentum(), G4VSplitableHadron::SetDefinition(), G4VSplitableHadron::SetPosition(), G4VSplitableHadron::SetStatus(), and G4VSplitableHadron::SetTimeOfCreation().
00078 { 00079 //G4cout<<G4endl<<"ExciteParticipants --------------"<<G4endl; 00080 // -------------------- Projectile parameters ----------------------- 00081 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 00082 00083 if(Pprojectile.z() < 0.) 00084 { 00085 target->SetStatus(2); 00086 return false; 00087 } 00088 00089 G4double ProjectileRapidity = Pprojectile.rapidity(); 00090 00091 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding(); 00092 G4int absProjectilePDGcode=std::abs(ProjectilePDGcode); 00093 00094 G4bool PutOnMassShell(false); 00095 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation 00096 G4double M0projectile = Pprojectile.mag(); // Without de-excitation 00097 //G4cout<<"M0projectile "<<M0projectile<<" "<<ProjectileRapidity<<G4endl; 00098 00099 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 00100 { 00101 PutOnMassShell=true; 00102 M0projectile=projectile->GetDefinition()->GetPDGMass(); 00103 } 00104 00105 G4double M0projectile2 = M0projectile * M0projectile; 00106 00107 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass(); 00108 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass(); 00109 G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff(); 00110 00111 // -------------------- Target parameters ------------------------- 00112 G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding(); 00113 G4int absTargetPDGcode=std::abs(TargetPDGcode); 00114 //G4cout<<"Entry to QE or Excit "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl; 00115 00116 G4LorentzVector Ptarget=target->Get4Momentum(); 00117 //G4cout<<"Pproj "<<Pprojectile<<G4endl; 00118 //G4cout<<"Ptarget "<<Ptarget<<G4endl; 00119 G4double M0target = Ptarget.mag(); 00120 00121 // G4double TargetRapidity = Ptarget.rapidity(); 00122 //G4cout<<"M0target "<<M0target<<" "<<TargetRapidity<<G4endl; 00123 if(M0target < target->GetDefinition()->GetPDGMass()) 00124 { 00125 PutOnMassShell=true; 00126 M0target=target->GetDefinition()->GetPDGMass(); 00127 } 00128 00129 G4double M0target2 = M0target * M0target; 00130 00131 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass(); 00132 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass(); 00133 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff(); 00134 00135 G4double AveragePt2=theParameters->GetAveragePt2(); 00136 G4double ProbLogDistr=theParameters->GetProbLogDistr(); // Uzhi 21.05.2012 00137 00138 // G4double ProbOfDiffraction=ProbProjectileDiffraction + 00139 // ProbTargetDiffraction; 00140 00141 // G4double SumMasses=M0projectile+M0target+220.*MeV; // 200->220 7 June 2011 00142 G4double SumMasses=M0projectile+M0target+220.*MeV; // 200->220 7 June 2011 00143 00144 // Kinematical properties of the interactions -------------- 00145 G4LorentzVector Psum; // 4-momentum in CMS 00146 Psum=Pprojectile+Ptarget; 00147 G4double S=Psum.mag2(); 00148 00149 //Uzhi_SqrtS=std::sqrt(S); 00150 00151 // Transform momenta to cms and then rotate parallel to z axis; 00152 G4LorentzRotation toCms(-1*Psum.boostVector()); 00153 00154 G4LorentzVector Ptmp=toCms*Pprojectile; 00155 00156 if ( Ptmp.pz() <= 0. ) 00157 { 00158 target->SetStatus(2); 00159 // "String" moving backwards in CMS, abort collision !! 00160 return false; 00161 } 00162 00163 toCms.rotateZ(-1*Ptmp.phi()); 00164 toCms.rotateY(-1*Ptmp.theta()); 00165 00166 G4LorentzRotation toLab(toCms.inverse()); 00167 00168 Pprojectile.transform(toCms); 00169 Ptarget.transform(toCms); 00170 00171 G4double PZcms2, PZcms; 00172 00173 G4double SqrtS=std::sqrt(S); 00174 00175 /* 00176 G4cout<<"Proj "<<Pprojectile<<G4endl; 00177 G4cout<<"Targ "<<Ptarget<<G4endl; 00178 G4cout<<"SqrtS "<<SqrtS<<G4endl; 00179 G4cout<<"M0pr M0tr "<<M0projectile<<" "<<M0target<<" "<<SumMasses<<G4endl; 00180 */ 00181 //G4cout<<"SqrtS < 2300*MeV Bary "<<SqrtS<<G4endl; 00182 // if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses)) // 7.06.11 00183 if(absProjectilePDGcode > 1000 && (SqrtS < SumMasses)) 00184 {target->SetStatus(2); return false;} // The model cannot work for 00185 // p+p-interactions 00186 // at Plab < 1.62 GeV/c. 00187 00188 //G4cout<<"SqrtS < 1600*MeV Pion "<<SqrtS<<G4endl; 00189 00190 // if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) && // 7.06.11 00191 // ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 00192 if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) && 00193 (SqrtS < SumMasses)) 00194 {target->SetStatus(2); return false;} // The model cannot work for 00195 // Pi+p-interactions 00196 // at Plab < 1. GeV/c. 00197 //G4cout<<"SqrtS < 1600*MeV "<<SqrtS<<G4endl; 00198 //SumMasses=M0projectile+M0target+20.*MeV; 00199 if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) || 00200 (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) || 00201 (absProjectilePDGcode == 310)) && (SqrtS < SumMasses)) 00202 // (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 00203 {target->SetStatus(2); return false;} // The model cannot work for 00204 // K+p-interactions 00205 // at Plab < ??? GeV/c. ??? 00206 00207 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- 00208 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) 00209 /4./S; 00210 //G4cout<<"PZcms2 "<<PZcms2<<G4endl; 00211 if(PZcms2 < 0) 00212 {target->SetStatus(2); return false;} // It can be in an interaction with 00213 // off-shell nuclear nucleon 00214 00215 PZcms = std::sqrt(PZcms2); 00216 00217 if(PutOnMassShell) 00218 { 00219 if(Pprojectile.z() > 0.) 00220 { 00221 Pprojectile.setPz( PZcms); 00222 Ptarget.setPz( -PZcms); 00223 } 00224 else 00225 { 00226 Pprojectile.setPz(-PZcms); 00227 Ptarget.setPz( PZcms); 00228 }; 00229 00230 Pprojectile.setE(std::sqrt(M0projectile2 + 00231 Pprojectile.x()*Pprojectile.x()+ 00232 Pprojectile.y()*Pprojectile.y()+ 00233 PZcms2)); 00234 Ptarget.setE(std::sqrt(M0target2 + 00235 Ptarget.x()*Ptarget.x()+ 00236 Ptarget.y()*Ptarget.y()+ 00237 PZcms2)); 00238 } 00239 00240 00241 //G4cout<<"Proj "<<Pprojectile<<G4endl; 00242 //G4cout<<"Targ "<<Ptarget<<G4endl; 00243 G4double maxPtSquare; // = PZcms2; 00244 /* 00245 Uzhi_targetdiffraction=0; 00246 Uzhi_projectilediffraction=0; 00247 Uzhi_Mx2=1.0; 00248 Uzhi_modT=0.; 00249 G4int Uzhi_QE=0; 00250 */ 00251 /* 00252 G4cout<<"Start --------------------"<<G4endl; 00253 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<G4endl; 00254 G4cout<<"Targ "<<M0target <<" "<<TargetDiffStateMinMass <<" "<<TargetNonDiffStateMinMass<<G4endl; 00255 G4cout<<"SqrtS "<<SqrtS<<G4endl; 00256 G4cout<<"Rapid "<<ProjectileRapidity<<G4endl; //" "<<TargetRapidity<<G4endl; 00257 */ 00258 // Charge exchange can be possible for baryons ----------------- 00259 00260 // Getting the values needed for exchange ---------------------- 00261 G4double MagQuarkExchange =theParameters->GetMagQuarkExchange();//0.045; // 00262 G4double SlopeQuarkExchange =theParameters->GetSlopeQuarkExchange();//0.; // 00263 // 777777777777777777777777777777777777777777777777777777777777777777777777777777 00264 G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange(); 00265 00266 // G4double NucleonMass= 00267 // (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass(); 00268 G4double DeltaMass= 00269 (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass(); 00270 00271 //G4double TargetRapidity(0.); 00272 //G4cout<<"Prob Q Exch "<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*ProjectileRapidity)<<G4endl; 00273 00274 //G4cout<<"Q exc Mag Slop Wdelta"<<MagQuarkExchange<<" "<<SlopeQuarkExchange<<" "<<DeltaProbAtQuarkExchange<<G4endl; 00275 //G4cout<<"ProjectileRapidity "<<ProjectileRapidity<<G4endl; 00276 //G4cout<<"Prob Exc "<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity))<<G4endl; 00277 //G4int Uzhi; G4cin>>Uzhi; 00278 // Check for possible quark exchange ----------------------------------- 00279 00280 if(G4UniformRand() < MagQuarkExchange* 00281 std::exp(-SlopeQuarkExchange*ProjectileRapidity)) //TargetRapidity))) 1.45 00282 { 00283 // std::exp(-SlopeQuarkExchange*(ProjectileRapidity - 1.36))) //TargetRapidity))) 1.45 00284 //G4cout<<"Q exchange"<<G4endl; 00285 //Uzhi_QE=1; 00286 G4int NewProjCode(0), NewTargCode(0); 00287 00288 G4int ProjQ1(0), ProjQ2(0), ProjQ3(0); 00289 00290 // Projectile unpacking -------------------------- 00291 if(absProjectilePDGcode < 1000 ) 00292 { // projectile is meson ----------------- 00293 UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2); 00294 } else 00295 { // projectile is baryon ---------------- 00296 UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3); 00297 } // End of the hadron's unpacking ---------- 00298 00299 // Target unpacking ------------------------------ 00300 G4int TargQ1(0), TargQ2(0), TargQ3(0); 00301 UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3); 00302 00303 //G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl; 00304 //G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl; 00305 // Sampling of exchanged quarks ------------------- 00306 G4int ProjExchangeQ(0); 00307 G4int TargExchangeQ(0); 00308 00309 if(absProjectilePDGcode < 1000 ) 00310 { // projectile is meson ----------------- 00311 00312 if(ProjQ1 > 0 ) // ProjQ1 is quark 00313 { 00314 G4int Navailable=0; 00315 ProjExchangeQ = ProjQ1; 00316 if(ProjExchangeQ != TargQ1) Navailable++; 00317 if(ProjExchangeQ != TargQ2) Navailable++; 00318 if(ProjExchangeQ != TargQ3) Navailable++; 00319 00320 G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1; 00321 //G4cout<<"Navailable Nsampled "<<Navailable<<" "<<Nsampled<<G4endl; 00322 Navailable=0; 00323 if(ProjExchangeQ != TargQ1) 00324 { 00325 Navailable++; 00326 if(Navailable == Nsampled) 00327 {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;} 00328 } 00329 00330 if(ProjExchangeQ != TargQ2) 00331 { 00332 Navailable++; 00333 if(Navailable == Nsampled) 00334 {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;} 00335 } 00336 00337 if(ProjExchangeQ != TargQ3) 00338 { 00339 Navailable++; 00340 if(Navailable == Nsampled) 00341 {TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;} 00342 } 00343 } else // ProjQ2 is quark 00344 { 00345 G4int Navailable=0; 00346 ProjExchangeQ = ProjQ2; 00347 if(ProjExchangeQ != TargQ1) Navailable++; 00348 if(ProjExchangeQ != TargQ2) Navailable++; 00349 if(ProjExchangeQ != TargQ3) Navailable++; 00350 00351 G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1; 00352 //G4cout<<"Navailable Nsampled "<<Navailable<<" "<<Nsampled<<G4endl; 00353 Navailable=0; 00354 if(ProjExchangeQ != TargQ1) 00355 { 00356 Navailable++; 00357 if(Navailable == Nsampled) 00358 {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;} 00359 } 00360 00361 if(ProjExchangeQ != TargQ2) 00362 { 00363 Navailable++; 00364 if(Navailable == Nsampled) 00365 {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;} 00366 } 00367 00368 if(ProjExchangeQ != TargQ3) 00369 { 00370 Navailable++; 00371 if(Navailable == Nsampled) 00372 {TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;} 00373 } 00374 } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark 00375 00376 //G4cout<<"Exch Pr Tr "<<ProjExchangeQ<<" "<<TargExchangeQ<<G4endl; 00377 //G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl; 00378 //G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl; 00379 00380 G4int aProjQ1=std::abs(ProjQ1); 00381 G4int aProjQ2=std::abs(ProjQ2); 00382 if(aProjQ1 == aProjQ2) {NewProjCode = 111;} // Pi0-meson 00383 else // |ProjQ1| # |ProjQ2| 00384 { 00385 if(aProjQ1 > aProjQ2) {NewProjCode = aProjQ1*100+aProjQ2*10+1;} 00386 else {NewProjCode = aProjQ2*100+aProjQ1*10+1;} 00387 // NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode); 00388 } 00389 00390 G4bool ProjExcited=false; 00391 //G4cout<<"NewProjCode "<<NewProjCode<<G4endl; 00392 if(G4UniformRand() < 0.5) 00393 { 00394 NewProjCode +=2; // Excited Pi0-meson 00395 ProjExcited=true; 00396 } 00397 if(aProjQ1 != aProjQ2) NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode); // Uzhi 00398 //G4cout<<"NewProjCode +2 or 0 "<<NewProjCode<<G4endl; 00399 00400 G4ParticleDefinition* TestParticle=0; 00401 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode); 00402 //G4cout<<"TestParticle ? "<<TestParticle<<G4endl; 00403 00404 if(TestParticle) 00405 { 00406 G4double MtestPart= // 31.05.2012 00407 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass(); 00408 /* 00409 G4cout<<"TestParticle Name "<<NewProjCode<<" "<<TestParticle->GetParticleName()<<G4endl; 00410 G4cout<<"MtestPart M0projectile projectile->GetDefinition()->GetPDGMass() "<<MtestPart<<" "<<M0projectile<<" "<<projectile->GetDefinition()->GetPDGMass()<<G4endl; 00411 G4bool Test =M0projectile <= projectile->GetDefinition()->GetPDGMass(); 00412 G4cout<<"M0projectile <= projectile->GetDefinition()->GetPDGMass() "<<Test<<G4endl; 00413 */ 00414 00415 if(MtestPart > M0projectile) 00416 {M0projectile = MtestPart;} 00417 else 00418 { 00419 if(std::abs(M0projectile - projectile->GetDefinition()->GetPDGMass()) < 140.*MeV) 00420 {M0projectile = MtestPart;} 00421 } 00422 //G4cout<<"M0projectile After check "<<M0projectile<<G4endl; 00423 M0projectile2 = M0projectile * M0projectile; 00424 00425 ProjectileDiffStateMinMass =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV 00426 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV 00427 } else 00428 {return false;} 00429 00430 //G4cout<<"New TrQ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl; 00431 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); 00432 //G4cout<<"NewTargCode "<<NewTargCode<<G4endl; 00433 00434 // if( (TargQ1 != TargQ2) && (TargQ1 != TargQ3) && (TargQ2 != TargQ3) // Lambda or Sigma0 00435 // {if(G4UniformRand() < 0.5) NewTargCode= 00436 00437 00438 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) && 00439 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2; //Create Delta isobar 00440 ProjExcited=true;} 00441 else if( target->GetDefinition()->GetPDGiIsospin() == 3 ) //Delta was the target 00442 { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2; //Save Delta isobar 00443 ProjExcited=true;} 00444 else {} // De-excite initial Delta isobar 00445 } 00446 00447 // else if((!CreateDelta) && 00448 else if((!ProjExcited) && 00449 (G4UniformRand() < DeltaProbAtQuarkExchange) && //Nucleon was the target 00450 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar 00451 // else if( CreateDelta) {NewTargCode +=2;} 00452 else {} //Save initial nucleon 00453 00454 // target->SetDefinition( // Fix 15.12.09 00455 // G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09 00456 00457 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode); 00458 //G4cout<<"New targ "<<NewTargCode<<" "<<TestParticle->GetParticleName()<<G4endl; 00459 if(TestParticle) 00460 { 00461 G4double MtestPart= // 31.05.2012 00462 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass(); 00463 00464 if(MtestPart > M0target) 00465 {M0target=MtestPart;} 00466 else 00467 { 00468 if(std::abs(M0target - target->GetDefinition()->GetPDGMass()) < 140.*MeV) 00469 {M0target=MtestPart;} 00470 } 00471 00472 TargetDiffStateMinMass =M0target+220.*MeV; //220 MeV=m_pi+80 MeV; 00473 TargetNonDiffStateMinMass=M0target+220.*MeV; //220 MeV=m_pi+80 MeV; 00474 } else 00475 {return false;} 00476 } else 00477 { // projectile is baryon ---------------- 00478 //========================================================================= 00479 G4double Same=theParameters->GetProbOfSameQuarkExchange(); //0.3; //0.5; 0. 00480 G4bool ProjDeltaHasCreated(false); 00481 G4bool TargDeltaHasCreated(false); 00482 00483 G4double Ksi=G4UniformRand(); 00484 if(G4UniformRand() < 0.5) // Sampling exchange quark from proj. or targ. 00485 { // Sampling exchanged quark from the projectile --- 00486 if( Ksi < 0.333333 ) 00487 {ProjExchangeQ = ProjQ1;} 00488 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 00489 {ProjExchangeQ = ProjQ2;} 00490 else 00491 {ProjExchangeQ = ProjQ3;} 00492 00493 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 00494 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) 00495 { 00496 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 00497 } else 00498 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) 00499 { 00500 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 00501 } else 00502 { 00503 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 00504 } 00505 00506 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 00507 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl; 00508 if( Ksi < 0.333333 ) 00509 {ProjQ1=ProjExchangeQ;} 00510 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 00511 {ProjQ2=ProjExchangeQ;} 00512 else 00513 {ProjQ3=ProjExchangeQ;} 00514 00515 } else 00516 { // Sampling exchanged quark from the target ------- 00517 if( Ksi < 0.333333 ) 00518 {TargExchangeQ = TargQ1;} 00519 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 00520 {TargExchangeQ = TargQ2;} 00521 else 00522 {TargExchangeQ = TargQ3;} 00523 00524 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) 00525 { 00526 ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 00527 } else 00528 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) 00529 { 00530 ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 00531 } else 00532 { 00533 ProjExchangeQ = ProjQ3; ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 00534 } 00535 00536 if( Ksi < 0.333333 ) 00537 {TargQ1=TargExchangeQ;} 00538 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 00539 {TargQ2=TargExchangeQ;} 00540 else 00541 {TargQ3=TargExchangeQ;} 00542 00543 } // End of sampling baryon 00544 00545 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // ***************************** 00546 00547 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;} 00548 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta 00549 { if(G4UniformRand() > DeltaProbAtQuarkExchange) 00550 {NewProjCode +=2; ProjDeltaHasCreated=true;} 00551 else {NewProjCode +=0; ProjDeltaHasCreated=false;} 00552 } 00553 else // Projectile was Nucleon 00554 { 00555 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) 00556 {NewProjCode +=2; ProjDeltaHasCreated=true;} 00557 else {NewProjCode +=0; ProjDeltaHasCreated=false;} 00558 } 00559 00560 //G4ParticleDefinition* NewTestParticle= 00561 // G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode); 00562 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl; 00563 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewProjCode=TestParticleID;} 00564 00565 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= 00566 00567 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // ***************************** 00568 00569 //G4cout<<"TargQ1, TargQ2, TargQ3 "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<" "<<NewTargCode<<G4endl; 00570 00571 //TestParticleID=NewTargCode; 00572 //TestParticleMass=DBL_MAX; 00573 00574 //TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode); 00575 //if(TestParticle) TestParticleMass=TestParticle->GetPDGMass(); 00576 00577 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 00578 else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta 00579 { if(G4UniformRand() > DeltaProbAtQuarkExchange) 00580 {NewTargCode +=2; TargDeltaHasCreated=true;} 00581 else {NewTargCode +=0; TargDeltaHasCreated=false;} 00582 } 00583 else // Target was Nucleon 00584 { 00585 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) 00586 {NewTargCode +=2; TargDeltaHasCreated=true;} 00587 else {NewTargCode +=0; TargDeltaHasCreated=false;} 00588 } 00589 00590 //NewTestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode); 00591 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl; 00592 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewTargCode=TestParticleID;} 00593 00594 //G4cout<<"NewProjCode NewTargCode "<<NewProjCode<<" "<<NewTargCode<<G4endl; 00595 //G4int Uzhi; G4cin>>Uzhi; 00596 00597 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode)) 00598 { // Nothing was changed! It is not right!? 00599 } 00600 // Forming baryons -------------------------------------------------- 00601 if(ProjDeltaHasCreated) {ProbProjectileDiffraction=1.; ProbTargetDiffraction=0.;} 00602 if(TargDeltaHasCreated) {ProbProjectileDiffraction=0.; ProbTargetDiffraction=1.;} 00603 if(ProjDeltaHasCreated) 00604 { 00605 G4double MtestPart= // 31.05.2012 00606 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass(); 00607 00608 if(MtestPart >= M0projectile) // 31.05.2012 00609 { // 31.05.2012 00610 M0projectile = MtestPart; // 31.05.2012 00611 M0projectile2 = M0projectile * M0projectile; // 31.05.2012 00612 } // 31.05.2012 00613 00614 ProjectileDiffStateMinMass =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV 00615 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV 00616 } 00617 00618 // if(M0target < 00619 // (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass()) 00620 if(TargDeltaHasCreated) 00621 { 00622 G4double MtestPart= // 31.05.2012 00623 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass(); 00624 00625 if(MtestPart >=M0target) // 31.05.2012 00626 { // 31.05.2012 00627 M0target=MtestPart; // 31.05.2012 00628 M0target2 = M0target * M0target; // 31.05.2012 00629 } // 31.05.2012 00630 00631 TargetDiffStateMinMass =M0target+210.*MeV; //210 MeV=m_pi+70 MeV; 00632 TargetNonDiffStateMinMass=M0target+210.*MeV; //210 MeV=m_pi+70 MeV; 00633 } 00634 } // End of if projectile is baryon --------------------------- 00635 00636 //G4cout<<"At end// NewProjCode "<<NewProjCode<<G4endl; 00637 //G4cout<<"At end// NewTargCode "<<NewTargCode<<G4endl; 00638 00639 // If we assume that final state hadrons after the charge exchange will be 00640 // in the ground states, we have to put ---------------------------------- 00641 //G4cout<<"M0pr M0tr SqS "<<M0projectile<<" "<<M0target<<" "<<SqrtS<<G4endl; 00642 00643 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- 00644 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) 00645 /4./S; 00646 //G4cout<<"PZcms2 1 "<<PZcms2<<G4endl<<G4endl; 00647 if(PZcms2 < 0) {return false;} // It can be if energy is not sufficient for Delta 00648 //---------------------------------------------------------- 00649 projectile->SetDefinition( 00650 G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); 00651 00652 target->SetDefinition( 00653 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 00654 //---------------------------------------------------------- 00655 PZcms = std::sqrt(PZcms2); 00656 00657 Pprojectile.setPz( PZcms); 00658 Pprojectile.setE(std::sqrt(M0projectile2+PZcms2)); 00659 00660 Ptarget.setPz( -PZcms); 00661 Ptarget.setE(std::sqrt(M0target2+PZcms2)); 00662 00663 // ---------------------------------------------------------- 00664 00665 if(absProjectilePDGcode < 1000) 00666 { // For projectile meson 00667 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity); 00668 Wexcit=0.; 00669 if(G4UniformRand() > Wexcit) 00670 { // Make elastic scattering 00671 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl; 00672 Pprojectile.transform(toLab); 00673 Ptarget.transform(toLab); 00674 00675 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); 00676 projectile->SetPosition(target->GetPosition()); 00677 00678 projectile->Set4Momentum(Pprojectile); 00679 target->Set4Momentum(Ptarget); 00680 00681 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters); 00682 //G4cout<<"Result of el. scatt "<<Result<<G4endl; 00683 return Result; 00684 } // end of if(Make elastic scattering for projectile meson?) 00685 } else 00686 { // For projectile baryon 00687 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity); 00688 //Wexcit=0.; 00689 if(G4UniformRand() > Wexcit) 00690 { // Make elastic scattering 00691 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl; 00692 Pprojectile.transform(toLab); 00693 Ptarget.transform(toLab); 00694 00695 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); 00696 projectile->SetPosition(target->GetPosition()); 00697 00698 projectile->Set4Momentum(Pprojectile); 00699 target->Set4Momentum(Ptarget); 00700 00701 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters); 00702 return Result; 00703 } // end of if(Make elastic scattering for projectile baryon?) 00704 } 00705 //G4cout<<"Make excitation of new hadrons"<<G4endl; 00706 } // End of charge exchange part ------------------------------ 00707 00708 // ------------------------------------------------------------------ 00709 G4double ProbOfDiffraction=ProbProjectileDiffraction + ProbTargetDiffraction; 00710 /* 00711 G4cout<<"Excite --------------------"<<G4endl; 00712 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<G4endl; 00713 G4cout<<"Targ "<<M0target <<" "<<TargetDiffStateMinMass <<" "<<TargetNonDiffStateMinMass<<G4endl; 00714 G4cout<<"SqrtS "<<SqrtS<<G4endl; 00715 00716 G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl; 00717 G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl; 00718 //G4int Uzhi; G4cin>>Uzhi; 00719 */ 00720 /* 00721 if(ProjectileNonDiffStateMinMass + TargetNonDiffStateMinMass > SqrtS) // 24.07.10 00722 { 00723 if(ProbOfDiffraction!=0.) 00724 { 00725 ProbProjectileDiffraction/=ProbOfDiffraction; 00726 ProbOfDiffraction=1.; 00727 } else {return false;} 00728 } 00729 00730 */ 00731 00732 if(ProbOfDiffraction!=0.) 00733 { 00734 ProbProjectileDiffraction/=ProbOfDiffraction; 00735 } 00736 else 00737 { 00738 ProbProjectileDiffraction=0.; 00739 } 00740 00741 //G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl; 00742 00743 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass * 00744 ProjectileDiffStateMinMass; 00745 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass * 00746 ProjectileNonDiffStateMinMass; 00747 00748 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass * 00749 TargetDiffStateMinMass; 00750 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass * 00751 TargetNonDiffStateMinMass; 00752 00753 G4double Pt2; 00754 G4double ProjMassT2, ProjMassT; 00755 G4double TargMassT2, TargMassT; 00756 G4double PMinusMin, PMinusMax; 00757 // G4double PPlusMin , PPlusMax; 00758 G4double TPlusMin , TPlusMax; 00759 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew; 00760 00761 G4LorentzVector Qmomentum; 00762 G4double Qminus, Qplus; 00763 00764 G4int whilecount=0; 00765 00766 // Choose a process --------------------------- 00767 //ProbOfDiffraction=1.; // Uzhi Difr 00768 //ProbProjectileDiffraction=1.; // Uzhi 00769 if(G4UniformRand() < ProbOfDiffraction) 00770 { 00771 if(G4UniformRand() < ProbProjectileDiffraction) 00772 { //-------- projectile diffraction --------------- 00773 //G4cout<<"projectile diffraction"<<G4endl; 00774 00775 do { 00776 /* 00777 Uzhi_projectilediffraction=1; 00778 Uzhi_targetdiffraction=0; 00779 Uzhi_Mx2=1.; 00780 */ 00781 // Generate pt 00782 // if (whilecount++ >= 500 && (whilecount%100)==0) 00783 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 00784 // << ", loop count/ maxPtSquare : " 00785 // << whilecount << " / " << maxPtSquare << G4endl; 00786 00787 // whilecount++; 00788 if (whilecount > 1000 ) 00789 { 00790 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 00791 target->SetStatus(2); return false; // Ignore this interaction 00792 }; 00793 00794 // --------------- Check that the interaction is possible ----------- 00795 ProjMassT2=ProjectileDiffStateMinMass2; 00796 ProjMassT =ProjectileDiffStateMinMass; 00797 00798 TargMassT2=M0target2; 00799 TargMassT =M0target; 00800 //G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl; 00801 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 00802 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 00803 /4./S; 00804 00805 //G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl; 00806 if(PZcms2 < 0 ) 00807 { 00808 target->SetStatus(2); 00809 return false; 00810 } 00811 maxPtSquare=PZcms2; 00812 00813 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 00814 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 00815 00816 ProjMassT2=ProjectileDiffStateMinMass2+Pt2; 00817 ProjMassT =std::sqrt(ProjMassT2); 00818 00819 TargMassT2=M0target2+Pt2; 00820 TargMassT =std::sqrt(TargMassT2); 00821 00822 //G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl; 00823 00824 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 00825 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 00826 /4./S; 00827 //G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl; 00828 if(PZcms2 < 0 ) continue; 00829 PZcms =std::sqrt(PZcms2); 00830 00831 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 00832 PMinusMax=SqrtS-TargMassT; 00833 00834 PMinusNew=ChooseP(PMinusMin, PMinusMax); 00835 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax)); 00836 //PMinusNew=1./sqr(1./std::sqrt(PMinusMin)-G4UniformRand()*(1./std::sqrt(PMinusMin)-1./std::sqrt(PMinusMax))); 00837 00838 TMinusNew=SqrtS-PMinusNew; 00839 Qminus=Ptarget.minus()-TMinusNew; 00840 TPlusNew=TargMassT2/TMinusNew; 00841 Qplus=Ptarget.plus()-TPlusNew; 00842 00843 Qmomentum.setPz( (Qplus-Qminus)/2 ); 00844 Qmomentum.setE( (Qplus+Qminus)/2 ); 00845 00846 } while ((Pprojectile+Qmomentum).mag2() < ProjectileDiffStateMinMass2); //|| 00847 //Repeat the sampling because there was not any excitation 00848 //((Ptarget -Qmomentum).mag2() < M0target2 )) ); 00849 projectile->SetStatus(1*projectile->GetStatus()); // VU 10.04.2012 00850 } 00851 else 00852 { // -------------- Target diffraction ---------------- 00853 00854 //G4cout<<"Target diffraction"<<G4endl; 00855 do { 00856 /* 00857 Uzhi_projectilediffraction=0; 00858 Uzhi_targetdiffraction=1; 00859 Uzhi_Mx2=1.; 00860 */ 00861 // Generate pt 00862 // if (whilecount++ >= 500 && (whilecount%100)==0) 00863 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 00864 // << ", loop count/ maxPtSquare : " 00865 // << whilecount << " / " << maxPtSquare << G4endl; 00866 00867 // whilecount++; 00868 if (whilecount > 1000 ) 00869 { 00870 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 00871 target->SetStatus(2); return false; // Ignore this interaction 00872 }; 00873 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl; 00874 // --------------- Check that the interaction is possible ----------- 00875 ProjMassT2=M0projectile2; 00876 ProjMassT =M0projectile; 00877 00878 TargMassT2=TargetDiffStateMinMass2; 00879 TargMassT =TargetDiffStateMinMass; 00880 00881 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 00882 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 00883 /4./S; 00884 00885 //G4cout<<"PZcms2 TrD <0 "<<PZcms2<<" return"<<G4endl; 00886 if(PZcms2 < 0 ) 00887 { 00888 target->SetStatus(2); 00889 return false; 00890 } 00891 maxPtSquare=PZcms2; 00892 00893 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 00894 00895 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl; 00896 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 00897 00898 ProjMassT2=M0projectile2+Pt2; 00899 ProjMassT =std::sqrt(ProjMassT2); 00900 00901 TargMassT2=TargetDiffStateMinMass2+Pt2; 00902 TargMassT =std::sqrt(TargMassT2); 00903 00904 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 00905 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 00906 /4./S; 00907 00908 //G4cout<<"PZcms2 <0 "<<PZcms2<<" continue"<<G4endl; 00909 if(PZcms2 < 0 ) continue; 00910 PZcms =std::sqrt(PZcms2); 00911 00912 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 00913 TPlusMax=SqrtS-ProjMassT; 00914 00915 TPlusNew=ChooseP(TPlusMin, TPlusMax); 00916 //TPlusNew=1./sqr(1./std::sqrt(TPlusMin)-G4UniformRand()*(1./std::sqrt(TPlusMin)-1./std::sqrt(TPlusMax))); 00917 00918 //TPlusNew=TPlusMax; 00919 00920 PPlusNew=SqrtS-TPlusNew; 00921 Qplus=PPlusNew-Pprojectile.plus(); 00922 PMinusNew=ProjMassT2/PPlusNew; 00923 Qminus=PMinusNew-Pprojectile.minus(); 00924 00925 Qmomentum.setPz( (Qplus-Qminus)/2 ); 00926 Qmomentum.setE( (Qplus+Qminus)/2 ); 00927 00928 /* 00929 G4cout<<(Pprojectile+Qmomentum).mag()<<" "<<M0projectile<<G4endl; 00930 G4bool First=(Pprojectile+Qmomentum).mag2() < M0projectile2; 00931 G4cout<<First<<G4endl; 00932 00933 G4cout<<(Ptarget -Qmomentum).mag()<<" "<<TargetDiffStateMinMass<<" "<<TargetDiffStateMinMass2<<G4endl; 00934 G4bool Seco=(Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2; 00935 G4cout<<Seco<<G4endl; 00936 */ 00937 00938 } while ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2); 00939 // Repeat the sampling because there was not any excitation 00940 // (((Pprojectile+Qmomentum).mag2() < M0projectile2 ) || //No without excitation 00941 // ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2)) ); 00942 //G4cout<<"Go out"<<G4endl; 00943 target->SetStatus(1*target->GetStatus()); // VU 10.04.2012 00944 } // End of if(G4UniformRand() < ProbProjectileDiffraction) 00945 } 00946 else //----------- Non-diffraction process ------------ 00947 { 00948 00949 //G4cout<<"Non-diffraction process"<<G4endl; 00950 do { 00951 /* 00952 Uzhi_projectilediffraction=0; 00953 Uzhi_targetdiffraction=0; 00954 Uzhi_Mx2=1.; 00955 */ 00956 // Generate pt 00957 // if (whilecount++ >= 500 && (whilecount%100)==0) 00958 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 00959 // << ", loop count/ maxPtSquare : " 00960 // << whilecount << " / " << maxPtSquare << G4endl; 00961 00962 // whilecount++; 00963 if (whilecount > 1000 ) 00964 { 00965 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 00966 target->SetStatus(2); return false; // Ignore this interaction 00967 }; 00968 // --------------- Check that the interaction is possible ----------- 00969 ProjMassT2=ProjectileNonDiffStateMinMass2; 00970 ProjMassT =ProjectileNonDiffStateMinMass; 00971 00972 TargMassT2=TargetNonDiffStateMinMass2; 00973 TargMassT =TargetNonDiffStateMinMass; 00974 00975 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 00976 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 00977 /4./S; 00978 00979 if(PZcms2 < 0 ) 00980 { 00981 target->SetStatus(2); 00982 return false; 00983 } 00984 maxPtSquare=PZcms2; 00985 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 00986 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 00987 00988 ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2; 00989 ProjMassT =std::sqrt(ProjMassT2); 00990 00991 TargMassT2=TargetNonDiffStateMinMass2+Pt2; 00992 TargMassT =std::sqrt(TargMassT2); 00993 00994 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 00995 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 00996 /4./S; 00997 //G4cout<<"PZcms2 ND"<<PZcms2<<G4endl; 00998 00999 if(PZcms2 < 0 ) continue; 01000 PZcms =std::sqrt(PZcms2); 01001 01002 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 01003 PMinusMax=SqrtS-TargMassT; 01004 01005 if(G4UniformRand() < ProbLogDistr) // Uzhi 25.04.2012 01006 { PMinusNew=ChooseP(PMinusMin, PMinusMax);} // 12.06.11 01007 else {PMinusNew=(PMinusMax-PMinusMin)*G4UniformRand() + PMinusMin;} 01008 Qminus=PMinusNew-Pprojectile.minus(); 01009 01010 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 01011 TPlusMax=SqrtS-PMinusNew; // Why is it closed??? 01012 // TPlusMax=SqrtS-ProjMassT; 01013 01014 if(G4UniformRand() < 0.5) //ProbLogDistr) // Uzhi 29.05.2012 0.5) 01015 { TPlusNew=ChooseP(TPlusMin, TPlusMax);} // 12.06.11 01016 else {TPlusNew=(TPlusMax-TPlusMin)*G4UniformRand() +TPlusMin;} 01017 01018 Qplus=-(TPlusNew-Ptarget.plus()); 01019 01020 Qmomentum.setPz( (Qplus-Qminus)/2 ); 01021 Qmomentum.setE( (Qplus+Qminus)/2 ); 01022 /* 01023 G4cout<<(Pprojectile+Qmomentum).mag2()<<" "<<ProjectileNonDiffStateMinMass2<<G4endl; 01024 G4cout<<(Ptarget -Qmomentum).mag2()<<" "<<TargetNonDiffStateMinMass2<<G4endl; 01025 G4int Uzhi; G4cin>>Uzhi; 01026 */ 01027 } while ( 01028 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction 01029 ((Ptarget -Qmomentum).mag2() < TargetNonDiffStateMinMass2 )); 01030 01031 projectile->SetStatus(0*projectile->GetStatus()); // VU 10.04.2012 01032 target->SetStatus(0*target->GetStatus()); // VU 10.04.2012 01033 } 01034 01035 Pprojectile += Qmomentum; 01036 Ptarget -= Qmomentum; 01037 01038 //G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl; 01039 01040 // Transform back and update SplitableHadron Participant. 01041 Pprojectile.transform(toLab); 01042 Ptarget.transform(toLab); 01043 01044 // Calculation of the creation time --------------------- 01045 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); 01046 projectile->SetPosition(target->GetPosition()); 01047 // Creation time and position of target nucleon were determined at 01048 // ReggeonCascade() of G4FTFModel 01049 // ------------------------------------------------------ 01050 /* 01051 if(Uzhi_projectilediffraction != 0) 01052 {Uzhi_Mx2=Pprojectile.mag2(); Uzhi_modT=(target->Get4Momentum()-Ptarget).mag2();} 01053 01054 if(Uzhi_targetdiffraction != 0) 01055 {Uzhi_Mx2=Ptarget.mag2(); Uzhi_modT=(projectile->Get4Momentum()-Pprojectile).mag2();} 01056 01057 if(Uzhi_QE!= 0) 01058 { 01059 Uzhi_projectilediffraction=0; 01060 Uzhi_targetdiffraction =0; 01061 Uzhi_Mx2 =1.; 01062 } 01063 */ 01064 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl; 01065 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl; 01066 projectile->Set4Momentum(Pprojectile); 01067 target->Set4Momentum(Ptarget); 01068 01069 projectile->IncrementCollisionCount(1); 01070 target->IncrementCollisionCount(1); 01071 01072 return true; 01073 }