00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include "G4VLongitudinalStringDecay.hh"
00036 #include "G4PhysicalConstants.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4ios.hh"
00039 #include "Randomize.hh"
00040 #include "G4FragmentingString.hh"
00041
00042 #include "G4ParticleDefinition.hh"
00043 #include "G4ParticleTypes.hh"
00044 #include "G4ParticleChange.hh"
00045 #include "G4VShortLivedParticle.hh"
00046 #include "G4ShortLivedConstructor.hh"
00047 #include "G4ParticleTable.hh"
00048 #include "G4ShortLivedTable.hh"
00049 #include "G4PhaseSpaceDecayChannel.hh"
00050 #include "G4VDecayChannel.hh"
00051 #include "G4DecayTable.hh"
00052
00053 #include "G4DiQuarks.hh"
00054 #include "G4Quarks.hh"
00055 #include "G4Gluons.hh"
00056
00057
00058
00059
00060
00061
00062
00063
00064 G4VLongitudinalStringDecay::G4VLongitudinalStringDecay()
00065 {
00066 MassCut = 0.35*GeV;
00067 ClusterMass = 0.15*GeV;
00068
00069 SmoothParam = 0.9;
00070 StringLoopInterrupt = 1000;
00071 ClusterLoopInterrupt = 500;
00072
00073
00074 SigmaQT = 0.5 * GeV;
00075
00076 StrangeSuppress = 0.44;
00077 DiquarkSuppress = 0.07;
00078 DiquarkBreakProb = 0.1;
00079
00080
00081 pspin_meson = 0.5;
00082
00083
00084 pspin_barion = 0.5;
00085
00086
00087 vectorMesonMix.resize(6);
00088 vectorMesonMix[0] = 0.5;
00089 vectorMesonMix[1] = 0.0;
00090 vectorMesonMix[2] = 0.5;
00091 vectorMesonMix[3] = 0.0;
00092 vectorMesonMix[4] = 1.0;
00093 vectorMesonMix[5] = 1.0;
00094
00095
00096 scalarMesonMix.resize(6);
00097 scalarMesonMix[0] = 0.5;
00098 scalarMesonMix[1] = 0.25;
00099 scalarMesonMix[2] = 0.5;
00100 scalarMesonMix[3] = 0.25;
00101 scalarMesonMix[4] = 1.0;
00102 scalarMesonMix[5] = 0.5;
00103
00104
00105 PastInitPhase=false;
00106 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
00107 scalarMesonMix,vectorMesonMix);
00108 Kappa = 1.0 * GeV/fermi;
00109
00110
00111 }
00112
00113
00114 G4VLongitudinalStringDecay::~G4VLongitudinalStringDecay()
00115 {
00116 delete hadronizer;
00117 }
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129 int G4VLongitudinalStringDecay::operator==(const G4VLongitudinalStringDecay &) const
00130 {
00131 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator== forbidden");
00132 return false;
00133 }
00134
00135
00136
00137 int G4VLongitudinalStringDecay::operator!=(const G4VLongitudinalStringDecay &) const
00138 {
00139 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator!= forbidden");
00140 return true;
00141 }
00142
00143
00144
00145
00146 void G4VLongitudinalStringDecay::SetMassCut(G4double aValue){MassCut=aValue;}
00147
00148
00149
00150
00151
00152 G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const
00153 G4ExcitedString * const string)
00154 {
00155
00156
00157 G4KineticTrackVector * result=0;
00158
00159 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0);
00160
00161 G4FragmentingString aString(*string);
00162
00163 if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) {
00164 return 0;
00165 }
00166
00167
00168
00169 result=new G4KineticTrackVector;
00170
00171 if ( hadrons.second ==0 )
00172 {
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 G4ThreeVector Mom3 = string->Get4Momentum().vect();
00184 G4LorentzVector Mom(Mom3,
00185 std::sqrt(Mom3.mag2() +
00186 sqr(hadrons.first->GetPDGMass())));
00187 result->push_back(new G4KineticTrack(hadrons.first, 0,
00188 string->GetPosition(),
00189 Mom));
00190 } else
00191 {
00192
00193
00194 #ifdef DEBUG_LightFragmentationTest
00195 G4cout << "VlongSF Warning replacing qq-qqbar string by TWO hadrons "
00196 << hadrons.first->GetParticleName() << " / "
00197 << hadrons.second->GetParticleName()
00198 << "string .. " << string->Get4Momentum() << " "
00199 << string->Get4Momentum().m() << G4endl;
00200 #endif
00201
00202 G4LorentzVector Mom1, Mom2;
00203 Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(),
00204 &Mom2,hadrons.second->GetPDGMass(),
00205 string->Get4Momentum().mag());
00206
00207 result->push_back(new G4KineticTrack(hadrons.first, 0,
00208 string->GetPosition(),
00209 Mom1));
00210 result->push_back(new G4KineticTrack(hadrons.second, 0,
00211 string->GetPosition(),
00212 Mom2));
00213
00214 G4ThreeVector Velocity = string->Get4Momentum().boostVector();
00215 result->Boost(Velocity);
00216 }
00217
00218 return result;
00219
00220 }
00221
00222
00223
00224 G4double G4VLongitudinalStringDecay::FragmentationMass(
00225 const G4FragmentingString * const string,
00226 Pcreate build, pDefPair * pdefs )
00227 {
00228
00229 G4double mass;
00230 static G4bool NeedInit(true);
00231 static std::vector<double> nomix;
00232 static G4HadronBuilder * minMassHadronizer;
00233 if ( NeedInit )
00234 {
00235 NeedInit = false;
00236 nomix.resize(6);
00237 for ( G4int i=0; i<6 ; i++ ) nomix[i]=0;
00238
00239
00240 minMassHadronizer=hadronizer;
00241 }
00242
00243 if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin;
00244
00245 G4ParticleDefinition *Hadron1, *Hadron2=0;
00246
00247 if (!string->FourQuarkString() )
00248 {
00249
00250
00251
00252 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),
00253 string->GetRightParton());
00254
00255 mass= (Hadron1)->GetPDGMass();
00256 } else
00257 {
00258
00259
00260 G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
00261 if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;
00262
00263
00264 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),
00265 FindParticle(iflc) );
00266 Hadron2 = (minMassHadronizer->*build)(string->GetRightParton(),
00267 FindParticle(-iflc) );
00268 mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();
00269 }
00270
00271 if ( pdefs != 0 )
00272 {
00273 pdefs->first = Hadron1;
00274 pdefs->second = Hadron2;
00275 }
00276
00277 return mass;
00278 }
00279
00280
00281
00282 G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding)
00283 {
00284 G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding);
00285 if (ptr == NULL)
00286 {
00287 G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl;
00288 throw G4HadronicException(__FILE__, __LINE__, "Check your particle table");
00289 }
00290 return ptr;
00291 }
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in)
00316 {
00317 G4Parton *Left=new G4Parton(*in.GetLeftParton());
00318 G4Parton *Right=new G4Parton(*in.GetRightParton());
00319 return new G4ExcitedString(Left,Right,in.GetDirection());
00320 }
00321
00322
00323
00324 G4KineticTrack * G4VLongitudinalStringDecay::Splitup(
00325 G4FragmentingString *string,
00326 G4FragmentingString *&newString)
00327 {
00328
00329
00330 G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
00331 if (SideOfDecay < 0)
00332 {
00333 string->SetLeftPartonStable();
00334 } else
00335 {
00336 string->SetRightPartonStable();
00337 }
00338
00339 G4ParticleDefinition *newStringEnd;
00340 G4ParticleDefinition * HadronDefinition;
00341 if (string->DecayIsQuark())
00342 {
00343 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
00344 } else {
00345 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
00346 }
00347
00348
00349
00350
00351 newString=new G4FragmentingString(*string,newStringEnd);
00352
00353
00354 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
00355
00356 delete newString; newString=0;
00357
00358 G4KineticTrack * Hadron =0;
00359 if ( HadronMomentum != 0 ) {
00360
00361 G4ThreeVector Pos;
00362 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
00363
00364 newString=new G4FragmentingString(*string,newStringEnd,
00365 HadronMomentum);
00366
00367 delete HadronMomentum;
00368 }
00369
00370 return Hadron;
00371 }
00372
00373
00374
00375 G4ParticleDefinition *
00376 G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition*
00377 decay, G4ParticleDefinition *&created)
00378 {
00379 G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1;
00380
00381
00382 pDefPair QuarkPair = CreatePartonPair(IsParticle);
00383 created = QuarkPair.second;
00384 return hadronizer->Build(QuarkPair.first, decay);
00385
00386 }
00387
00388
00389
00390 G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup(
00391 G4ParticleDefinition* decay,
00392 G4ParticleDefinition *&created)
00393 {
00394
00395 if (G4UniformRand() < DiquarkBreakProb ){
00396
00397
00398 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
00399 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
00400 if (G4UniformRand() < 0.5)
00401 {
00402 G4int Swap = stableQuarkEncoding;
00403 stableQuarkEncoding = decayQuarkEncoding;
00404 decayQuarkEncoding = Swap;
00405 }
00406
00407 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
00408
00409 pDefPair QuarkPair = CreatePartonPair(IsParticle,false);
00410
00411 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
00412 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
00413 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
00414 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
00415 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
00416 created = FindParticle(NewDecayEncoding);
00417 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
00418 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
00419 return had;
00420
00421
00422 } else {
00423
00424
00425 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;
00426
00427 pDefPair QuarkPair = CreatePartonPair(IsParticle,false);
00428 created = QuarkPair.second;
00429
00430 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
00431 return had;
00432
00433 }
00434 }
00435
00436
00437
00438 G4int G4VLongitudinalStringDecay::SampleQuarkFlavor(void)
00439 {
00440 return (1 + (int)(G4UniformRand()/StrangeSuppress));
00441 }
00442
00443
00444
00445 G4VLongitudinalStringDecay::pDefPair G4VLongitudinalStringDecay::CreatePartonPair(G4int NeedParticle,G4bool AllowDiquarks)
00446 {
00447
00448
00449 if ( AllowDiquarks && G4UniformRand() < DiquarkSuppress )
00450 {
00451
00452 G4int q1 = SampleQuarkFlavor();
00453 G4int q2 = SampleQuarkFlavor();
00454 G4int spin = (q1 != q2 && G4UniformRand() <= 0.5)? 1 : 3;
00455
00456 G4int PDGcode = (std::max(q1,q2) * 1000 + std::min(q1,q2) * 100 + spin) * NeedParticle;
00457 return pDefPair (FindParticle(-PDGcode),FindParticle(PDGcode));
00458
00459
00460 } else {
00461
00462 G4int PDGcode=SampleQuarkFlavor()*NeedParticle;
00463 return pDefPair (FindParticle(PDGcode),FindParticle(-PDGcode));
00464 }
00465
00466 }
00467
00468
00469 G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt(G4double ptMax)
00470 {
00471 G4double Pt;
00472 if ( ptMax < 0 ) {
00473
00474 Pt = -std::log(G4UniformRand());
00475 } else {
00476
00477 Pt = -std::log(CLHEP::RandFlat::shoot(std::exp(-sqr(ptMax)/sqr(SigmaQT)), 1.));
00478 }
00479 Pt = SigmaQT * std::sqrt(Pt);
00480 G4double phi = 2.*pi*G4UniformRand();
00481 return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
00482 }
00483
00484
00485
00486 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons)
00487 {
00488
00489
00490
00491 G4double kappa = GetStringTensionParameter();
00492 for(size_t c1 = 0; c1 < Hadrons->size(); c1++)
00493 {
00494 G4double SumPz = 0;
00495 G4double SumE = 0;
00496 for(size_t c2 = 0; c2 < c1; c2++)
00497 {
00498 SumPz += Hadrons->operator[](c2)->Get4Momentum().pz();
00499 SumE += Hadrons->operator[](c2)->Get4Momentum().e();
00500 }
00501 G4double HadronE = Hadrons->operator[](c1)->Get4Momentum().e();
00502 G4double HadronPz = Hadrons->operator[](c1)->Get4Momentum().pz();
00503 Hadrons->operator[](c1)->SetFormationTime(
00504 (theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa)/c_light);
00505
00506 G4ThreeVector aPosition(0, 0,
00507 (theInitialStringMass - 2.*SumE - HadronE + HadronPz)/(2.*kappa));
00508 Hadrons->operator[](c1)->SetPosition(aPosition);
00509
00510 }
00511 }
00512
00513
00514
00515 void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue)
00516 {
00517 if ( PastInitPhase ) {
00518 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed");
00519 } else {
00520 SigmaQT = aValue;
00521 }
00522 }
00523
00524
00525
00526 void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue)
00527 {
00528 if ( PastInitPhase ) {
00529 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed");
00530 } else {
00531 StrangeSuppress = aValue;
00532 }
00533 }
00534
00535
00536
00537 void G4VLongitudinalStringDecay::SetDiquarkSuppression(G4double aValue)
00538 {
00539 if ( PastInitPhase ) {
00540 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkSuppression after FragmentString() not allowed");
00541 } else {
00542 DiquarkSuppress = aValue;
00543 }
00544 }
00545
00546
00547
00548 void G4VLongitudinalStringDecay::SetDiquarkBreakProbability(G4double aValue)
00549 {
00550 if ( PastInitPhase ) {
00551 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkBreakProbability after FragmentString() not allowed");
00552 } else {
00553 DiquarkBreakProb = aValue;
00554 }
00555 }
00556
00557
00558
00559 void G4VLongitudinalStringDecay::SetVectorMesonProbability(G4double aValue)
00560 {
00561 if ( PastInitPhase ) {
00562 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonProbability after FragmentString() not allowed");
00563 } else {
00564 pspin_meson = aValue;
00565 delete hadronizer;
00566 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
00567 scalarMesonMix,vectorMesonMix);
00568 }
00569 }
00570
00571
00572
00573 void G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability(G4double aValue)
00574 {
00575 if ( PastInitPhase ) {
00576 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability after FragmentString() not allowed");
00577 } else {
00578 pspin_barion = aValue;
00579 delete hadronizer;
00580 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
00581 scalarMesonMix,vectorMesonMix);
00582 }
00583 }
00584
00585
00586
00587 void G4VLongitudinalStringDecay::SetScalarMesonMixings(std::vector<G4double> aVector)
00588 {
00589 if ( PastInitPhase ) {
00590 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetScalarMesonMixings after FragmentString() not allowed");
00591 } else {
00592 if ( aVector.size() < 6 )
00593 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetScalarMesonMixings( argument Vector too small");
00594 scalarMesonMix[0] = aVector[0];
00595 scalarMesonMix[1] = aVector[1];
00596 scalarMesonMix[2] = aVector[2];
00597 scalarMesonMix[3] = aVector[3];
00598 scalarMesonMix[4] = aVector[4];
00599 scalarMesonMix[5] = aVector[5];
00600 delete hadronizer;
00601 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
00602 scalarMesonMix,vectorMesonMix);
00603 }
00604 }
00605
00606
00607
00608 void G4VLongitudinalStringDecay::SetVectorMesonMixings(std::vector<G4double> aVector)
00609 {
00610 if ( PastInitPhase ) {
00611 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonMixings after FragmentString() not allowed");
00612 } else {
00613 if ( aVector.size() < 6 )
00614 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonMixings( argument Vector too small");
00615 vectorMesonMix[0] = aVector[0];
00616 vectorMesonMix[1] = aVector[1];
00617 vectorMesonMix[2] = aVector[2];
00618 vectorMesonMix[3] = aVector[3];
00619 vectorMesonMix[4] = aVector[4];
00620 vectorMesonMix[5] = aVector[5];
00621 delete hadronizer;
00622 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
00623 scalarMesonMix,vectorMesonMix);
00624
00625 }
00626 }
00627
00628
00629 void G4VLongitudinalStringDecay::SetStringTensionParameter(G4double aValue)
00630 {
00631 Kappa = aValue * GeV/fermi;
00632 }
00633
00634