- Timestamp:
- Nov 25, 2009, 5:13:58 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/parton_string/diffraction/src
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc
r962 r1196 25 25 // 26 26 // 27 // $Id: G4DiffractiveExcitation.cc,v 1. 7 2008/12/18 13:01:58 gunterExp $27 // $Id: G4DiffractiveExcitation.cc,v 1.16 2009/10/06 10:10:36 vuzhinsk Exp $ 28 28 // ------------------------------------------------------------ 29 29 // GEANT 4 class implemetation file … … 48 48 49 49 #include "G4DiffractiveExcitation.hh" 50 #include "G4FTFParameters.hh" 51 #include "G4ElasticHNScattering.hh" 52 50 53 #include "G4LorentzRotation.hh" 54 #include "G4RotationMatrix.hh" 51 55 #include "G4ThreeVector.hh" 52 #include "G4ParticleDefinition.hh" 56 #include "G4ParticleDefinition.hh" 53 57 #include "G4VSplitableHadron.hh" 54 58 #include "G4ExcitedString.hh" 55 #include "G4FTFParameters.hh" // Uzhi 19.04.08 59 #include "G4ParticleTable.hh" 60 #include "G4Neutron.hh" 61 #include "G4ParticleDefinition.hh" 62 56 63 //#include "G4ios.hh" 57 64 //#include "UZHI_diffraction.hh" … … 65 72 ExciteParticipants(G4VSplitableHadron *projectile, 66 73 G4VSplitableHadron *target, 67 G4FTFParameters *theParameters) const 74 G4FTFParameters *theParameters, 75 G4ElasticHNScattering *theElastic) const // Uzhi 03.09.09 68 76 { 69 G4bool PutOnMassShell=0;70 71 77 // -------------------- Projectile parameters ----------------------- 72 73 78 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 79 //G4cout<<"Excite P "<<Pprojectile<<G4endl; 80 81 if(Pprojectile.z() < 0.) 82 { 83 target->SetStatus(2); 84 return false; 85 } 86 87 G4double ProjectileRapidity = Pprojectile.rapidity(); 88 89 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding(); 90 // G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition(); 91 G4int absProjectilePDGcode=std::abs(ProjectilePDGcode); 92 93 G4bool PutOnMassShell(false); 74 94 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation 75 95 G4double M0projectile = Pprojectile.mag(); // Without de-excitation 76 /* 77 G4cout<<"ExciteParticipants-------------------"<<G4endl; 78 G4cout<<"Mom "<<Pprojectile<<" mass "<<M0projectile<<G4endl; 79 */ 96 80 97 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 81 98 { 82 PutOnMassShell= 1;99 PutOnMassShell=true; 83 100 M0projectile=projectile->GetDefinition()->GetPDGMass(); 84 101 } 85 102 86 103 G4double M0projectile2 = M0projectile * M0projectile; 87 88 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();89 G4int absPDGcode=std::abs(PDGcode);90 104 91 105 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass(); 92 106 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass(); 93 107 G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff(); 94 /* 95 G4cout<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<" "<<ProbProjectileDiffraction<<G4endl; 96 */ 97 // -------------------- Target paraExciteParticipantsmeters ------------------------- 108 109 // -------------------- Target parameters ------------------------- 110 G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding(); 111 G4int absTargetPDGcode=std::abs(TargetPDGcode); 112 98 113 G4LorentzVector Ptarget=target->Get4Momentum(); 114 99 115 G4double M0target = Ptarget.mag(); 100 116 101 //G4cout<<"Mom "<<Ptarget<<" mass "<<M0target<<G4endl; 117 //G4cout<<"Excite T "<<Ptarget<<G4endl; 118 //G4cout<<"PDGs "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl; 119 //G4cout<<"Masses "<<M0projectile<<" "<<M0target<<G4endl; 120 121 G4double TargetRapidity = Ptarget.rapidity(); 102 122 103 123 if(M0target < target->GetDefinition()->GetPDGMass()) 104 124 { 105 PutOnMassShell= 1;125 PutOnMassShell=true; 106 126 M0target=target->GetDefinition()->GetPDGMass(); 107 127 } 108 128 109 G4double M0target2 = M0target * M0target; //Ptarget.mag2();110 // for AA-inter.129 G4double M0target2 = M0target * M0target; 130 111 131 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass(); 112 132 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass(); 113 133 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff(); 114 /* 115 G4cout<<TargetDiffStateMinMass<<" "<<TargetNonDiffStateMinMass<<" "<<ProbTargetDiffraction<<G4endl; 116 */ 134 117 135 G4double AveragePt2=theParameters->GetAveragePt2(); 136 137 G4double ProbOfDiffraction=ProbProjectileDiffraction + 138 ProbTargetDiffraction; 139 140 G4double SumMasses=M0projectile+M0target+200.*MeV; 118 141 119 142 // Kinematical properties of the interactions -------------- … … 122 145 G4double S=Psum.mag2(); 123 146 124 //G4cout<<" sqrt(s) "<<std::sqrt(S)<<G4endl;125 126 // ------------------------------------------------------------------127 128 //ProbProjectileDiffraction=1.;129 //ProbTargetDiffraction =1.;130 G4double ProbOfDiffraction=ProbProjectileDiffraction +131 ProbTargetDiffraction;132 133 if(ProbOfDiffraction!=0.)134 {135 ProbProjectileDiffraction/=ProbOfDiffraction;136 }137 else138 {139 ProbProjectileDiffraction=0.;140 }141 // ProbTargetDiffraction /=ProbOfDiffraction;142 143 //G4cout<<"ProbOfDiffraction "<<ProbOfDiffraction<<"ProbProjectileDiffraction "<<ProbProjectileDiffraction<<G4endl; // Vova144 145 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass *146 ProjectileDiffStateMinMass;147 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *148 ProjectileNonDiffStateMinMass;149 150 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass *151 TargetDiffStateMinMass;152 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass *153 TargetNonDiffStateMinMass;154 155 147 // Transform momenta to cms and then rotate parallel to z axis; 156 157 // G4LorentzVector Psum;158 // Psum=Pprojectile+Ptarget;159 160 148 G4LorentzRotation toCms(-1*Psum.boostVector()); 161 149 162 150 G4LorentzVector Ptmp=toCms*Pprojectile; 163 151 if ( Ptmp.pz() <= 0. ) 164 165 // "String" moving backwards in CMS, abort collision !! 166 //G4cout << " abort Collision!! " << G4endl;167 168 152 { 153 target->SetStatus(2); 154 // "String" moving backwards in CMS, abort collision !! 155 return false; 156 } 169 157 170 158 toCms.rotateZ(-1*Ptmp.phi()); … … 176 164 Ptarget.transform(toCms); 177 165 178 G4double Pt2;179 G4double ProjMassT2, ProjMassT;180 G4double TargMassT2, TargMassT;181 166 G4double PZcms2, PZcms; 182 G4double PMinusMin, PMinusMax; 183 // G4double PPlusMin , PPlusMax; 184 G4double TPlusMin , TPlusMax; 185 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew; 186 187 // G4double S=Psum.mag2(); 167 188 168 G4double SqrtS=std::sqrt(S); 189 190 if(absP DGcode > 1000 && SqrtS < 2200*MeV)191 { return false;}// The model cannot work for169 //G4cout<<" SqrtS "<<SqrtS<<" 2200 "<<G4endl; 170 if(absProjectilePDGcode > 1000 && SqrtS < 2300*MeV && SqrtS < SumMasses) 171 {target->SetStatus(2); return false;} // The model cannot work for 192 172 // p+p-interactions 193 // at Plab < 1.3 GeV/c. 194 195 if(( absPDGcode == 211 || PDGcode == 111) && SqrtS < 1600*MeV) 196 {return false;} // The model cannot work for 173 // at Plab < 1.62 GeV/c. 174 175 if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) && 176 (SqrtS < 1600*MeV) && (SqrtS < SumMasses)) 177 {target->SetStatus(2); return false;} // The model cannot work for 197 178 // Pi+p-interactions 198 179 // at Plab < 1. GeV/c. 199 180 200 if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV) 201 {return false;} // The model cannot work for 181 if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) || 182 (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) || 183 (absProjectilePDGcode == 310)) && (SqrtS < 1600*MeV) && (SqrtS < SumMasses)) 184 {target->SetStatus(2); return false;} // The model cannot work for 202 185 // K+p-interactions 203 186 // at Plab < ??? GeV/c. ??? … … 208 191 209 192 if(PZcms2 < 0) 210 {return false;} // It can be in an interaction with off-shell nuclear nucleon 193 {target->SetStatus(2); return false;} // It can be in an interaction with 194 // off-shell nuclear nucleon 211 195 212 196 PZcms = std::sqrt(PZcms2); … … 236 220 237 221 G4double maxPtSquare; // = PZcms2; 238 /* 239 G4cout << "Pprojectile aft boost : " << Pprojectile <<" "<<Pprojectile.mag()<< G4endl; 240 G4cout << "Ptarget aft boost : " << Ptarget <<" "<<Ptarget.mag()<< G4endl; 241 G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl; 242 G4cout << " Projectile Xplus / Xminus : " << 243 Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl; 244 G4cout << " Target Xplus / Xminus : " << Ptarget.plus() << " / " << Ptarget.minus() << G4endl; 245 G4cout<<"maxPtSquare "<<maxPtSquare<<G4endl; 246 */ 222 223 // Charge exchange can be possible for baryons ----------------- 224 225 //G4cout<<1./std::exp(-2.0*(ProjectileRapidity - TargetRapidity))<<G4endl; 226 //G4int Uzhi; G4cin>>Uzhi; 227 // Getting the values needed for exchange ---------------------- 228 G4double MagQuarkExchange =theParameters->GetMagQuarkExchange(); 229 G4double SlopeQuarkExchange =theParameters->GetSlopeQuarkExchange(); 230 G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange(); 231 232 // G4double NucleonMass= 233 (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass(); 234 G4double DeltaMass= 235 (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass(); 236 237 // Check for possible quark excjane ----------------------------------- 238 if(G4UniformRand() < MagQuarkExchange* 239 std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity))) 240 { 241 242 //G4cout<<"Exchange "<<G4endl; 243 244 G4int NewProjCode(0), NewTargCode(0); 245 246 //G4bool StandardExcitation(false); // ================================= 247 G4int ProjQ1(0), ProjQ2(0), ProjQ3(0); 248 249 // Projectile unpacking -------------------------- 250 if(absProjectilePDGcode < 1000 ) 251 { // projectile is meson ----------------- 252 UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2); 253 } else 254 { // projectile is baryon ---------------- 255 UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3); 256 } // End of the hadron's unpacking ---------- 257 258 // Target unpacking ------------------------------ 259 G4int TargQ1(0), TargQ2(0), TargQ3(0); 260 UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3); 261 262 // Sampling of exchanged quarks ------------------- 263 G4int ProjExchangeQ(0); 264 G4int TargExchangeQ(0); 265 266 //G4cout<<"Targ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl; 267 if(absProjectilePDGcode < 1000 ) 268 { // projectile is meson ----------------- 269 //G4cout<<"Proj Q1 Q2 "<<ProjQ1<<" "<<ProjQ2<<G4endl; 270 271 if(ProjQ1 > 0 ) // ProjQ1 is quark 272 { 273 ProjExchangeQ = ProjQ1; 274 if(ProjExchangeQ != TargQ1) 275 { 276 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ; 277 } else 278 if(ProjExchangeQ != TargQ2) 279 { 280 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ; 281 } else 282 { 283 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ; 284 } 285 //G4cout<<" Pr Tr "<<ProjQ1<<" "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl; 286 } else // ProjQ2 is quark 287 { 288 ProjExchangeQ = ProjQ2; 289 if(ProjExchangeQ != TargQ1) 290 { 291 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ; 292 } else 293 if(ProjExchangeQ != TargQ2) 294 { 295 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ; 296 } else 297 { 298 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ; 299 } 300 301 } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark 302 303 G4int aProjQ1=std::abs(ProjQ1); 304 G4int aProjQ2=std::abs(ProjQ2); 305 if(aProjQ1 == aProjQ2) {NewProjCode = 111;} // Pi0-meson 306 else // |ProjQ1| # |ProjQ2| 307 { 308 if(aProjQ1 > aProjQ2) {NewProjCode = aProjQ1*100+aProjQ2*10+1;} 309 else {NewProjCode = aProjQ2*100+aProjQ1*10+1;} 310 } 311 312 // projectile->SetDefinition( 313 // G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); 314 315 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); 316 317 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) && 318 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar 319 320 else if( target->GetDefinition()->GetPDGiIsospin() == 3 ) //Delta was the target 321 { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;} //Save Delta isobar 322 else {} // De-excite initial Delta isobar 323 } 324 325 else if((G4UniformRand() < DeltaProbAtQuarkExchange) && //Nucleon was the target 326 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar 327 else {} //Save initial nucleon 328 329 target->SetDefinition( 330 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 331 //G4cout<<"New PDGs "<<NewProjCode<<" "<<NewTargCode<<G4endl; 332 //G4int Uzhi; G4cin>>Uzhi; 333 //} 334 } else 335 { // projectile is baryon ---------------- 336 G4double Same=0.; //0.3; //0.5; 337 G4bool ProjDeltaHasCreated(false); 338 G4bool TargDeltaHasCreated(false); 339 340 G4double Ksi=G4UniformRand(); 341 if(G4UniformRand() < 0.5) // Sampling exchange quark from proj. or targ. 342 { // Sampling exchanged quark from the projectile --- 343 //G4cout<<"Proj Exc"<<G4endl; 344 if( Ksi < 0.333333 ) 345 {ProjExchangeQ = ProjQ1;} 346 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 347 {ProjExchangeQ = ProjQ2;} 348 else 349 {ProjExchangeQ = ProjQ3;} 350 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 351 352 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) // Vova 353 { 354 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 355 } else 356 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) // Vova 357 { 358 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 359 } else 360 { 361 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 362 } 363 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 364 if( Ksi < 0.333333 ) 365 {ProjQ1=ProjExchangeQ;} 366 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 367 {ProjQ2=ProjExchangeQ;} 368 else 369 {ProjQ3=ProjExchangeQ;} 370 371 /* 372 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // ***************************** 373 374 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;} 375 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta 376 { if(G4UniformRand() > DeltaProbAtQuarkExchange) 377 {NewProjCode +=2; ProjDeltaHasCreated=true;} 378 else {NewProjCode +=0; ProjDeltaHasCreated=false;} 379 } 380 else // Projectile was Nucleon 381 { 382 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) 383 {NewProjCode +=2; ProjDeltaHasCreated=true;} 384 else {NewProjCode +=0; ProjDeltaHasCreated=false;} 385 } 386 387 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // ***************************** 388 389 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 390 else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta 391 { if(G4UniformRand() > DeltaProbAtQuarkExchange) 392 {NewTargCode +=2; TargDeltaHasCreated=true;} 393 else {NewTargCode +=0; TargDeltaHasCreated=false;} 394 } 395 else // Target was Nucleon 396 { 397 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) 398 {NewTargCode +=2; TargDeltaHasCreated=true;} 399 else {NewTargCode +=0; TargDeltaHasCreated=false;} 400 } 401 402 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode)) 403 { // Nothing was changed! It is not right!? 404 } 405 */ 406 /*----------------------------------------------------------------------------------------- 407 if(!ProjDeltaHasCreated && !TargDeltaHasCreated) // No Deltas were created then 408 { if( G4UniformRand() < 0.5) {NewProjCode +=2; ProjDeltaHasCreated = true;} 409 else {NewTargCode +=2; TargDeltaHasCreated = true;} 410 } 411 412 if(!(ProjDeltaHasCreated && TargDeltaHasCreated)) 413 { 414 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.)) 415 { 416 G4cout<<"2 Deltas "<<G4endl; 417 if(!ProjDeltaHasCreated) 418 {NewProjCode +=2; ProjDeltaHasCreated = true;} 419 else {NewTargCode +=2; TargDeltaHasCreated = true;} // For Delta isobar 420 } 421 } 422 */ 423 } else 424 { // Sampling exchanged quark from the target ------- 425 //G4cout<<"Targ Exc"<<G4endl; 426 if( Ksi < 0.333333 ) 427 {TargExchangeQ = TargQ1;} 428 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 429 {TargExchangeQ = TargQ2;} 430 else 431 {TargExchangeQ = TargQ3;} 432 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl; 433 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) // Vova 434 { 435 ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 436 } else 437 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) // Vova 438 { 439 ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 440 } else 441 { 442 ProjExchangeQ = ProjQ3; ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 443 } 444 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl; 445 if( Ksi < 0.333333 ) 446 {TargQ1=TargExchangeQ;} 447 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 448 {TargQ2=TargExchangeQ;} 449 else 450 {TargQ3=TargExchangeQ;} 451 /* 452 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); 453 if( (ProjQ1 == ProjQ2) && (ProjQ1 == ProjQ3) ) 454 {NewProjCode +=2; ProjDeltaHasCreated = true;} 455 456 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); 457 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) ) 458 {NewTargCode +=2; TargDeltaHasCreated = true;} 459 460 if(!ProjDeltaHasCreated && !TargDeltaHasCreated) 461 {NewTargCode +=2; TargDeltaHasCreated = true;} 462 463 if(!(ProjDeltaHasCreated && TargDeltaHasCreated)) 464 { 465 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.)) 466 { 467 G4cout<<"2 Deltas "<<G4endl; 468 if(!ProjDeltaHasCreated) 469 {NewProjCode +=2; ProjDeltaHasCreated = true;} 470 else {NewTargCode +=2; TargDeltaHasCreated = true;} 471 } 472 } 473 */ 474 } // End of sampling baryon 475 476 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // ***************************** 477 478 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;} 479 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta 480 { if(G4UniformRand() > DeltaProbAtQuarkExchange) 481 {NewProjCode +=2; ProjDeltaHasCreated=true;} 482 else {NewProjCode +=0; ProjDeltaHasCreated=false;} 483 } 484 else // Projectile was Nucleon 485 { 486 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) 487 {NewProjCode +=2; ProjDeltaHasCreated=true;} 488 else {NewProjCode +=0; ProjDeltaHasCreated=false;} 489 } 490 491 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // ***************************** 492 493 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 494 else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta 495 { if(G4UniformRand() > DeltaProbAtQuarkExchange) 496 {NewTargCode +=2; TargDeltaHasCreated=true;} 497 else {NewTargCode +=0; TargDeltaHasCreated=false;} 498 } 499 else // Target was Nucleon 500 { 501 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) 502 {NewTargCode +=2; TargDeltaHasCreated=true;} 503 else {NewTargCode +=0; TargDeltaHasCreated=false;} 504 } 505 506 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode)) 507 { // Nothing was changed! It is not right!? 508 } 509 // Forming baryons -------------------------------------------------- 510 511 // if( ProjDeltaHasCreated && TargDeltaHasCreated ) {StandardExcitation = true;} 512 513 //G4cout<<"New PDG "<<NewProjCode<<" "<<NewTargCode<<G4endl; 514 //G4int Uzhi; G4cin>>Uzhi; 515 } // End of if projectile is baryon --------------------------- 516 517 518 // If we assume that final state hadrons after the charge exchange will be 519 // in the ground states, we have to put ---------------------------------- 520 if(M0projectile < 521 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass()) 522 { 523 M0projectile= 524 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass(); 525 M0projectile2 = M0projectile * M0projectile; 526 } 527 528 if(M0target < 529 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass()) 530 { 531 M0target= 532 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass(); 533 M0target2 = M0target * M0target; 534 } 535 536 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- 537 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) 538 /4./S; 539 //G4cout<<"New Masses "<<M0projectile<<" "<<M0target<<G4endl; 540 //G4cout<<"PZcms2 "<<PZcms2<<G4endl; 541 542 if(PZcms2 < 0) {return false;} // It can be if energy is not sufficient for Delta 543 //---------------------------------------------------------- 544 projectile->SetDefinition( 545 G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); 546 547 target->SetDefinition( 548 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 549 //---------------------------------------------------------- 550 PZcms = std::sqrt(PZcms2); 551 552 Pprojectile.setPz( PZcms); 553 Pprojectile.setE(std::sqrt(M0projectile2+PZcms2)); 554 555 Ptarget.setPz( -PZcms); 556 Ptarget.setE(std::sqrt(M0target2+PZcms2)); 557 558 //G4cout<<"Proj "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl; 559 //G4cout<<"Targ "<<Ptarget<<" "<<Ptarget.mag()<<G4endl; 560 561 // if(!StandardExcitation) 562 { 563 Pprojectile.transform(toLab); 564 Ptarget.transform(toLab); 565 566 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); 567 projectile->SetPosition(target->GetPosition()); 568 569 projectile->Set4Momentum(Pprojectile); 570 target->Set4Momentum(Ptarget); 571 572 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters); 573 //G4cout<<"1 Delta result "<<Result<<"**********"<<G4endl; 574 return Result; 575 } 576 /* 577 else 578 { 579 if(M0projectile > ProjectileDiffStateMinMass) 580 { ProjectileDiffStateMinMass +=200.*MeV; ProjectileNonDiffStateMinMass +=200.*MeV;} 581 582 if(M0target > TargetDiffStateMinMass) 583 { TargetDiffStateMinMass +=200.*MeV; TargetNonDiffStateMinMass +=200.*MeV;} 584 585 ProbOfDiffraction = 1.; 586 ProbProjectileDiffraction =0.5; 587 ProbTargetDiffraction =0.5; 588 G4cout<<"Exc DD "<<M0projectile<<" "<<M0target<<" --------------------"<<G4endl; 589 // After that standard FTF excitation 590 } 591 */ 592 } // End of charge exchange part ------------------------------ 593 594 // ------------------------------------------------------------------ 595 // G4double ProbOfDiffraction=ProbProjectileDiffraction + 596 // ProbTargetDiffraction; 597 //G4cout<<ProbOfDiffraction<<" "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<G4endl; 598 //G4int Uzhi; G4cin>>Uzhi; 599 if(ProbOfDiffraction!=0.) 600 { 601 ProbProjectileDiffraction/=ProbOfDiffraction; 602 } 603 else 604 { 605 ProbProjectileDiffraction=0.; 606 } 607 608 //ProbOfDiffraction=1.; //0.5; // Vova 609 610 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass * 611 ProjectileDiffStateMinMass; 612 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass * 613 ProjectileNonDiffStateMinMass; 614 615 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass * 616 TargetDiffStateMinMass; 617 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass * 618 TargetNonDiffStateMinMass; 619 620 621 G4double Pt2; 622 G4double ProjMassT2, ProjMassT; 623 G4double TargMassT2, TargMassT; 624 G4double PMinusMin, PMinusMax; 625 // G4double PPlusMin , PPlusMax; 626 G4double TPlusMin , TPlusMax; 627 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew; 628 247 629 G4LorentzVector Qmomentum; 248 630 G4double Qminus, Qplus; 249 631 250 632 G4int whilecount=0; 251 // Choose a process633 // Choose a process --------------------------- 252 634 253 635 if(G4UniformRand() < ProbOfDiffraction) … … 255 637 if(G4UniformRand() < ProbProjectileDiffraction) 256 638 { //-------- projectile diffraction --------------- 257 //G4cout<<" Projectile diffraction"<<G4endl; 258 //Uzhi_projectilediffraction++; 639 //G4cout<<"Proj difr"<<G4endl; 259 640 do { 260 641 // Generate pt … … 266 647 { 267 648 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 268 return false; // Ignore this interaction649 target->SetStatus(2); return false; // Ignore this interaction 269 650 }; 270 651 // --------------- Check that the interaction is possible ----------- … … 278 659 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 279 660 /4./S; 280 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 281 282 if(PZcms2 < 0 ) 283 { 284 /* 285 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 286 G4int Uzhi; G4cin>>Uzhi; 287 */ 288 return false; 289 }; 661 //G4cout<<"PZcms2 < 0 false "<<PZcms2<<G4endl; 662 if(PZcms2 < 0 ) 663 { 664 target->SetStatus(2); 665 return false; 666 } 290 667 maxPtSquare=PZcms2; 291 668 … … 302 679 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 303 680 /4./S; 304 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 305 306 // if(PZcms2 < 0 ) {PZcms2=0;}; 681 307 682 if(PZcms2 < 0 ) continue; 308 683 PZcms =std::sqrt(PZcms2); … … 310 685 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 311 686 PMinusMax=SqrtS-TargMassT; 312 //G4cout<<" SqrtS P+mim max "<<SqrtS<<" "<<PMinusMin<<" "<<PMinusMax<<G4endl;313 687 314 688 PMinusNew=ChooseP(PMinusMin, PMinusMax); … … 328 702 else 329 703 { // -------------- Target diffraction ---------------- 330 //G4cout<<" Target difraction"<<G4endl; 331 //Uzhi_targetdiffraction++; 704 do { 705 //G4cout<<"Targ difr"<<G4endl; 706 // Generate pt 707 // if (whilecount++ >= 500 && (whilecount%100)==0) 708 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 709 // << ", loop count/ maxPtSquare : " 710 // << whilecount << " / " << maxPtSquare << G4endl; 711 if (whilecount > 1000 ) 712 { 713 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 714 target->SetStatus(2); return false; // Ignore this interaction 715 }; 716 // --------------- Check that the interaction is possible ----------- 717 ProjMassT2=M0projectile2; 718 ProjMassT =M0projectile; 719 720 TargMassT2=TargetDiffStateMinMass2; 721 TargMassT =TargetDiffStateMinMass; 722 723 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 724 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 725 /4./S; 726 727 if(PZcms2 < 0 ) 728 { 729 target->SetStatus(2); 730 return false; 731 } 732 maxPtSquare=PZcms2; 733 734 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 735 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 736 737 ProjMassT2=M0projectile2+Pt2; 738 ProjMassT =std::sqrt(ProjMassT2); 739 740 TargMassT2=TargetDiffStateMinMass2+Pt2; 741 TargMassT =std::sqrt(TargMassT2); 742 743 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 744 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 745 /4./S; 746 747 if(PZcms2 < 0 ) continue; 748 PZcms =std::sqrt(PZcms2); 749 750 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 751 TPlusMax=SqrtS-ProjMassT; 752 753 TPlusNew=ChooseP(TPlusMin, TPlusMax); 754 755 //TPlusNew=TPlusMax; 756 757 PPlusNew=SqrtS-TPlusNew; 758 Qplus=PPlusNew-Pprojectile.plus(); 759 PMinusNew=ProjMassT2/PPlusNew; 760 Qminus=PMinusNew-Pprojectile.minus(); 761 762 Qmomentum.setPz( (Qplus-Qminus)/2 ); 763 Qmomentum.setE( (Qplus+Qminus)/2 ); 764 765 } while ( 766 ((Pprojectile+Qmomentum).mag2() < M0projectile2 ) || //No without excitation 767 ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2)); 768 } 769 } 770 else //----------- Non-diffraction process ------------ 771 { 772 //G4cout<<"Non difr"<<G4endl; 332 773 do { 333 774 // Generate pt … … 339 780 { 340 781 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 341 return false; // Ignore this interaction 342 }; 343 // --------------- Check that the interaction is possible ----------- 344 ProjMassT2=M0projectile2; 345 ProjMassT =M0projectile; 346 347 TargMassT2=TargetDiffStateMinMass2; 348 TargMassT =TargetDiffStateMinMass; 349 350 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 351 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 352 /4./S; 353 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 354 355 if(PZcms2 < 0 ) 356 { 357 /* 358 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 359 G4int Uzhi; G4cin>>Uzhi; 360 */ 361 return false; 362 }; 363 maxPtSquare=PZcms2; 364 365 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 366 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 367 368 ProjMassT2=M0projectile2+Pt2; 369 ProjMassT =std::sqrt(ProjMassT2); 370 371 TargMassT2=TargetDiffStateMinMass2+Pt2; 372 TargMassT =std::sqrt(TargMassT2); 373 374 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 375 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 376 /4./S; 377 /* 378 if(PZcms2 < 0 ) 379 { 380 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 381 G4int Uzhi; G4cin>>Uzhi; 382 return false; 383 }; 384 */ 385 if(PZcms2 < 0 ) continue; 386 PZcms =std::sqrt(PZcms2); 387 388 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 389 TPlusMax=SqrtS-ProjMassT; 390 391 //G4cout<<" Tmin max "<<TPlusMin<<" "<<TPlusMax<<G4endl; 392 393 TPlusNew=ChooseP(TPlusMin, TPlusMax); 394 395 //TPlusNew=TPlusMax; 396 //G4cout<<"T+new "<<TPlusNew<<G4endl; 397 398 PPlusNew=SqrtS-TPlusNew; 399 Qplus=PPlusNew-Pprojectile.plus(); 400 PMinusNew=ProjMassT2/PPlusNew; 401 Qminus=PMinusNew-Pprojectile.minus(); 402 403 Qmomentum.setPz( (Qplus-Qminus)/2 ); 404 Qmomentum.setE( (Qplus+Qminus)/2 ); 405 406 } while ( 407 ((Pprojectile+Qmomentum).mag2() < M0projectile2 ) || //No without excitation 408 ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2)); 409 } 410 } 411 else //----------- Non-diffraction process ------------ 412 { 413 //G4cout<<" Non-difraction"<<G4endl; 414 do { 415 // Generate pt 416 // if (whilecount++ >= 500 && (whilecount%100)==0) 417 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 418 // << ", loop count/ maxPtSquare : " 419 // << whilecount << " / " << maxPtSquare << G4endl; 420 if (whilecount > 1000 ) 421 { 422 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 423 return false; // Ignore this interaction 782 target->SetStatus(2); return false; // Ignore this interaction 424 783 }; 425 784 // --------------- Check that the interaction is possible ----------- … … 433 792 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 434 793 /4./S; 435 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 436 437 if(PZcms2 < 0 ) 438 { 439 /* 440 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 441 G4int Uzhi; G4cin>>Uzhi; 442 */ 443 return false; 444 }; 794 795 if(PZcms2 < 0 ) 796 { 797 target->SetStatus(2); 798 return false; 799 } 445 800 maxPtSquare=PZcms2; 446 801 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); … … 456 811 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 457 812 /4./S; 458 /* 459 G4cout<<"ProjectileNonDiffStateMinMass2 "<<ProjectileNonDiffStateMinMass2<<G4endl; 460 G4cout<<"TargetNonDiffStateMinMass2 "<<TargetNonDiffStateMinMass2<<G4endl; 461 G4cout<<"Mt "<<ProjMassT<<" "<<TargMassT<<" "<<Pt2<<" "<<PZcms2<<G4endl<<G4endl; 462 */ 463 // if(PZcms2 < 0 ) {PZcms2=0;}; 813 464 814 if(PZcms2 < 0 ) continue; 465 815 PZcms =std::sqrt(PZcms2); … … 469 819 470 820 PMinusNew=ChooseP(PMinusMin, PMinusMax); 471 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));472 473 //G4cout<<"Proj "<<PMinusMin<<" "<<PMinusMax<<" "<<PMinusNew<<G4endl;474 475 //PMinusNew=PMinusMax; //+++++++++++++++++++++++++++++++++++ Vova476 821 477 822 Qminus=PMinusNew-Pprojectile.minus(); 478 823 479 824 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 480 // TPlusMax=SqrtS-PMinusNew; // Vova481 TPlusMax=SqrtS-ProjMassT; // Vova825 // TPlusMax=SqrtS-PMinusNew; 826 TPlusMax=SqrtS-ProjMassT; 482 827 483 828 TPlusNew=ChooseP(TPlusMin, TPlusMax); 484 485 //G4cout<<"Targ "<<TPlusMin<<" "<<TPlusMax<<" "<<TPlusNew<<G4endl;486 //G4cout<<PMinusNew<<" "<<TPlusNew<<G4endl;487 829 488 830 Qplus=-(TPlusNew-Ptarget.plus()); … … 490 832 Qmomentum.setPz( (Qplus-Qminus)/2 ); 491 833 Qmomentum.setE( (Qplus+Qminus)/2 ); 492 /* 493 G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl; 494 G4cout << "pt2" << pt2 << G4endl; 495 G4cout << "Qmomentum " << Qmomentum << G4endl; 496 G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() << 497 " / " << (Ptarget-Qmomentum).mag() << G4endl; // mag() 498 G4cout<<"Mprojectile "<<std::sqrt(M0projectile2)<<G4endl; 499 G4cout<<"Mtarget "<<std::sqrt(M0target2 )<<G4endl; 500 G4cout<<"ProjectileDiffStateMinMass "<<std::sqrt(ProjectileDiffStateMinMass2)<<G4endl; 501 G4cout<<"TargetDiffStateMinMass "<<std::sqrt(TargetDiffStateMinMass2)<<G4endl; 502 */ 834 503 835 } while ( 504 836 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction … … 506 838 } 507 839 508 //G4int Uzhiinp; G4cin>>Uzhiinp; // Vova509 510 840 Pprojectile += Qmomentum; 511 841 Ptarget -= Qmomentum; 512 /* 513 G4cout << "Pprojectile with Q : " << Pprojectile << G4endl; 514 G4cout << "Ptarget with Q : " << Ptarget << G4endl; 515 G4cout << "Target mass " << Ptarget.mag() << G4endl; 516 G4cout << "Projectile mass " << Pprojectile.mag() << G4endl; 517 // 518 //G4cout << "Projectile back: " << toLab * Pprojectile << G4endl; 519 //G4cout << "Target back: " << toLab * Ptarget << G4endl; 520 */ 521 //-------------- Flip if projectale moves in backward direction ------------ 522 //G4bool Flip=Pprojectile.pz()< 0.; 523 842 843 //G4cout<<"Masses "<<Pprojectile.mag()<<" "<<Ptarget.mag()<<G4endl; 844 //G4int Uzhi; G4cin>>Uzhi; 524 845 525 846 // Transform back and update SplitableHadron Participant. … … 527 848 Ptarget.transform(toLab); 528 849 529 // G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;530 //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;531 //G4cout << "Target mass " << Ptarget.mag() << G4endl;532 // G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;533 534 / *535 if(!Flip){ 850 // Calculation of the creation time --------------------- 851 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); 852 projectile->SetPosition(target->GetPosition()); 853 // Creation time and position of target nucleon were determined at 854 // ReggeonCascade() of G4FTFModel 855 // ------------------------------------------------------ 856 536 857 projectile->Set4Momentum(Pprojectile); 537 858 target->Set4Momentum(Ptarget); 538 }539 else {540 G4ParticleDefinition * t_Definition=projectile->GetDefinition();541 projectile->SetDefinition(target->GetDefinition());542 projectile->Set4Momentum(Ptarget);543 target->SetDefinition(t_Definition);544 target->Set4Momentum(Pprojectile);545 }546 */547 //548 /*549 if(G4UniformRand() < 1.) {550 G4ParticleDefinition * t_Definition=projectile->GetDefinition();551 projectile->SetDefinition(target->GetDefinition());552 target->SetDefinition(t_Definition);553 }554 */ // For flip, for HARP555 556 G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z();557 // It is assumed that nucleon z-coordinates are ordered on increasing -----------558 559 G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e();560 561 G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z();562 if(projectile->GetSoftCollisionCount()==0) {563 projectile->SetTimeOfCreation(0.);564 target->SetTimeOfCreation(0.);565 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;566 }567 568 G4ThreeVector thePosition(projectile->GetPosition().x(),569 projectile->GetPosition().y(),570 ZcoordinateOfCurrentInteraction);571 projectile->SetPosition(thePosition);572 573 G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation();574 G4double TimeOfCurrentCollision=TimeOfPreviousCollision+575 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;576 577 projectile->SetTimeOfCreation(TimeOfCurrentCollision);578 target->SetTimeOfCreation(TimeOfCurrentCollision);579 580 projectile->Set4Momentum(Pprojectile);581 target->Set4Momentum(Ptarget);582 859 583 860 projectile->IncrementCollisionCount(1); 584 861 target->IncrementCollisionCount(1); 585 862 586 //587 //G4cout<<"Out of Excitation --------------------"<<G4endl;588 //G4int Uzhiinp; G4cin>>Uzhiinp; // Vova589 590 863 return true; 591 864 } 592 865 593 866 // --------------------------------------------------------------------- 594 G4ExcitedString * G4DiffractiveExcitation:: 595 String(G4VSplitableHadron * hadron, G4bool isProjectile) const 867 //G4ExcitedString * G4DiffractiveExcitation:: 868 // String(G4VSplitableHadron * hadron, G4bool isProjectile) const 869 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, 870 G4bool isProjectile, 871 G4ExcitedString * &FirstString, 872 G4ExcitedString * &SecondString, 873 G4FTFParameters *theParameters) const 596 874 { 597 598 //G4cout<<"G4DiffractiveExcitation::String isProj"<<isProjectile<<G4endl;599 600 875 hadron->SplitUp(); 601 876 G4Parton *start= hadron->GetNextParton(); 602 877 if ( start==NULL) 603 878 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; 604 return NULL; 879 FirstString=0; SecondString=0; 880 return; 605 881 } 606 882 G4Parton *end = hadron->GetNextParton(); 607 883 if ( end==NULL) 608 884 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; 609 return NULL; 885 FirstString=0; SecondString=0; 886 return; 610 887 } 611 612 G4ExcitedString * string; 613 if ( isProjectile ) 614 { 615 string= new G4ExcitedString(end,start, +1); 616 } else { 617 string= new G4ExcitedString(start,end, -1); 618 } 619 // Uzhi 620 //G4cout<<"G4ExcitedString * G4DiffractiveExcitation::String"<<G4endl; 621 //G4cout<<hadron->GetTimeOfCreation()<<" "<<hadron->GetPosition()/fermi<<G4endl; 622 623 string->SetTimeOfCreation(hadron->GetTimeOfCreation()); 624 string->SetPosition(hadron->GetPosition()); 888 G4LorentzVector Phadron=hadron->Get4Momentum(); 889 /* 890 G4cout<<"Create strings had "<<hadron->GetDefinition()->GetParticleName()<<" "<<Phadron<<" "<<Phadron.mag()<<G4endl; 891 G4cout<<"isProjectile "<<isProjectile<<G4endl; 892 G4cout<<"start Q "<<start->GetDefinition()->GetPDGEncoding()<<G4endl; 893 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<G4endl; 894 */ 895 G4LorentzVector Pstart(0.,0.,0.,0.); 896 G4LorentzVector Pend(0.,0.,0.,0.); 897 G4LorentzVector Pkink(0.,0.,0.,0.); 898 G4LorentzVector PkinkQ1(0.,0.,0.,0.); 899 G4LorentzVector PkinkQ2(0.,0.,0.,0.); 900 901 G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding()); 902 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding()); 903 904 //-------------------------------------------------------------------------------- 905 G4double Wmin(0.); 906 if(isProjectile) 907 { 908 Wmin=theParameters->GetProjMinDiffMass(); 909 } else 910 { 911 Wmin=theParameters->GetTarMinDiffMass(); 912 } // end of if(isProjectile) 913 914 G4double W = hadron->Get4Momentum().mag(); 915 G4double W2=W*W; 916 //G4cout<<"W Wmin "<<W<<" "<<Wmin<<G4endl; 917 G4double Pt(0.), x1(0.), x2(0.), x3(0.); 918 G4bool Kink=false; 919 920 if(W > Wmin) 921 { // Kink is possible 922 G4double Pt2kink=theParameters->GetPt2Kink(); 923 Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.)); 924 // Pt=0.; 925 //G4cout<<"Pt2kink Pt Pt2 "<<Pt2kink<<" "<<Pt<<" "<<Pt*Pt<<G4endl; 926 if(Pt > 500.*MeV) 927 { 928 G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.)); 929 G4double Y=Ymax*(1.- 2.*G4UniformRand()); 930 //G4cout<<"Ymax Y "<<Ymax<<" "<<Y<<G4endl; 931 x1=1.-Pt/W*std::exp( Y); 932 x3=1.-Pt/W*std::exp(-Y); 933 x2=2.-x1-x3; 934 //G4cout<<"X1 X2 X3 "<<x1<<" "<<x2<<" "<<x3<<G4endl; 935 G4double Mass_startQ = 650.*MeV; 936 if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV; 937 if(PDGcode_startQ == 3) Mass_startQ = 500.*MeV; 938 if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV; 939 940 G4double Mass_endQ = 650.*MeV; 941 if(PDGcode_endQ < 3) Mass_endQ = 325.*MeV; 942 if(PDGcode_endQ == 3) Mass_endQ = 500.*MeV; 943 if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV; 944 945 G4double P2_1=W2*x1*x1/4.-Mass_endQ *Mass_endQ; 946 G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ; 947 948 G4double P2_2=sqr((2.-x1-x3)*W/2.); 949 //G4cout<<"P2_1 P2_2 P2_3 "<<P2_1<<" "<<P2_2<<" "<<P2_3<<G4endl; 950 if((P2_1 < 0.) || (P2_3 <0.)) 951 { Kink=false;} 952 else 953 { 954 G4double P_1=std::sqrt(P2_1); 955 G4double P_2=std::sqrt(P2_2); 956 G4double P_3=std::sqrt(P2_3); 957 //G4cout<<"P_1 P_2 P_3 "<<P_1<<" "<<P_2<<" "<<P_3<<G4endl; 958 G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2); 959 G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3); 960 //G4cout<<"CosT12 CosT13 "<<CosT12<<" "<<CosT13<<" "<<std::acos(CosT12)*180./3.14159<<" "<<std::acos(CosT13)*180./3.14159<<G4endl; 961 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) 962 { Kink=false;} 963 else 964 { 965 Kink=true; 966 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); 967 Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1); 968 Pkink.setPx( Pt); Pkink.setPy( 0.); Pkink.setPz( P_2*CosT12); 969 Pstart.setE(x3*W/2.); 970 Pkink.setE(Pkink.vect().mag()); 971 Pend.setE(x1*W/2.); 972 973 G4double XkQ=GetQuarkFractionOfKink(0.,1.); 974 if(Pkink.getZ() > 0.) 975 { 976 if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;} 977 } else { 978 if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;} 979 } 980 981 PkinkQ2=Pkink - PkinkQ1; 982 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------ 983 /* 984 G4cout<<"Pstart "<<Pstart<<G4endl; 985 G4cout<<"Pkink "<<Pkink <<G4endl; 986 G4cout<<"Pkink1 "<<PkinkQ1<<G4endl; 987 G4cout<<"Pkink2 "<<PkinkQ2<<G4endl; 988 G4cout<<"Pend "<<Pend <<G4endl; 989 */ 990 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/ 991 std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13)); 992 G4double Psi=std::acos(Cos2Psi); 993 //G4cout<<"Psi "<<Psi<<" "<<Psi*180./twopi<<G4endl; 994 995 G4LorentzRotation Rotate; 996 if(isProjectile) {Rotate.rotateY(Psi);} 997 else {Rotate.rotateY(pi-Psi);} 998 Rotate.rotateZ(twopi*G4UniformRand()); 999 1000 //G4cout<<"Rotate"<<G4endl; 1001 1002 Pstart*=Rotate; 1003 Pkink*=Rotate; 1004 PkinkQ1*=Rotate; 1005 PkinkQ2*=Rotate; 1006 Pend*=Rotate; 1007 /* 1008 G4cout<<"Pstart "<<Pstart<<G4endl; 1009 G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl; 1010 G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl; 1011 G4cout<<"Pend "<<Pend <<G4endl; 1012 */ 1013 } 1014 } // end of if((P2_1 < 0.) || (P2_3 <0.)) 1015 } // end of if(Pt > 500.*MeV) 1016 } // end of if(W > Wmin) Check for a kink 1017 1018 //-------------------------------------------------------------------------------- 1019 1020 //G4cout<<"Kink "<<Kink<<G4endl; 1021 1022 if(Kink) 1023 { // Kink is possible 1024 std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp = 1025 theParameters->GetQuarkProbabilitiesAtGluonSplitUp(); 1026 1027 G4int QuarkInGluon(1); G4double Ksi=G4UniformRand(); 1028 for(unsigned int Iq=0; Iq <3; Iq++) 1029 { 1030 //G4cout<<Iq<<" "<<QuarkProbabilitiesAtGluonSplitUp[Iq]<<G4endl; 1031 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;} 1032 1033 //G4cout<<"Gquark "<<QuarkInGluon<<G4endl; 1034 G4Parton * Gquark = new G4Parton(QuarkInGluon); 1035 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon); 1036 /* 1037 G4cout<<Gquark->GetDefinition()->GetParticleName()<<" "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->GetDefinition()->GetPDGMass()<<G4endl; 1038 G4cout<<Ganti_quark->GetDefinition()->GetParticleName()<<" "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->GetDefinition()->GetPDGMass()<<G4endl; 1039 */ 1040 1041 //G4int Uzhi; G4cin>>Uzhi; 1042 1043 //------------------------------------------------------------------------------- 1044 /* 1045 DefineMomentumInZ(G4double aLightConeMomentum, G4bool aDirection); 1046 void Set4Momentum(const G4LorentzVector & aMomentum); 1047 G4int PDGencoding; 1048 G4ParticleDefinition * theDefinition; 1049 G4LorentzVector theMomentum; 1050 G4ThreeVector thePosition; 1051 1052 G4int theColour; 1053 G4double theIsoSpinZ; 1054 G4double theSpinZ; 1055 1056 G4double theX; 1057 */ 1058 //------------------------------------------------------------------------------- 1059 1060 //G4cout<<"Phadron "<<Phadron<<" mass "<<Phadron.mag()<<G4endl; 1061 G4LorentzRotation toCMS(-1*Phadron.boostVector()); 1062 //G4cout<<"To lab"<<G4endl; 1063 G4LorentzRotation toLab(toCMS.inverse()); 1064 1065 Pstart.transform(toLab); start->Set4Momentum(Pstart); 1066 PkinkQ1.transform(toLab); 1067 PkinkQ2.transform(toLab); 1068 Pend.transform(toLab); end->Set4Momentum(Pend); 1069 /* 1070 G4cout<<"Pstart "<<start->GetDefinition()->GetPDGEncoding()<<Pstart<<G4endl; 1071 G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl; 1072 G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl; 1073 G4cout<<"Pend "<<end->GetDefinition()->GetPDGEncoding()<<Pend <<G4endl; 1074 */ 1075 // !!! G4ExcitedString * FirstString(0); G4ExcitedString * SecondString(0); 1076 G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding()); 1077 1078 //G4cout<<"absPDGcode "<<absPDGcode<<G4endl; 1079 1080 if(absPDGcode < 1000) 1081 { // meson 1082 if ( isProjectile ) 1083 { // Projectile 1084 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end 1085 { // Quark on the end 1086 FirstString = new G4ExcitedString(end ,Ganti_quark, +1); 1087 SecondString= new G4ExcitedString(Gquark,start ,+1); 1088 Ganti_quark->Set4Momentum(PkinkQ1); 1089 Gquark->Set4Momentum(PkinkQ2); 1090 /* 1091 G4cout<<" Proj Meson end Q"<<G4endl; 1092 G4cout<<"First string ============ "<<G4endl; 1093 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1094 G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1095 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1096 G4cout<<"Secondstring ============ "<<G4endl; 1097 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1098 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1099 1100 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1101 */ 1102 } else 1103 { // Anti_Quark on the end 1104 FirstString = new G4ExcitedString(end ,Gquark, +1); 1105 SecondString= new G4ExcitedString(Ganti_quark,start ,+1); 1106 Gquark->Set4Momentum(PkinkQ1); 1107 Ganti_quark->Set4Momentum(PkinkQ2); 1108 /* 1109 G4cout<<" Proj Meson end Qbar"<<G4endl; 1110 G4cout<<"First string ============ "<<G4endl; 1111 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1112 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1113 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1114 G4cout<<"Secondstring ============ "<<G4endl; 1115 G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1116 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1117 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1118 */ 1119 } // end of if(end->GetPDGcode() > 0) 1120 } else { // Target 1121 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end 1122 { // Quark on the end 1123 FirstString = new G4ExcitedString(Ganti_quark,end ,-1); 1124 SecondString= new G4ExcitedString(start ,Gquark,-1); 1125 Ganti_quark->Set4Momentum(PkinkQ2); 1126 Gquark->Set4Momentum(PkinkQ1); 1127 /* 1128 G4cout<<" Targ Meson end Q"<<G4endl; 1129 G4cout<<"First string ============ "<<G4endl; 1130 G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1131 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1132 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1133 G4cout<<"Secondstring ============ "<<G4endl; 1134 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1135 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1136 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1137 */ 1138 } else 1139 { // Anti_Quark on the end 1140 FirstString = new G4ExcitedString(Gquark,end ,-1); 1141 SecondString= new G4ExcitedString(start ,Ganti_quark,-1); 1142 Gquark->Set4Momentum(PkinkQ2); 1143 Ganti_quark->Set4Momentum(PkinkQ1); 1144 /* 1145 G4cout<<" Targ Meson end Qbar"<<G4endl; 1146 G4cout<<"First string ============ "<<G4endl; 1147 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1148 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1149 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1150 G4cout<<"Secondstring ============ "<<G4endl; 1151 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1152 G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1153 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1154 */ 1155 } // end of if(end->GetPDGcode() > 0) 1156 } // end of if ( isProjectile ) 1157 } else // if(absPDGCode < 1000) 1158 { // Baryon/AntiBaryon 1159 if ( isProjectile ) 1160 { // Projectile 1161 if((end->GetDefinition()->GetParticleType() == "diquarks") && 1162 (end->GetDefinition()->GetPDGEncoding() > 0 ) ) 1163 { // DiQuark on the end 1164 FirstString = new G4ExcitedString(end ,Gquark, +1); 1165 SecondString= new G4ExcitedString(Ganti_quark,start ,+1); 1166 Gquark->Set4Momentum(PkinkQ1); 1167 Ganti_quark->Set4Momentum(PkinkQ2); 1168 /* 1169 G4cout<<" Proj baryon end QQ"<<G4endl; 1170 G4cout<<"First string ============ "<<G4endl; 1171 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1172 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1173 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1174 G4cout<<"Secondstring ============ "<<G4endl; 1175 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1176 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1177 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1178 */ 1179 } else 1180 { // Anti_DiQuark on the end or quark 1181 FirstString = new G4ExcitedString(end ,Ganti_quark, +1); 1182 SecondString= new G4ExcitedString(Gquark,start ,+1); 1183 Ganti_quark->Set4Momentum(PkinkQ1); 1184 Gquark->Set4Momentum(PkinkQ2); 1185 /* 1186 G4cout<<" Proj baryon end Q"<<G4endl; 1187 G4cout<<"First string ============ "<<G4endl; 1188 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1189 G4cout<<"G antQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1190 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1191 G4cout<<"Secondstring ============ "<<G4endl; 1192 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1193 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1194 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1195 */ 1196 } // end of if(end->GetPDGcode() > 0) 1197 } else { // Target 1198 1199 if((end->GetDefinition()->GetParticleType() == "diquarks") && 1200 (end->GetDefinition()->GetPDGEncoding() > 0 ) ) 1201 { // DiQuark on the end 1202 FirstString = new G4ExcitedString(Gquark,end ,-1); 1203 1204 SecondString= new G4ExcitedString(start ,Ganti_quark,-1); 1205 Gquark->Set4Momentum(PkinkQ1); 1206 Ganti_quark->Set4Momentum(PkinkQ2); 1207 /* 1208 G4cout<<" Targ baryon end QQ"<<G4endl; 1209 G4cout<<"First string ============ "<<G4endl; 1210 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1211 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1212 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1213 G4cout<<"Secondstring ============ "<<G4endl; 1214 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1215 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1216 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1217 */ 1218 } else 1219 { // Anti_DiQuark on the end or Q 1220 FirstString = new G4ExcitedString(Ganti_quark,end ,-1); 1221 SecondString= new G4ExcitedString(start ,Gquark,-1); 1222 Gquark->Set4Momentum(PkinkQ2); 1223 Ganti_quark->Set4Momentum(PkinkQ1); 1224 /* 1225 G4cout<<" Targ baryon end Q"<<G4endl; 1226 G4cout<<"First string ============ "<<G4endl; 1227 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1228 G4cout<<"end O "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1229 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1230 G4cout<<"Secondstring ============ "<<G4endl; 1231 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1232 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1233 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1234 */ 1235 } // end of if(end->GetPDGcode() > 0) 1236 } // end of if ( isProjectile ) 1237 } // end of if(absPDGcode < 1000) 1238 1239 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); 1240 FirstString->SetPosition(hadron->GetPosition()); 1241 1242 SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation()); 1243 SecondString->SetPosition(hadron->GetPosition()); 1244 1245 // ------------------------------------------------------------------------- 1246 } else // End of kink is possible 1247 { // Kink is impossible 1248 if ( isProjectile ) 1249 { 1250 FirstString= new G4ExcitedString(end,start, +1); 1251 } else { 1252 FirstString= new G4ExcitedString(start,end, -1); 1253 } 1254 1255 SecondString=0; 1256 1257 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); 1258 FirstString->SetPosition(hadron->GetPosition()); 625 1259 626 1260 // momenta of string ends 627 1261 // 628 G4double Momentum=hadron->Get4Momentum().vect().mag(); 629 G4double Plus=hadron->Get4Momentum().e() + Momentum; 630 G4double Minus=hadron->Get4Momentum().e() - Momentum; 631 632 G4ThreeVector tmp; 633 if(Momentum > 0.) 634 { 635 tmp.set(hadron->Get4Momentum().px(), 636 hadron->Get4Momentum().py(), 637 hadron->Get4Momentum().pz()); 638 tmp/=Momentum; 639 } 640 else 641 { 642 tmp.set(0.,0.,1.); 643 }; 644 645 G4LorentzVector Pstart(tmp,0.); 646 G4LorentzVector Pend(tmp,0.); 647 648 if(isProjectile) 649 { 650 Pstart*=(-1.)*Minus/2.; 651 Pend *=(+1.)*Plus /2.; 652 } 653 else 654 { 655 Pstart*=(+1.)*Plus/2.; 656 Pend *=(-1.)*Minus/2.; 657 }; 658 659 Momentum=-Pstart.mag(); 660 Pstart.setT(Momentum); // It is assumed that quark has m=0. 661 662 Momentum=-Pend.mag(); 663 Pend.setT(Momentum); // It is assumed that di-quark has m=0. 664 // 665 /* Uzhi 666 G4double ptSquared= hadron->Get4Momentum().perp2(); 667 G4double transverseMassSquared= hadron->Get4Momentum().plus() 668 * hadron->Get4Momentum().minus(); 669 670 671 G4double maxAvailMomentumSquared= 672 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) ); 673 674 G4double widthOfPtSquare = 0.25*GeV*GeV; // Uzhi 11.07 <Pt^2>=0.25 ?????????????????? 675 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared); 676 677 G4LorentzVector Pstart(G4LorentzVector(pt,0.)); 678 G4LorentzVector Pend; 679 Pend.setPx(hadron->Get4Momentum().px() - pt.x()); 680 Pend.setPy(hadron->Get4Momentum().py() - pt.y()); 681 682 G4double tm1=hadron->Get4Momentum().minus() + 683 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus(); 684 685 G4double tm2= std::sqrt( std::max(0., sqr(tm1) - 686 4. * Pend.perp2() * hadron->Get4Momentum().minus() 687 / hadron->Get4Momentum().plus() )); 688 689 G4int Sign= isProjectile ? -1 : 1; 690 691 G4double endMinus = 0.5 * (tm1 + Sign*tm2); 692 G4double startMinus= hadron->Get4Momentum().minus() - endMinus; 693 694 G4double startPlus= Pstart.perp2() / startMinus; 695 G4double endPlus = hadron->Get4Momentum().plus() - startPlus; 696 697 Pstart.setPz(0.5*(startPlus - startMinus)); 698 Pstart.setE(0.5*(startPlus + startMinus)); 699 700 Pend.setPz(0.5*(endPlus - endMinus)); 701 Pend.setE(0.5*(endPlus + endMinus)); 702 */ // Uzhi 703 start->Set4Momentum(Pstart); 704 end->Set4Momentum(Pend); 705 /* 706 G4cout<<"G4DiffractiveExcitation::String hadro"<<hadron->Get4Momentum()<<" "<<hadron->Get4Momentum().mag2()<<G4endl; 707 708 G4cout<<"G4DiffractiveExcitation::String start"<<start->Get4Momentum()<<" "<<start->GetPDGcode()<<G4endl; 709 710 G4cout<<"G4DiffractiveExcitation::String end "<< end->Get4Momentum()<<" "<< end->GetPDGcode()<<G4endl; 711 G4int Uzhi; G4cin>>Uzhi; 712 */ 1262 G4double Momentum=hadron->Get4Momentum().vect().mag(); 1263 G4double Plus=hadron->Get4Momentum().e() + Momentum; 1264 G4double Minus=hadron->Get4Momentum().e() - Momentum; 1265 1266 G4ThreeVector tmp; 1267 if(Momentum > 0.) 1268 { 1269 tmp.set(hadron->Get4Momentum().px(), 1270 hadron->Get4Momentum().py(), 1271 hadron->Get4Momentum().pz()); 1272 tmp/=Momentum; 1273 } 1274 else 1275 { 1276 tmp.set(0.,0.,1.); 1277 } 1278 1279 G4LorentzVector Pstart(tmp,0.); 1280 G4LorentzVector Pend(tmp,0.); 1281 1282 if(isProjectile) 1283 { 1284 Pstart*=(-1.)*Minus/2.; 1285 Pend *=(+1.)*Plus /2.; 1286 } 1287 else 1288 { 1289 Pstart*=(+1.)*Plus/2.; 1290 Pend *=(-1.)*Minus/2.; 1291 } 1292 1293 Momentum=-Pstart.mag(); 1294 Pstart.setT(Momentum); // It is assumed that quark has m=0. 1295 1296 Momentum=-Pend.mag(); 1297 Pend.setT(Momentum); // It is assumed that di-quark has m=0. 1298 1299 start->Set4Momentum(Pstart); 1300 end->Set4Momentum(Pend); 1301 SecondString=0; 1302 } // End of kink is impossible 1303 713 1304 #ifdef G4_FTFDEBUG 714 G4cout << " generated string flavors "715 << start->GetPDGcode() << " / "716 << end->GetPDGcode() << G4endl;717 G4cout << " generated string momenta: quark "718 << start->Get4Momentum() << "mass : "719 <<start->Get4Momentum().mag() << G4endl;720 G4cout << " generated string momenta: Diquark "721 << end ->Get4Momentum()722 << "mass : " <<end->Get4Momentum().mag()<< G4endl;723 G4cout << " sum of ends " << Pstart+Pend << G4endl;724 G4cout << " Original " << hadron->Get4Momentum() << G4endl;1305 G4cout << " generated string flavors " 1306 << start->GetPDGcode() << " / " 1307 << end->GetPDGcode() << G4endl; 1308 G4cout << " generated string momenta: quark " 1309 << start->Get4Momentum() << "mass : " 1310 <<start->Get4Momentum().mag() << G4endl; 1311 G4cout << " generated string momenta: Diquark " 1312 << end ->Get4Momentum() 1313 << "mass : " <<end->Get4Momentum().mag()<< G4endl; 1314 G4cout << " sum of ends " << Pstart+Pend << G4endl; 1315 G4cout << " Original " << hadron->Get4Momentum() << G4endl; 725 1316 #endif 726 1317 727 return string; 1318 return; 1319 728 1320 } 729 1321 … … 732 1324 733 1325 // --------------------------------------------------------------------- 734 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi1326 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const 735 1327 { 736 1328 // choose an x between Xmin and Xmax with P(x) ~ 1/x 737 1329 // to be improved... 738 1330 739 G4double range=Pmax-Pmin; // Uzhi1331 G4double range=Pmax-Pmin; 740 1332 741 1333 if ( Pmin <= 0. || range <=0. ) … … 745 1337 } 746 1338 747 G4double P; 748 /* // Uzhi 749 do { 750 x=Xmin + G4UniformRand() * range; 751 } while ( Xmin/x < G4UniformRand() ); 752 */ // Uzhi 753 754 P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); // Uzhi 755 756 //debug-hpw cout << "DiffractiveX "<<x<<G4endl; 1339 G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); 757 1340 return P; 758 1341 } … … 760 1343 // --------------------------------------------------------------------- 761 1344 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, 762 G4double maxPtSquare) const // Uzhi1345 G4double maxPtSquare) const 763 1346 { // @@ this method is used in FTFModel as well. Should go somewhere common! 764 1347 765 1348 G4double Pt2; 766 /* // Uzhi767 do {768 pt2=widthSquare * std::log( G4UniformRand() );769 } while ( pt2 > maxPtSquare);770 */ // Uzhi771 772 1349 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 773 (std::exp(-maxPtSquare/AveragePt2)-1.)); // Uzhi1350 (std::exp(-maxPtSquare/AveragePt2)-1.)); 774 1351 775 1352 G4double Pt=std::sqrt(Pt2); 776 777 1353 G4double phi=G4UniformRand() * twopi; 778 779 1354 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 780 1355 } 781 1356 1357 // --------------------------------------------------------------------- 1358 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const 1359 { 1360 G4double z, yf; 1361 do { 1362 z = zmin + G4UniformRand()*(zmax-zmin); 1363 yf = z*z +sqr(1 - z); 1364 } 1365 while (G4UniformRand() > yf); 1366 return z; 1367 } 1368 // --------------------------------------------------------------------- 1369 void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09 1370 { 1371 G4int absIdPDG = std::abs(IdPDG); 1372 Q1= absIdPDG/ 100; 1373 Q2= (absIdPDG %100)/10; 1374 1375 G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 ); 1376 1377 if (IdPDG < 0 ) anti *=-1; 1378 Q1 *= anti; 1379 Q2 *= -1 * anti; 1380 return; 1381 } 1382 // --------------------------------------------------------------------- 1383 void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG, 1384 G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09 1385 { 1386 Q1 = IdPDG / 1000; 1387 Q2 = (IdPDG % 1000) / 100; 1388 Q3 = (IdPDG % 100) / 10; 1389 return; 1390 } 1391 // --------------------------------------------------------------------- 1392 G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09 1393 { 1394 G4int TmpQ(0); 1395 if( Q3 > Q2 ) 1396 { 1397 TmpQ = Q2; 1398 Q2 = Q3; 1399 Q3 = TmpQ; 1400 } else if( Q3 > Q1 ) 1401 { 1402 TmpQ = Q1; 1403 Q1 = Q3; 1404 Q3 = TmpQ; 1405 } 1406 1407 if( Q2 > Q1 ) 1408 { 1409 TmpQ = Q1; 1410 Q1 = Q2; 1411 Q2 = TmpQ; 1412 } 1413 1414 G4int NewCode = Q1*1000 + Q2* 100 + Q3* 10 + 2; 1415 return NewCode; 1416 } 782 1417 // --------------------------------------------------------------------- 783 1418 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &) -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveHHScatterer.cc
r962 r1196 38 38 {} 39 39 40 // ------------------------------------------------------------------- 40 void G4DiffractiveHHScatterer::CreateStrings() 41 /* 42 G4VSplitableHadron * aHadron, 43 G4bool isProjectile, 44 G4ExcitedString * FirstString, 45 G4ExcitedString * SecondString, 46 G4FTFParameters *theParameters) 47 */ 48 const 49 50 {} 51 52 /* ------------------------------------------------------------------- 41 53 G4KineticTrackVector * G4DiffractiveHHScatterer:: 42 54 Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack) … … 76 88 return result; 77 89 } 90 */ -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveSplitableHadron.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4DiffractiveSplitableHadron.cc,v 1. 7 2008/03/31 15:34:01vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4DiffractiveSplitableHadron.cc,v 1.8 2009/07/31 11:03:00 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 87 87 88 88 G4int PDGcode=GetDefinition()->GetPDGEncoding(); 89 89 90 G4int stringStart, stringEnd; 90 91 ChooseStringEnds(PDGcode, &stringStart,&stringEnd); 91 92 92 93 Parton[0] = new G4Parton(stringStart); 93 94 Parton[1] = new G4Parton(stringEnd); 94 95 PartonIndex=-1; 95 96 96 } 97 97 -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4ElasticHNScattering.cc
r968 r1196 25 25 // 26 26 // 27 // $Id: G4ElasticHNScattering.cc,v 1. 3 2008/05/19 12:56:36 vuzhinsk Exp $27 // $Id: G4ElasticHNScattering.cc,v 1.11 2009/10/06 10:10:36 vuzhinsk Exp $ 28 28 // ------------------------------------------------------------ 29 29 // GEANT 4 class implemetation file … … 46 46 #include "G4VSplitableHadron.hh" 47 47 #include "G4ExcitedString.hh" 48 #include "G4FTFParameters.hh" // Uzhi 29.03.0848 #include "G4FTFParameters.hh" 49 49 //#include "G4ios.hh" 50 50 … … 59 59 { 60 60 //G4cout<<"G4ElasticHNScattering::ElasticScattering"<<G4endl; 61 61 // -------------------- Projectile parameters ----------------------------------- 62 62 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 63 64 // -------------------- Projectile parameters ----------------------------------- 65 G4bool PutOnMassShell=0; 63 //G4cout<<"Elastic P "<<Pprojectile<<G4endl; 64 if(Pprojectile.z() < 0.) 65 { 66 target->SetStatus(2); 67 return false; 68 } 69 70 // G4double ProjectileRapidity = Pprojectile.rapidity(); 71 // G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding(); 72 // G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition(); 73 74 G4bool PutOnMassShell(false); 66 75 67 76 G4double M0projectile = Pprojectile.mag(); 68 69 77 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 70 71 PutOnMassShell= 1;78 { 79 PutOnMassShell=true; 72 80 M0projectile=projectile->GetDefinition()->GetPDGMass(); 73 81 } 74 82 75 83 G4double Mprojectile2 = M0projectile * M0projectile; 76 84 77 // G4double AveragePt2=theParameters->GetSlope(); // Uzhi ???78 // AveragePt2 = AveragePt2 * GeV*GeV;79 80 85 G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering(); 81 86 82 87 // -------------------- Target parameters ---------------------------------------------- 88 // G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding(); 89 83 90 G4LorentzVector Ptarget=target->Get4Momentum(); 84 91 //G4cout<<"Elastic T "<<Ptarget<<G4endl; 92 // G4double TargetRapidity = Ptarget.rapidity(); 85 93 G4double M0target = Ptarget.mag(); 86 //G4cout<<" Mp Mt Pt2 "<<M0projectile<<" "<<M0target<<" "<<AveragePt2/GeV/GeV<<G4endl;87 94 88 95 if(M0target < target->GetDefinition()->GetPDGMass()) 89 90 PutOnMassShell= 1;96 { 97 PutOnMassShell=true; 91 98 M0target=target->GetDefinition()->GetPDGMass(); 92 99 } 93 100 94 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2();95 // for AA-inter. 101 G4double Mtarget2 = M0target * M0target; 102 96 103 // Transform momenta to cms and then rotate parallel to z axis; 97 104 … … 103 110 G4LorentzVector Ptmp=toCms*Pprojectile; 104 111 105 if ( Ptmp.pz() <= 0. ) // Uzhi ???112 if ( Ptmp.pz() <= 0. ) 106 113 { 107 114 // "String" moving backwards in CMS, abort collision !! 108 115 //G4cout << " abort Collision!! " << G4endl; 116 target->SetStatus(2); 109 117 return false; 110 118 } … … 118 126 Ptarget.transform(toCms); 119 127 120 // ---------------------- Sampling of transfered Pt ------------------------ 128 // ---------------------- Putting on mass-on-shell, if needed ------------------------ 129 G4double PZcms2, PZcms; 130 131 G4double S=Psum.mag2(); 132 // G4double SqrtS=std::sqrt(S); 133 134 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2- 135 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S; 136 137 if(PZcms2 < 0.) 138 { // It can be in an interaction with off-shell nuclear nucleon 139 if(M0projectile > projectile->GetDefinition()->GetPDGMass()) 140 { // An attempt to de-excite the projectile 141 // It is assumed that the target is in the ground state 142 M0projectile = projectile->GetDefinition()->GetPDGMass(); 143 Mprojectile2=M0projectile*M0projectile; 144 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2- 145 2*S*Mprojectile2 - 2*S*Mtarget2 - 2*Mprojectile2*Mtarget2) 146 /4./S; 147 if(PZcms2 < 0.){ return false;} // Non succesful attempt after the de-excitation 148 } 149 else // if(M0projectile > projectile->GetDefinition()->GetPDGMass()) 150 { 151 target->SetStatus(2); // Uzhi 17.07.09 152 return false; // The projectile was not excited, 153 // but the energy was too low to put 154 // the target nucleon on mass-shell 155 } // end of if(M0projectile > projectile->GetDefinition()->GetPDGMass()) 156 } // end of if(PZcms2 < 0.) 157 158 PZcms = std::sqrt(PZcms2); 159 160 if(PutOnMassShell) 161 { 162 if(Pprojectile.z() > 0.) 163 { 164 Pprojectile.setPz( PZcms); 165 Ptarget.setPz( -PZcms); 166 } 167 else // if(Pprojectile.z() > 0.) 168 { 169 Pprojectile.setPz(-PZcms); 170 Ptarget.setPz( PZcms); 171 }; 172 173 Pprojectile.setE(std::sqrt(Mprojectile2+ 174 Pprojectile.x()*Pprojectile.x()+ 175 Pprojectile.y()*Pprojectile.y()+ 176 PZcms2)); 177 Ptarget.setE(std::sqrt( Mtarget2 + 178 Ptarget.x()*Ptarget.x()+ 179 Ptarget.y()*Ptarget.y()+ 180 PZcms2)); 181 } // end of if(PutOnMassShell) 182 183 G4double maxPtSquare = PZcms2; 184 185 //----- Charge exchange between the projectile and the target, if possible 186 // (G4UniformRand() < 0.5*std::exp(-0.5*(ProjectileRapidity - TargetRapidity)))) 187 /* 188 if((ProjectilePDGcode != TargetPDGcode) && 189 ((ProjectilePDGcode > 1000) && (TargetPDGcode > 1000)) && 190 (G4UniformRand() < 1.0*std::exp(-0.5*(ProjectileRapidity - TargetRapidity)))) 191 { 192 projectile->SetDefinition(target->GetDefinition()); 193 target->SetDefinition(ProjectileDefinition); 194 } 195 */ 196 // ------ Now we can calculate the transfered Pt -------------------------- 121 197 G4double Pt2; 122 198 G4double ProjMassT2, ProjMassT; 123 G4double TargMassT2, TargMassT; 124 G4double PZcms2, PZcms; 125 126 G4double S=Psum.mag2(); 127 // G4double SqrtS=std::sqrt(S); 128 129 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2- 130 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S; 131 if(PZcms2 < 0) 132 {return false;} // It can be in an interaction with off-shell nuclear nucleon 133 134 PZcms = std::sqrt(PZcms2); 135 136 if(PutOnMassShell) 137 { 138 if(Pprojectile.z() > 0.) 139 { 140 Pprojectile.setPz( PZcms); 141 Ptarget.setPz( -PZcms); 142 } 143 else 144 { 145 Pprojectile.setPz(-PZcms); 146 Ptarget.setPz( PZcms); 147 }; 148 149 Pprojectile.setE(std::sqrt(Mprojectile2+ 150 Pprojectile.x()*Pprojectile.x()+ 151 Pprojectile.y()*Pprojectile.y()+ 152 PZcms2)); 153 Ptarget.setE(std::sqrt( Mtarget2 + 154 Ptarget.x()*Ptarget.x()+ 155 Ptarget.y()*Ptarget.y()+ 156 PZcms2)); 157 } 158 159 G4double maxPtSquare = PZcms2; 199 G4double TargMassT2, TargMassT; 200 160 201 G4LorentzVector Qmomentum; 161 202 … … 163 204 164 205 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 165 //G4cout<<"Pt2 GeV^2 "<<(Pt2)/GeV/GeV<<G4endl;166 206 167 207 ProjMassT2=Mprojectile2+Pt2; … … 175 215 2.*S*ProjMassT2-2.*S*TargMassT2- 176 216 2.*ProjMassT2*TargMassT2)/4./S; 177 if(PZcms2 < 0 ) {PZcms2=0;}; 217 if(PZcms2 < 0 ) {PZcms2=0;};// to avoid the exactness problem 178 218 PZcms =std::sqrt(PZcms2); 179 219 180 Pprojectile.setPz( PZcms); // Uzhi Proj can move backward 181 Ptarget.setPz( -PZcms); // Uzhi Proj can move backward 182 183 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl; 184 // G4cout << "pt2" << pt2 << G4endl; 185 // G4cout << "Qmomentum " << Qmomentum << G4endl; 186 // G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() << 187 // " / " << (Ptarget-Qmomentum).mag() << G4endl; 220 Pprojectile.setPz( PZcms); 221 Ptarget.setPz( -PZcms); 188 222 189 223 Pprojectile += Qmomentum; 190 224 Ptarget -= Qmomentum; 191 192 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;193 //G4cout << "Ptarget with Q : " << Ptarget << G4endl;194 195 // G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;196 // G4cout << "Target back: " << toLab * Ptarget << G4endl;197 225 198 226 // Transform back and update SplitableHadron Participant. 199 227 Pprojectile.transform(toLab); 200 228 Ptarget.transform(toLab); 201 202 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl; 203 //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl; 204 205 //G4cout << "Target mass " << Ptarget.mag() << G4endl; 206 207 //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl; 208 209 G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z(); 210 // It is assumed that nucleon z-coordinates are ordered on increasing ----------- 211 212 G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e(); 213 214 G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z(); 215 if(projectile->GetSoftCollisionCount()==0) { 216 projectile->SetTimeOfCreation(0.); 217 target->SetTimeOfCreation(0.); 218 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction; 219 } 220 221 G4ThreeVector thePosition(projectile->GetPosition().x(), 222 projectile->GetPosition().y(), 223 ZcoordinateOfCurrentInteraction); 224 projectile->SetPosition(thePosition); 225 226 G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation(); 227 G4double TimeOfCurrentCollision=TimeOfPreviousCollision+ 228 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z; 229 230 projectile->SetTimeOfCreation(TimeOfCurrentCollision); 231 target->SetTimeOfCreation(TimeOfCurrentCollision); 229 /* // Maybe it will be needed for an exact calculations-------------------- 230 G4double TargetMomentum=std::sqrt(Ptarget.x()*Ptarget.x()+ 231 Ptarget.y()*Ptarget.y()+ 232 Ptarget.z()*Ptarget.z()); 233 */ 234 235 // Calculation of the creation time --------------------- 236 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); 237 projectile->SetPosition(target->GetPosition()); 238 // Creation time and position of target nucleon were determined at 239 // ReggeonCascade() of G4FTFModel 240 // ------------------------------------------------------ 232 241 233 242 projectile->Set4Momentum(Pprojectile); -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4FTFModel.cc,v 1. 13 2008/12/09 10:40:52vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4FTFModel.cc,v 1.28 2009/10/29 14:55:33 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 38 38 39 39 #include "G4FTFModel.hh" 40 #include "G4FTFParameters.hh" // Uzhi 29.03.0840 #include "G4FTFParameters.hh" 41 41 #include "G4FTFParticipants.hh" 42 #include "G4DiffractiveSplitableHadron.hh" 42 43 #include "G4InteractionContent.hh" 43 44 #include "G4LorentzRotation.hh" 44 45 #include "G4ParticleDefinition.hh" 46 #include "G4ParticleTable.hh" 45 47 #include "G4ios.hh" 46 #include <utility> // Uzhi 29.03.08 48 #include <utility> 49 #include "G4IonTable.hh" 47 50 48 51 // Class G4FTFModel 49 52 50 53 G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()), 51 theElastic(new G4ElasticHNScattering()) // Uzhi 29.03.0854 theElastic(new G4ElasticHNScattering()) 52 55 { 53 56 G4VPartonStringModel::SetThisPointer(this); 54 theParameters=0; // Uzhi 9.12.08 55 } 56 57 /* 58 G4FTFModel::G4FTFModel(G4double , G4double , G4double ):theExcitation(new // Uzhi 9.12.08 G4DiffractiveExcitation()) 59 { 60 G4VPartonStringModel::SetThisPointer(this); 61 } 62 63 G4FTFModel::G4FTFModel(G4DiffractiveExcitation * anExcitation) 64 : 65 theExcitation(anExcitation) 66 { 67 G4VPartonStringModel::SetThisPointer(this); 68 } 69 */ 57 theParameters=0; 58 } 70 59 71 60 72 61 G4FTFModel::~G4FTFModel() 73 62 { 74 if( theParameters != 0 ) delete theParameters; // Uzhi 5.12.0863 if( theParameters != 0 ) delete theParameters; 75 64 // Because FTF model can be called for various particles 76 65 // theParameters must be erased at the end of each call. 77 // Thus the delete is olso in G4FTFModel::GetStrings() method 78 if( theExcitation != 0 ) delete theExcitation; // Uzhi 5.12.08 79 if( theElastic != 0 ) delete theElastic; // Uzhi 5.12.08 80 } 81 66 // Thus the delete is also in G4FTFModel::GetStrings() method 67 if( theExcitation != 0 ) delete theExcitation; 68 if( theElastic != 0 ) delete theElastic; 69 } 82 70 83 71 const G4FTFModel & G4FTFModel::operator=(const G4FTFModel &) … … 87 75 } 88 76 89 90 77 int G4FTFModel::operator==(const G4FTFModel &right) const 91 78 { … … 102 89 { 103 90 theProjectile = aProjectile; 104 //G4cout<<"G4FTFModel::Init "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;105 91 theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); 106 // Uzhi N-mass number Z-charge ------------------------- Uzhi 29.03.0892 // ----------- N-mass number Z-charge ------------------------- 107 93 108 94 // --- cms energy … … 112 98 2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass(); 113 99 /* 100 G4cout<<"----------------------------------------"<<G4endl; 114 101 G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl; 115 102 G4cout << " primary Mass (GeV): " << theProjectile.GetMass() /GeV << G4endl; 103 G4cout << " primary 3Mom " << theProjectile.GetMomentum() << G4endl; 104 G4cout << " primary space position " << theProjectile.GetPositionInNucleus() << G4endl; 116 105 G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl; 117 106 */ 118 107 119 if( theParameters != 0 ) delete theParameters; // Uzhi 9.12.08108 if( theParameters != 0 ) delete theParameters; 120 109 theParameters = new G4FTFParameters(theProjectile.GetDefinition(), 121 110 aNucleus.GetN(),aNucleus.GetZ(), 122 s);// ------------------------- Uzhi 19.04.08 123 //theParameters->SetProbabilityOfElasticScatt(0.); // To turn on/off (1/0) elastic scattering 111 s); 112 //theParameters->SetProbabilityOfElasticScatt(0.); 113 // To turn on/off (1/0) elastic scattering 114 124 115 } 125 116 … … 127 118 G4ExcitedStringVector * G4FTFModel::GetStrings() 128 119 { 129 //G4cout<<"theParticipants.GetList"<<G4endl; 120 //G4int Uzhi; G4cin>>Uzhi; 121 122 //G4cout<<"GetList"<<G4endl; 130 123 theParticipants.GetList(theProjectile,theParameters); 131 //G4cout<<"ExciteParticipants()"<<G4endl; 132 if (! ExciteParticipants()) return NULL;; 133 //G4cout<<"theStrings = BuildStrings()"<<G4endl; 124 //G4cout<<"Regge"<<G4endl; 125 ReggeonCascade(); // Uzhi 26 July 09 126 //G4cout<<"On mass shell"<<G4endl; 127 if (! PutOnMassShell() ) return NULL; // Uzhi 26 July 09 128 //G4cout<<"Excite"<<G4endl; 129 if (! ExciteParticipants()) return NULL; 130 //G4cout<<"Strings"<<G4endl; 134 131 G4ExcitedStringVector * theStrings = BuildStrings(); 135 //G4cout<<"Return to theStrings "<<G4endl; 136 if( theParameters != 0 ) // Uzhi 9.12.08 137 { // Uzhi 9.12.08 138 delete theParameters; // Uzhi 9.12.08 139 theParameters=0; // Uzhi 9.12.08 140 } // Uzhi 9.12.08 132 //G4cout<<"Out FTF N strings "<<theStrings->size()<<G4endl; 133 //G4cout<<"GetResidualNucleus"<<G4endl; 134 GetResidualNucleus(); 135 //G4cout<<"Out of FTF"<<G4endl; 136 if( theParameters != 0 ) 137 { 138 delete theParameters; 139 theParameters=0; 140 } 141 141 return theStrings; 142 142 } 143 143 144 144 // ------------------------------------------------------------ 145 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){delete aH;} }; 145 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} }; 146 147 // ------------------------------------------------------------ 148 void G4FTFModel::ReggeonCascade() // Uzhi 26 July 2009 149 { //--- Implementation of reggeon theory inspired model------- 150 NumberOfInvolvedNucleon=0; 151 152 //G4int PrimInt(0); 153 theParticipants.StartLoop(); 154 while (theParticipants.Next()) 155 { 156 //PrimInt++; 157 const G4InteractionContent & collision=theParticipants.GetInteraction(); 158 G4Nucleon * TargetNucleon=collision.GetTargetNucleon(); 159 //G4cout<<"Prim Nnucl "<<TargetNucleon->Get4Momentum()<<G4endl; 160 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon; 161 NumberOfInvolvedNucleon++; 162 163 G4double XofWoundedNucleon = TargetNucleon->GetPosition().x(); 164 G4double YofWoundedNucleon = TargetNucleon->GetPosition().y(); 165 166 theParticipants.theNucleus->StartLoop(); 167 G4Nucleon * Neighbour; 168 //G4int NucleonNum(0); 169 while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) ) 170 { 171 if(!Neighbour->AreYouHit()) 172 { 173 G4double impact2= sqr(XofWoundedNucleon - Neighbour->GetPosition().x()) + 174 sqr(YofWoundedNucleon - Neighbour->GetPosition().y()); 175 176 if(G4UniformRand() < theParameters->GetCofNuclearDestruction()* 177 std::exp(-impact2/theParameters->GetR2ofNuclearDestruction())) 178 { // The neighbour nucleon is involved in the reggeon cascade 179 //G4cout<<" involved "<<NucleonNum<<" "<<Neighbour->Get4Momentum()<<G4endl; 180 TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour; 181 NumberOfInvolvedNucleon++; 182 //G4cout<<" PrimInt"<<" "<<NumberOfInvolvedNucleon<<G4endl; 183 184 G4VSplitableHadron *targetSplitable; 185 targetSplitable = new G4DiffractiveSplitableHadron(*Neighbour); 186 187 Neighbour->Hit(targetSplitable); 188 // Neighbour->SetBindingEnergy(3.*Neighbour->GetBindingEnergy()); // Uzhi 5.10.09 189 targetSplitable->SetStatus(2); 190 } 191 } // end of if(!Neighbour->AreYouHit()) 192 //NucleonNum++; 193 } // end of while (theParticipant.theNucleus->GetNextNucleon()) 194 //G4cout<<"Prim Int N Ninvolv "<<PrimInt<<" "<<NumberOfInvolvedNucleon<<G4endl; 195 } // end of while (theParticipants.Next()) 196 //G4cout<<"At end "<<PrimInt<<" "<<NumberOfInvolvedNucleon<<G4endl; 197 198 // ---------------- Calculation of creation time for each target nucleon ----------- 199 theParticipants.StartLoop(); // restart a loop 200 theParticipants.Next(); 201 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile(); 202 G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e(); 203 primary->SetTimeOfCreation(0.); 204 205 G4double ZcoordinateOfPreviousCollision(0.); 206 G4double ZcoordinateOfCurrentInteraction(0.); 207 G4double TimeOfPreviousCollision(0.); 208 G4double TimeOfCurrentCollision(0); 209 210 theParticipants.theNucleus->StartLoop(); 211 G4Nucleon * aNucleon; 212 G4bool theFirstInvolvedNucleon(true); 213 while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) ) 214 { 215 if(aNucleon->AreYouHit()) 216 { 217 if(theFirstInvolvedNucleon) 218 { 219 ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z(); 220 theFirstInvolvedNucleon=false; 221 } 222 223 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z(); 224 TimeOfCurrentCollision=TimeOfPreviousCollision+ 225 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z; 226 // It is assumed that the nucleons are ordered on increasing z-coordinate ------------ 227 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision); 228 229 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction; 230 TimeOfPreviousCollision=TimeOfCurrentCollision; 231 } // end of if(aNucleon->AreYouHit()) 232 } // end of while (theParticipant.theNucleus->GetNextNucleon()) 233 // 234 // The algorithm can be improved, but it will be more complicated, and will require 235 // changes in G4DiffractiveExcitation.cc and G4ElasticHNScattering.cc 236 } // Uzhi 26 July 2009 237 238 // ------------------------------------------------------------ 239 G4bool G4FTFModel::PutOnMassShell() 240 { 241 // -------------- Properties of the projectile ---------------- 242 //G4cout<<"Put on Mass-shell"<<G4endl; 243 theParticipants.StartLoop(); // restart a loop 244 theParticipants.Next(); 245 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile(); 246 G4LorentzVector Pprojectile=primary->Get4Momentum(); 247 //G4cout<<"Proj "<<Pprojectile<<G4endl; 248 if(Pprojectile.z() < 0.){return false;} 249 250 G4double Mprojectile = Pprojectile.mag(); 251 G4double M2projectile = Pprojectile.mag2(); 252 //------------------------------------------------------------- 253 G4LorentzVector Psum = Pprojectile; 254 G4double SumMasses = Mprojectile; 255 256 //--------------- Target nucleus ------------------------------ 257 G4V3DNucleus *theNucleus = GetWoundedNucleus(); 258 G4Nucleon * aNucleon; 259 G4int ResidualMassNumber=theNucleus->GetMassNumber(); 260 G4int ResidualCharge =theNucleus->GetCharge(); 261 //G4cout<<"Nucleus "<<ResidualMassNumber<<" "<<ResidualCharge<<G4endl; 262 ResidualExcitationEnergy=0.; 263 G4LorentzVector PnuclearResidual(0.,0.,0.,0.); 264 265 G4double ExcitationEnergyPerWoundedNucleon= 266 theParameters->GetExcitationEnergyPerWoundedNucleon(); 267 268 theNucleus->StartLoop(); 269 //G4int NucleonNum(0); 270 while ((aNucleon = theNucleus->GetNextNucleon())) 271 { 272 if(aNucleon->AreYouHit()) 273 { // Involved nucleons 274 //G4cout<<" "<<NucleonNum<<" "<<aNucleon->Get4Momentum()<<G4endl; 275 Psum += aNucleon->Get4Momentum(); 276 SumMasses += aNucleon->GetDefinition()->GetPDGMass(); 277 ResidualMassNumber--; 278 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge(); 279 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon; 280 } 281 else 282 { // Spectator nucleons 283 PnuclearResidual += aNucleon->Get4Momentum(); 284 } // end of if(!aNucleon->AreYouHit()) 285 //NucleonNum++; 286 } // end of while (theNucleus->GetNextNucleon()) 287 288 //G4cout<<"Nucleus "<<ResidualMassNumber<<" "<<ResidualCharge<<G4endl; 289 //G4cout<<"PResid "<<PnuclearResidual<<G4endl; 290 291 Psum += PnuclearResidual; 292 G4double ResidualMass(0.); 293 if(ResidualMassNumber == 0) 294 { 295 ResidualMass=0.; 296 ResidualExcitationEnergy=0.; 297 } 298 else 299 { 300 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()-> 301 GetIonMass(ResidualCharge ,ResidualMassNumber); 302 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;} 303 } 304 305 // ResidualMass +=ResidualExcitationEnergy; // Will be given after checks 306 SumMasses += ResidualMass; 307 308 //------------------------------------------------------------- 309 310 G4double SqrtS=Psum.mag(); 311 G4double S=Psum.mag2(); 312 //G4cout<<" SqrtS SumMasses "<<SqrtS <<" "<< SumMasses<<G4endl; 313 //G4cout<<"Res M E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl; 314 if(SqrtS < SumMasses) {return false;} // It is impossible to simulate 315 // after putting nuclear nucleons 316 // on mass-shell 317 if(SqrtS < SumMasses+ResidualExcitationEnergy) {ResidualExcitationEnergy=0.;} 318 319 ResidualMass +=ResidualExcitationEnergy; 320 SumMasses +=ResidualExcitationEnergy; 321 //G4cout<<"Res M E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl; 322 323 //------------------------------------------------------------- 324 // Sampling of nucleons what are transfered to delta-isobars -- 325 G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV)); 326 G4int NumberOfDeltas(0); 327 //SumMasses=Mprojectile; 328 if(theNucleus->GetMassNumber() != 1) 329 { 330 G4double ProbDeltaIsobar(0.); // 1. ******************************* 331 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 332 { 333 if((G4UniformRand() < ProbDeltaIsobar)&&(NumberOfDeltas < MaxNumberOfDeltas)) 334 { 335 NumberOfDeltas++; 336 G4VSplitableHadron * targetSplitable=TheInvolvedNucleon[i]->GetSplitableHadron(); 337 SumMasses-=targetSplitable->GetDefinition()->GetPDGMass(); 338 339 G4int PDGcode = targetSplitable->GetDefinition()->GetPDGEncoding(); 340 G4int newPDGcode = PDGcode/10; newPDGcode=newPDGcode*10+4; // Delta 341 G4ParticleDefinition* ptr = 342 G4ParticleTable::GetParticleTable()->FindParticle(newPDGcode); 343 targetSplitable->SetDefinition(ptr); 344 SumMasses+=targetSplitable->GetDefinition()->GetPDGMass(); 345 //G4cout<<" i "<<i<<" Delta "<<targetSplitable->GetDefinition()->GetPDGMass()<<G4endl; 346 } 347 else 348 { 349 // SumMasses+=TheInvolvedNucleon[i]->GetSplitableHadron()-> 350 // GetDefinition()->GetPDGMass(); 351 //G4cout<<" i "<<i<<" Nuclo "<<TheInvolvedNucleon[i]->GetSplitableHadron()->GetDefinition()->GetPDGMass()<<G4endl; 352 } 353 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 354 } // end of if(theNucleus.GetMassNumber() != 1) 355 //G4cout<<"MaxNumberOfDeltas NumberOfDeltas "<<MaxNumberOfDeltas<<" "<<NumberOfDeltas<<G4endl; 356 //G4cout<<" SqrtS SumMasses "<<SqrtS <<" "<< SumMasses<<G4endl; 357 //------------------------------------------------------------- 358 359 G4LorentzRotation toCms(-1*Psum.boostVector()); 360 G4LorentzVector Ptmp=toCms*Pprojectile; 361 362 if ( Ptmp.pz() <= 0. ) 363 { // "String" moving backwards in CMS, abort collision !! 364 //G4cout << " abort ColliDeleteVSplitableHadronsion!! " << G4endl; 365 return false; 366 } 367 368 toCms.rotateZ(-1*Ptmp.phi()); 369 toCms.rotateY(-1*Ptmp.theta()); 370 371 G4LorentzRotation toLab(toCms.inverse()); 372 373 // Pprojectile.transform(toCms); 374 //G4cout<<"Proj in CMS "<<Pprojectile<<G4endl; 375 376 //G4cout<<"Main work"<<G4endl; 377 //------------------------------------------------------------- 378 //------- Ascribing of the involved nucleons Pt and Xminus ---- 379 G4double Dcor = theParameters->GetDofNuclearDestruction()/ 380 theNucleus->GetMassNumber(); 381 // NumberOfInvolvedNucleon; 382 G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction(); 383 G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction(); 384 //G4cout<<" D Pt2 "<<Dcor<<" "<<AveragePt2<<G4endl; 385 386 G4double M2target(0.); 387 G4int NumberOfTries(0); 388 G4double ScaleFactor(1.); 389 do // while (SqrtS < Mprojectile + std::sqrt(M2target)) 390 { // while (DecayMomentum < 0.) 391 392 NumberOfTries++; 393 //G4cout<<"NumberOfTries "<<NumberOfTries<<G4endl; 394 if(NumberOfTries == 100*(NumberOfTries/100)) 395 { // At large number of tries it would be better to reduce the values 396 ScaleFactor/=2.; 397 Dcor *=ScaleFactor; 398 AveragePt2 *=ScaleFactor; 399 //G4cout<<"NumberOfTries "<<NumberOfTries<<" "<<Dcor<<" "<<AveragePt2<<G4endl; 400 //G4int Uzhi; G4cin>>Uzhi; 401 } 402 //G4cout<<"Start Decay "<<G4endl; G4int Uzhi; G4cin>>Uzhi; 403 G4ThreeVector PtSum(0.,0.,0.); 404 G4double XminusSum(0.); 405 G4double Xminus(0.); 406 G4bool Success=true; 407 408 do // while(Success == false); 409 { 410 //G4cout<<"Sample Pt and X"<<" Ninv "<<NumberOfInvolvedNucleon<<G4endl; // G4int Uzhi1; G4cin>>Uzhi1; 411 Success=true; 412 413 PtSum =G4ThreeVector(0.,0.,0.); 414 XminusSum=0.; 415 416 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 417 { 418 G4Nucleon * aNucleon = TheInvolvedNucleon[i]; 419 420 G4ThreeVector tmpPt = GaussianPt(AveragePt2, maxPtSquare); 421 PtSum += tmpPt; 422 423 G4ThreeVector tmpX=GaussianPt(Dcor*Dcor, 1.); 424 Xminus=tmpX.x(); 425 XminusSum+=Xminus; 426 427 G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,0.); 428 aNucleon->SetMomentum(tmp); 429 //G4cout<<"Nucl mom "<<aNucleon->GetMomentum()<<G4endl; 430 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 431 432 //--------------------------------------------------------------------------- 433 G4double DeltaX(0.); 434 G4double DeltaY(0.); 435 G4double DeltaXminus(0.); 436 437 if(ResidualMassNumber == 0) 438 { 439 DeltaX = PtSum.x()/NumberOfInvolvedNucleon; 440 DeltaY = PtSum.y()/NumberOfInvolvedNucleon; 441 DeltaXminus = (XminusSum-1.)/NumberOfInvolvedNucleon; 442 } 443 else 444 { 445 DeltaXminus = -1./theNucleus->GetMassNumber(); 446 } 447 //G4cout<<"Deltas "<<DeltaX<<" "<<DeltaY<<" "<<DeltaXminus<<G4endl; 448 449 XminusSum=1.; 450 M2target =0.; 451 452 453 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 454 { 455 G4Nucleon * aNucleon = TheInvolvedNucleon[i]; 456 457 Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus; 458 //G4cout<<i<<" "<<Xminus<<" "<<XminusSum<<G4endl; 459 XminusSum-=Xminus; 460 if((Xminus <= 0.) || (Xminus > 1.) || 461 (XminusSum <=0.) || (XminusSum > 1.)) {Success=false; break;} 462 463 G4double Px=aNucleon->Get4Momentum().px() - DeltaX; 464 G4double Py=aNucleon->Get4Momentum().py() - DeltaY; 465 466 M2target +=(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()* 467 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() + 468 Px*Px + Py*Py)/Xminus; 469 470 G4LorentzVector tmp(Px,Py,Xminus,0.); 471 aNucleon->SetMomentum(tmp); 472 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 473 474 if(Success && (ResidualMassNumber != 0)) 475 { 476 M2target +=(ResidualMass*ResidualMass + PtSum.mag2())/XminusSum; 477 } 478 //G4cout<<"Success "<<Success<<G4endl; 479 //G4int Uzhi; G4cin>>Uzhi; 480 } while(Success == 0); 481 482 //------------------------------------------------------------- 483 //G4cout<<"SqrtS > Mprojectile + std::sqrt(M2target) "<<SqrtS<<" "<<Mprojectile<<" "<< std::sqrt(M2target)<<" "<<Mprojectile + std::sqrt(M2target)<<G4endl; 484 //G4int Uzhi3; G4cin>>Uzhi3; 485 } while (SqrtS < Mprojectile + std::sqrt(M2target)); 486 487 //------------------------------------------------------------- 488 G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target 489 -2.*S*M2projectile - 2.*S*M2target 490 -2.*M2projectile*M2target; 491 492 G4double WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS; 493 G4double WplusProjectile=SqrtS - M2target/WminusTarget; 494 //G4cout<<"DecM W+ W- "<<DecayMomentum2<<" "<<WplusProjectile<<" "<<WminusTarget<<G4endl; 495 //------------------------------------------------------------- 496 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile; 497 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile; 498 Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile); 499 500 Pprojectile.transform(toLab); // The work with the projectile 501 primary->Set4Momentum(Pprojectile); // is finished at the moment. 502 //G4cout<<"Proj After "<<Pprojectile<<G4endl; 503 //------------------------------------------------------------- 504 //Ninvolv=0; 505 506 G4ThreeVector Residual3Momentum(0.,0.,1.); 507 508 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 509 { 510 G4Nucleon * aNucleon = TheInvolvedNucleon[i]; 511 G4LorentzVector tmp=aNucleon->Get4Momentum(); 512 Residual3Momentum-=tmp.vect(); 513 514 G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+ 515 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()* 516 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass(); 517 G4double Xminus=tmp.z(); 518 519 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus); 520 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus); 521 522 tmp.setPz(Pz); 523 tmp.setE(E); 524 tmp.transform(toLab); 525 //G4cout<<"invol "<<Ninvolv<<" "<<tmp<<G4endl; 526 //Ninvolv++; 527 aNucleon->SetMomentum(tmp); 528 G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron(); 529 targetSplitable->Set4Momentum(tmp); 530 //G4cout<<"nucleon M "<<aNucleon->Get4Momentum()<<G4endl; 531 532 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 533 534 G4double Mt2Residual=sqr(ResidualMass) + 535 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.x()); 536 537 G4double PzResidual=-WminusTarget*Residual3Momentum.z()/2. + 538 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z()); 539 G4double EResidual = WminusTarget*Residual3Momentum.z()/2. + 540 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z()); 541 542 // G4LorentzVector Residual4Momentum(0.,0.,0.,0.); 543 Residual4Momentum.setPx(Residual3Momentum.x()); 544 Residual4Momentum.setPy(Residual3Momentum.y()); 545 Residual4Momentum.setPz(PzResidual); 546 Residual4Momentum.setE(EResidual); 547 Residual4Momentum.transform(toLab); 548 549 //G4cout<<"Return Nucleus"<<G4endl; 550 //------------------------------------------------------------- 551 //------------------------------------------------------------- 552 //------------------------------------------------------------- 553 //G4int Uzhi2; G4cin>>Uzhi2; 554 555 return true; 556 } 146 557 147 558 // ------------------------------------------------------------ 148 559 G4bool G4FTFModel::ExciteParticipants() 149 560 { 150 / *// Uzhi 29.03.08 For elastic Scatt.151 G4cout<<" In ExciteParticipants() "<<theParticipants.theInteractions.size()<<G4endl;152 G4cout<<" test Params Tot "<<theParameters->GetTotalCrossSection()<<G4endl;153 G4cout<<" test Params Ela "<<theParameters->GetElasticCrossSection()<<G4endl;561 // // Uzhi 29.03.08 For elastic Scatt. 562 //G4cout<<" In ExciteParticipants() "<<theParticipants.theInteractions.size()<<G4endl; 563 //G4cout<<" test Params Tot "<<theParameters->GetTotalCrossSection()<<G4endl; 564 //G4cout<<" test Params Ela "<<theParameters->GetElasticCrossSection()<<G4endl; 154 565 155 G4int counter=0; 156 */ // Uzhi 29.03.08 157 158 159 566 //G4int counter=0; 567 // // Uzhi 29.03.08 568 569 570 //G4int InterNumber=0; // Vova 571 572 G4bool Successfull=false; 573 theParticipants.StartLoop(); //Uzhi 26 July 09 574 575 // while (theParticipants.Next()&& (InterNumber < 3)) // Vova 160 576 while (theParticipants.Next()) 161 577 { 162 578 const G4InteractionContent & collision=theParticipants.GetInteraction(); 163 / *164 counter++;165 G4cout<<" Inter #"<<counter<<G4endl;166 */579 // 580 //counter++; 581 //G4cout<<" Int num "<<counter<<G4endl; 582 // 167 583 G4VSplitableHadron * projectile=collision.GetProjectile(); 168 584 G4VSplitableHadron * target=collision.GetTarget(); 169 170 // // Uzhi 29.03.08 171 G4bool Successfull; 585 // G4Nucleon * TargetNucleon=collision.GetTargetNucleon(); // Uzhi 16.07.09 586 // Uzhi 16.07.09 ---------------------------- 172 587 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt()) 173 { 588 { // Elastic scattering ------------------------- 174 589 //G4cout<<"Elastic"<<G4endl; 175 Successfull=theElastic->ElasticScattering(projectile, target, theParameters); 590 if(theElastic->ElasticScattering(projectile, target, theParameters)) 591 { 592 Successfull = Successfull || true; 593 } else 594 { 595 //G4cout<<"Elastic Not succes"<<G4endl; 596 Successfull = Successfull || false; 597 target->SetStatus(2); 598 /* 599 if(target->GetStatus() == 0) // Uzhi 17.07.09 600 { 601 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09 602 TargetNucleon->Hit(aHit); // Uzhi 16.07.09 603 }; 604 */ 605 }; 176 606 } 177 607 else 178 { 179 //G4cout<<"Inelastic"<<G4endl; 180 Successfull=theExcitation->ExciteParticipants(projectile, target, theParameters); 608 { // Inelastic scattering ---------------------- 609 //G4cout<<"InElastic"<<G4endl; 610 if(theExcitation->ExciteParticipants(projectile, target, 611 theParameters, theElastic)) 612 { 613 Successfull = Successfull || true; 614 } else 615 { 616 //G4cout<<"InElastic Non succes"<<G4endl; 617 Successfull = Successfull || false; 618 target->SetStatus(2); 619 /* 620 if(target->GetStatus() == 0) // Uzhi 16.06.09 621 { 622 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09 623 TargetNucleon->Hit(aHit); // Uzhi 16.07.09 624 }; 625 */ 626 }; 181 627 } 182 // if(!Successfull) 183 // // Uzhi 29.03.08184 185 // if ( ! theExcitation->ExciteParticipants(projectile, target) ) 186 if(!Successfull) 187 { 628 } // end of the loop Uzhi 9.07.09 629 // Uzhi 16.07.09 ---------------------------- 630 631 if(!Successfull) 632 { 633 //G4cout<<"Process not successfull"<<G4endl; 188 634 // give up, clean up 189 190 std::vector<G4VSplitableHadron *> targets; 191 192 193 194 635 std::vector<G4VSplitableHadron *> primaries; 636 // std::vector<G4VSplitableHadron *> targets; // Uzhi 31.07.09 637 theParticipants.StartLoop(); // restart a loop 638 while ( theParticipants.Next() ) 639 { 640 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 195 641 // do not allow for duplicates ... 196 197 198 199 200 642 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 643 interaction.GetProjectile()) ) 644 primaries.push_back(interaction.GetProjectile()); 645 /* // Uzhi 31.07.09 646 if ( targets.end() == std::find(targets.begin(), targets.end(), 201 647 interaction.GetTarget()) ) 202 targets.push_back(interaction.GetTarget()); 203 } 204 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 205 primaries.clear(); 206 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron()); 207 targets.clear(); 208 209 return false; 210 } // End of the loop Uzhi 211 } 648 targets.push_back(interaction.GetTarget()); 649 */ // Uzhi 31.07.09 650 } 651 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 652 primaries.clear(); 653 /* // Uzhi 31.07.09 654 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron()); 655 targets.clear(); 656 */ // Uzhi 31.07.09 657 // theParticipants.theNucleus->StartLoop(); 658 659 //G4cout<<"NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl; 660 G4VSplitableHadron * aNucleon = 0; 661 for(G4int i=0; i < NumberOfInvolvedNucleon; i++) 662 { 663 aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron(); 664 if(aNucleon) 665 { 666 if(aNucleon->GetStatus() != 0) delete aNucleon; 667 // if(aNucleon->GetStatus() == 2) DeleteVSplitableHadron()(aNucleon); 668 } 669 } 670 671 NumberOfInvolvedNucleon=0; 672 //G4cout<<"Out of Excit"<<G4endl; G4int Uzhi; G4cin>>Uzhi; 673 return false; 674 } // End of if(!Successfull) 675 212 676 return true; 213 677 } … … 217 681 // Loop over all collisions; find all primaries, and all target ( targets may 218 682 // be duplicate in the List ( to unique G4VSplitableHadrons) 683 684 //G4cout<<"In build string"<<G4endl; 219 685 220 686 G4ExcitedStringVector * strings; … … 223 689 std::vector<G4VSplitableHadron *> primaries; 224 690 std::vector<G4VSplitableHadron *> targets; 691 std::vector<G4Nucleon *> TargetNucleons; // Uzhi 16.07.09 225 692 693 G4ExcitedString * FirstString(0); // If there will be a kink, 694 G4ExcitedString * SecondString(0); // two strings will be prodused. 695 226 696 theParticipants.StartLoop(); // restart a loop 697 //G4int InterCount(0); // Uzhi 227 698 while ( theParticipants.Next() ) 228 699 { 229 700 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 230 701 // do not allow for duplicates ... 702 231 703 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 232 704 interaction.GetProjectile()) ) … … 236 708 interaction.GetTarget()) ) 237 709 targets.push_back(interaction.GetTarget()); 710 711 if ( TargetNucleons.end() == std::find(TargetNucleons.begin(), 712 TargetNucleons.end(), 713 interaction.GetTargetNucleon()) ) 714 TargetNucleons.push_back(interaction.GetTargetNucleon()); 715 //InterCount++; 238 716 } 239 717 240 718 241 // G4cout << "BuildStrings prim/targ " << primaries.size() << " , " << 242 // targets.size() << G4endl; 719 //G4cout << "BuildStrings prim/targ/woundN " << primaries.size() << " , " <<targets.size() <<", "<<TargetNucleons.size()<< G4endl; 243 720 244 721 unsigned int ahadron; 245 722 // Only for hA-interactions Uzhi ------------------------------------- 723 //G4int StringN(0); 724 //G4cout<<"Proj strings -----------------------"<<G4endl; 246 725 for ( ahadron=0; ahadron < primaries.size() ; ahadron++) 247 726 { 248 //G4ThreeVector aPosition=primaries[ahadron]->GetPosition(); 249 //G4cout<<"Proj Build "<<aPosition<<" "<<primaries[ahadron]->GetTimeOfCreation()<<G4endl; 250 G4bool isProjectile=true; 251 strings->push_back(theExcitation->String(primaries[ahadron], isProjectile)); 727 //G4cout<<" string# "<<ahadron<<" "<<primaries[ahadron]->Get4Momentum()<<G4endl; 728 G4bool isProjectile(0); 729 if(primaries[ahadron]->GetStatus() == 1) {isProjectile=true; } 730 if(primaries[ahadron]->GetStatus() == 3) {isProjectile=false;} 731 732 FirstString=0; SecondString=0; 733 theExcitation->CreateStrings(primaries[ahadron], isProjectile, 734 FirstString, SecondString, 735 theParameters); 736 //G4cout<<"1str 2str "<<FirstString<<" "<<SecondString<<G4endl; 737 if(FirstString != 0) strings->push_back(FirstString); 738 if(SecondString != 0) strings->push_back(SecondString); 739 740 //StringN++; G4cout<<"Proj string "<<strings->size()<<G4endl; 252 741 } 253 742 //G4cout<<"Targ strings ------------------------------"<<G4endl; 254 743 for ( ahadron=0; ahadron < targets.size() ; ahadron++) 255 744 { 256 //G4ThreeVector aPosition=targets[ahadron]->GetPosition(); 257 //G4cout<<"Targ Build "<<aPosition<<" "<<targets[ahadron]->GetTimeOfCreation()<<G4endl; 258 G4bool isProjectile=false; 259 strings->push_back(theExcitation->String(targets[ahadron], isProjectile)); 745 //G4cout<<"targets[ahadron]->GetStatus() "<<targets[ahadron]->GetStatus()<<G4endl; 746 if(targets[ahadron]->GetStatus() == 1) // Uzhi 17.07.09 747 { 748 G4bool isProjectile=false; 749 FirstString=0; SecondString=0; 750 theExcitation->CreateStrings(targets[ahadron], isProjectile, 751 FirstString, SecondString, 752 theParameters); 753 if(FirstString != 0) strings->push_back(FirstString); 754 if(SecondString != 0) strings->push_back(SecondString); 755 756 //StringN++; G4cout<<"Targ string"<<StringN<<G4endl; 757 } else 758 { 759 if(targets[ahadron]->GetStatus() == 0)// Uzhi 17.07.09 Nucleon was rejected 760 { 761 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09 762 TargetNucleons[ahadron]->Hit(aHit); // Uzhi 16.07.09 763 } 764 } 260 765 } 766 767 //G4cout<<"Proj + Targ string "<<strings->size()<<G4endl; 768 //G4cout<<"Involv strings NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl; 769 for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++) 770 { 771 /* 772 G4cout<<"ahadron "<<ahadron<<" Status "<< 773 TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus()<< 774 TheInvolvedNucleon[ahadron]->GetSplitableHadron()->Get4Momentum()<<G4endl; 775 */ 776 if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() == 2) 777 { 778 //StringN++; G4cout<<"Invol string "<<StringN<<G4endl; 779 G4bool isProjectile=false; 780 FirstString=0; SecondString=0; 781 theExcitation->CreateStrings( 782 TheInvolvedNucleon[ahadron]->GetSplitableHadron(), 783 isProjectile, 784 FirstString, SecondString, 785 theParameters); 786 if(FirstString != 0) strings->push_back(FirstString); 787 if(SecondString != 0) strings->push_back(SecondString); 788 789 // strings->push_back(theExcitation->String( 790 // TheInvolvedNucleon[ahadron]->GetSplitableHadron(), 791 // isProjectile)); 792 } 793 //G4cout<<"Proj + Targ+Involved string "<<strings->size()<<G4endl; 794 /* 795 else 796 { 797 G4cout<<"Else ????????????"<<G4endl; 798 if(targets[ahadron]->GetStatus() == 0)// Uzhi 17.07.09 Nucleon was rejected 799 { 800 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09 801 TargetNucleons[ahadron]->Hit(aHit); // Uzhi 16.07.09 802 } 803 } 804 */ 805 } 806 807 //G4int Uzhi; G4cin>>Uzhi; 261 808 262 809 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 263 810 primaries.clear(); 811 264 812 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron()); 265 813 targets.clear(); 814 815 return strings; 816 } 817 // ------------------------------------------------------------ 818 void G4FTFModel::GetResidualNucleus() 819 { // This method is needed for the correct application of G4PrecompoundModelInterface 820 G4double DeltaExcitationE=ResidualExcitationEnergy/ 821 (G4double) NumberOfInvolvedNucleon; 822 G4LorentzVector DeltaPResidualNucleus = Residual4Momentum/ 823 (G4double) NumberOfInvolvedNucleon; 824 825 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 826 { 827 G4Nucleon * aNucleon = TheInvolvedNucleon[i]; 828 G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus; 829 aNucleon->SetMomentum(tmp); 830 aNucleon->SetBindingEnergy(DeltaExcitationE); 831 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 832 833 } 834 835 // ------------------------------------------------------------ 836 G4ThreeVector G4FTFModel::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const 837 { // @@ this method is used in FTFModel as well. Should go somewhere common! 266 838 267 return strings; 268 } 269 // ------------------------------------------------------------ 839 G4double Pt2; 840 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 841 (std::exp(-maxPtSquare/AveragePt2)-1.)); 842 843 G4double Pt=std::sqrt(Pt2); 844 G4double phi=G4UniformRand() * twopi; 845 846 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 847 } -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParameters.cc
r1007 r1196 1 // 2 // ******************************************************************** 1 //******************************************************************** 3 2 // * License and Disclaimer * 4 3 // * * … … 25 24 // 26 25 // 27 // $Id: G4FTFParameters.cc,v 1. 4 2008/12/18 13:02:00 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4FTFParameters.cc,v 1.11 2009/10/25 10:50:54 vuzhinsk Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 28 // 30 29 31 30 #include "G4FTFParameters.hh" 31 32 #include "G4ios.hh" 33 #include <utility> // Uzhi 29.03.08 32 34 33 35 G4FTFParameters::G4FTFParameters() … … 43 45 G4double theZ, 44 46 G4double s) 45 47 { 46 48 G4int PDGcode = particle->GetPDGEncoding(); 47 49 G4int absPDGcode = std::abs(PDGcode); 48 G4double Elab = (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV; 49 G4double Plab = std::sqrt(Elab * Elab - 0.88); 50 G4double ProjectileMass = particle->GetPDGMass(); 51 G4double TargetMass = G4Proton::Proton()->GetPDGMass(); 52 53 G4double Elab = (s - ProjectileMass*ProjectileMass - TargetMass*TargetMass)/ 54 (2*TargetMass); 55 G4double Plab = std::sqrt(Elab * Elab - ProjectileMass*ProjectileMass); 50 56 51 57 G4double LogPlab = std::log( Plab ); 52 58 G4double sqrLogPlab = LogPlab * LogPlab; 53 54 //G4cout<<"G4FTFParameters Plab "<<Plab<<G4endl;55 59 56 60 G4int NumberOfTargetProtons = (G4int) theZ; … … 144 148 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons; 145 149 } 146 else if( PDGcode == 311 ) //------Projectile is KaonZero ------150 else if((PDGcode == 311) || (PDGcode == 130) || (PDGcode == 310))//Projectile is KaonZero 147 151 { 148 152 G4double XtotKP =( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+ … … 181 185 SetInelasticCrossSection(Xtotal-Xelastic); 182 186 183 // // Interactions with elastic an sinelastic collisions187 // // Interactions with elastic and inelastic collisions 184 188 SetProbabilityOfElasticScatt(Xtotal, Xelastic); 185 189 SetRadiusOfHNinteractions2(Xtotal/pi/10.); … … 190 194 */ //======================================================= 191 195 192 //G4cout<<" Rnn "<<Xtotal/pi/10.<<" "<<Xtotal/pi/10.*fermi*fermi<<G4endl;193 //G4cout<<"G4FTFParameters Xt Xel MeV "<<Xtotal<<" "<<Xelastic<<" "<<GeV<<G4endl;194 195 196 //----------------------------------------------------------------------------------- 196 197 SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering 197 198 // (GeV/c)^(-2)) 198 //G4cout<<"G4FTFParameters Slope "<<GetSlope()<<G4endl;199 199 //----------------------------------------------------------------------------------- 200 200 SetGamma0( GetSlope()*Xtotal/10./2./pi ); … … 208 208 if( absPDGcode > 1000 ) //------Projectile is baryon -------- 209 209 { 210 SetMagQuarkExchange(3.4); //3.8); 211 SetSlopeQuarkExchange(1.2); 212 SetDeltaProbAtQuarkExchange(0.1); //(0.1*4.); 213 210 214 SetProjMinDiffMass(1.1); // GeV 211 215 SetProjMinNonDiffMass(1.1); // GeV 212 SetProbabilityOfProjDiff(0. 95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel216 SetProbabilityOfProjDiff(0.76*std::pow(s/GeV/GeV,-0.35)); 213 217 214 218 SetTarMinDiffMass(1.1); // GeV 215 219 SetTarMinNonDiffMass(1.1); // GeV 216 SetProbabilityOfTarDiff(0. 95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel220 SetProbabilityOfTarDiff(0.76*std::pow(s/GeV/GeV,-0.35)); 217 221 218 222 SetAveragePt2(0.3); // GeV^2 … … 220 224 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion ----------- 221 225 { 226 SetMagQuarkExchange(120.); // 210. 227 SetSlopeQuarkExchange(2.0); 228 SetDeltaProbAtQuarkExchange(0.6); 229 222 230 SetProjMinDiffMass(0.5); // GeV 223 231 SetProjMinNonDiffMass(0.3); // GeV 224 SetProbabilityOfProjDiff(0. 62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel225 232 SetProbabilityOfProjDiff(0.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel 233 //G4cout<<"Params "<<0.6*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl; 226 234 SetTarMinDiffMass(1.1); // GeV 227 235 SetTarMinNonDiffMass(1.1); // GeV 228 SetProbabilityOfTarDiff(0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel 236 //G4cout<<" "<<2.7*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl; // was 2 237 //G4int Uzhi; G4cin>>Uzhi; 238 SetProbabilityOfTarDiff(2.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel 229 239 230 240 /* … … 236 246 SetAveragePt2(0.3); // GeV^2 237 247 } 238 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon ----------- 248 else if( (absPDGcode == 321) || (PDGcode == 311) || 249 (PDGcode == 130) || (PDGcode == 310)) //Projectile is Kaon 239 250 { 251 // Must be corrected, taken from PiN 252 SetMagQuarkExchange(120.); 253 SetSlopeQuarkExchange(2.0); 254 SetDeltaProbAtQuarkExchange(0.6); 255 240 256 SetProjMinDiffMass(0.7); // GeV 1.1 241 257 SetProjMinNonDiffMass(0.7); // GeV … … 251 267 //------Nucleon assumed 252 268 { 269 SetMagQuarkExchange(3.5); 270 SetSlopeQuarkExchange(1.0); 271 SetDeltaProbAtQuarkExchange(0.1); 272 253 273 SetProjMinDiffMass((particle->GetPDGMass()+160.*MeV)/GeV); 254 274 SetProjMinNonDiffMass((particle->GetPDGMass()+160.*MeV)/GeV); … … 260 280 261 281 SetAveragePt2(0.3); // GeV^2 262 }; 263 264 265 //G4cout<<"G4FTFParameters Out"<<G4endl; 266 267 } 282 } 283 284 // ---------- Set parameters of a string kink ------------------------------- 285 SetPt2Kink(6.*GeV*GeV); 286 G4double Puubar(1./3.), Pddbar(1./3.), Pssbar(1./3.); // SU(3) symmetry 287 // G4double Puubar(0.41 ), Pddbar(0.41 ), Pssbar(0.18 ); // Broken SU(3) symmetry 288 SetQuarkProbabilitiesAtGluonSplitUp(Puubar, Pddbar, Pssbar); 289 290 // --------- Set parameters of nuclear destruction-------------------- 291 292 if( absPDGcode < 1000 ) 293 { 294 SetCofNuclearDestruction(1.); //1.0); // for meson projectile 295 } else if( theA > 20. ) 296 { 297 SetCofNuclearDestruction(0.2); //2); // for baryon projectile and heavy target 298 } else 299 { 300 SetCofNuclearDestruction(0.2); //1.0); // for baryon projectile and light target 301 } 302 303 SetR2ofNuclearDestruction(1.5*fermi*fermi); 304 305 SetExcitationEnergyPerWoundedNucleon(100*MeV); 306 307 SetDofNuclearDestruction(0.4); 308 SetPt2ofNuclearDestruction(0.17*GeV*GeV); 309 SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV); 310 } 268 311 //********************************************************************************************** -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4FTFParticipants.cc,v 1. 9 2008/06/13 12:49:23vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4FTFParticipants.cc,v 1.15 2009/10/25 10:50:54 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------ … … 87 87 xyradius =theNucleus->GetOuterRadius() + deltaxy; // Impact parameter sampling 88 88 // radius 89 G4bool nucleusNeedsShift = true; 89 // G4bool nucleusNeedsShift = true; // Uzhi 20 July 2009 90 90 91 91 while ( theInteractions.size() == 0 ) … … 96 96 G4double impactY = theImpactParameter.second; 97 97 98 G4ThreeVector thePosition(impactX, impactY, -DBL_MAX); //Uzhi 24 July 09 99 primarySplitable->SetPosition(thePosition); //Uzhi 24 July 09 100 98 101 theNucleus->StartLoop(); 99 102 G4Nucleon * nucleon; 100 103 //G4int InterNumber=0; // Uzhi 104 G4int NucleonNumber(0); // Uzhi 101 105 //while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi 102 106 while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi 103 107 { 108 //G4cout<<"Nucl# "<<NucleonNumber<<G4endl; // Vova 104 109 G4double impact2= sqr(impactX - nucleon->GetPosition().x()) + 105 110 sqr(impactY - nucleon->GetPosition().y()); … … 110 115 { 111 116 //InterNumber++; 117 /* // Uzhi 20 July 2009 112 118 if ( nucleusNeedsShift ) 113 119 { // on the first hit, shift nucleus … … 117 123 impactY=0; 118 124 } 119 G4VSplitableHadron * targetSplitable; 120 if ( (targetSplitable=nucleon->GetSplitableHadron()) == NULL ) 125 */ // Uzhi 20 July 2009 126 //G4cout<<" Interact"<<G4endl; 127 primarySplitable->SetStatus(1); // It takes part in the interaction 128 129 G4VSplitableHadron * targetSplitable=0; 130 if ( ! nucleon->AreYouHit() ) 121 131 { 132 //G4cout<<"Part "<<nucleon->Get4Momentum()<<G4endl; 122 133 targetSplitable= new G4DiffractiveSplitableHadron(*nucleon); 123 134 nucleon->Hit(targetSplitable); 135 nucleon->SetBindingEnergy(3.*nucleon->GetBindingEnergy()); // Uzhi 5.10.09 136 targetSplitable->SetStatus(1); // It takes part in the interaction 124 137 } 125 138 G4InteractionContent * aInteraction = 126 139 new G4InteractionContent(primarySplitable); 127 140 aInteraction->SetTarget(targetSplitable); 141 aInteraction->SetTargetNucleon(nucleon); // Uzhi 16.07.09 128 142 theInteractions.push_back(aInteraction); 129 143 } 144 NucleonNumber++; // Uzhi 130 145 } 131 132 // G4cout << "Number of Hit nucleons " << theInteractions.size() // entries()146 // // Uzhi 147 // G4cout << "Number of Hit nucleons " << theInteractions.size()<<G4endl; // entries() 133 148 // << "\t" << impactX/fermi << "\t"<<impactY/fermi 134 149 // << "\t" << std::sqrt(sqr(impactX)+sqr(impactY))/fermi <<G4endl; 135 150 // // Uzhi 136 151 } 137 152
Note: See TracChangeset
for help on using the changeset viewer.