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
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "globals.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "Randomize.hh"
00050
00051 #include "G4QGSDiffractiveExcitation.hh"
00052 #include "G4LorentzRotation.hh"
00053 #include "G4ThreeVector.hh"
00054 #include "G4ParticleDefinition.hh"
00055 #include "G4VSplitableHadron.hh"
00056 #include "G4ExcitedString.hh"
00057
00058
00059 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation()
00060 {
00061 }
00062
00063 G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation()
00064 {
00065 }
00066
00067
00068 G4bool G4QGSDiffractiveExcitation::
00069 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
00070 {
00071
00072 G4LorentzVector Pprojectile=projectile->Get4Momentum();
00073
00074
00075 G4bool PutOnMassShell=0;
00076
00077
00078 G4double M0projectile = Pprojectile.mag();
00079
00080 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
00081 {
00082 PutOnMassShell=1;
00083 M0projectile=projectile->GetDefinition()->GetPDGMass();
00084 }
00085
00086 G4double Mprojectile2 = M0projectile * M0projectile;
00087
00088 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();
00089 G4int absPDGcode=std::abs(PDGcode);
00090 G4double ProjectileDiffCut;
00091 G4double AveragePt2;
00092
00093 if( absPDGcode > 1000 )
00094 {
00095 ProjectileDiffCut = 1.1;
00096 AveragePt2 = 0.3;
00097 }
00098 else if( absPDGcode == 211 || PDGcode == 111)
00099 {
00100 ProjectileDiffCut = 1.0;
00101 AveragePt2 = 0.3;
00102 }
00103 else if( absPDGcode == 321 || PDGcode == -311)
00104 {
00105 ProjectileDiffCut = 1.1;
00106 AveragePt2 = 0.3;
00107 }
00108 else
00109 {
00110 ProjectileDiffCut = 1.1;
00111 AveragePt2 = 0.3;
00112 };
00113
00114 ProjectileDiffCut = ProjectileDiffCut * GeV;
00115 AveragePt2 = AveragePt2 * GeV*GeV;
00116
00117
00118 G4LorentzVector Ptarget=target->Get4Momentum();
00119
00120 G4double M0target = Ptarget.mag();
00121
00122 if(M0target < target->GetDefinition()->GetPDGMass())
00123 {
00124 PutOnMassShell=1;
00125 M0target=target->GetDefinition()->GetPDGMass();
00126 }
00127
00128 G4double Mtarget2 = M0target * M0target;
00129
00130 G4double NuclearNucleonDiffCut = 1.1*GeV;
00131
00132 G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
00133 G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
00134
00135
00136
00137 G4LorentzVector Psum;
00138 Psum=Pprojectile+Ptarget;
00139
00140 G4LorentzRotation toCms(-1*Psum.boostVector());
00141
00142 G4LorentzVector Ptmp=toCms*Pprojectile;
00143
00144 if ( Ptmp.pz() <= 0. )
00145 {
00146
00147
00148 return false;
00149 }
00150
00151 toCms.rotateZ(-1*Ptmp.phi());
00152 toCms.rotateY(-1*Ptmp.theta());
00153
00154 G4LorentzRotation toLab(toCms.inverse());
00155
00156 Pprojectile.transform(toCms);
00157 Ptarget.transform(toCms);
00158
00159 G4double Pt2;
00160 G4double ProjMassT2, ProjMassT;
00161 G4double TargMassT2, TargMassT;
00162 G4double PZcms2, PZcms;
00163 G4double PMinusNew, TPlusNew;
00164
00165 G4double S=Psum.mag2();
00166 G4double SqrtS=std::sqrt(S);
00167
00168 if(SqrtS < 2200*MeV) {return false;}
00169
00170
00171 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
00172 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
00173 if(PZcms2 < 0)
00174 {return false;}
00175
00176 PZcms = std::sqrt(PZcms2);
00177
00178 if(PutOnMassShell)
00179 {
00180 if(Pprojectile.z() > 0.)
00181 {
00182 Pprojectile.setPz( PZcms);
00183 Ptarget.setPz( -PZcms);
00184 }
00185 else
00186 {
00187 Pprojectile.setPz(-PZcms);
00188 Ptarget.setPz( PZcms);
00189 };
00190
00191 Pprojectile.setE(std::sqrt(Mprojectile2+
00192 Pprojectile.x()*Pprojectile.x()+
00193 Pprojectile.y()*Pprojectile.y()+
00194 PZcms2));
00195 Ptarget.setE(std::sqrt( Mtarget2 +
00196 Ptarget.x()*Ptarget.x()+
00197 Ptarget.y()*Ptarget.y()+
00198 PZcms2));
00199 }
00200
00201 G4double maxPtSquare = PZcms2;
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 G4LorentzVector Qmomentum;
00213 G4double Qminus, Qplus;
00214
00215
00216 G4int whilecount=0;
00217 do {
00218
00219
00220 if (whilecount++ >= 500 && (whilecount%100)==0)
00221
00222
00223
00224 if (whilecount > 1000 )
00225 {
00226 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00227 return false;
00228 }
00229
00230 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00252 ProjMassT2=Mprojectile2+Pt2;
00253 ProjMassT =std::sqrt(ProjMassT2);
00254
00255 TargMassT2=Mtarget2+Pt2;
00256 TargMassT =std::sqrt(TargMassT2);
00257
00258 PZcms2=(S*S+ProjMassT2*ProjMassT2+
00259 TargMassT2*TargMassT2-
00260 2.*S*ProjMassT2-2.*S*TargMassT2-
00261 2.*ProjMassT2*TargMassT2)/4./S;
00262 if(PZcms2 < 0 ) {PZcms2=0;};
00263 PZcms =std::sqrt(PZcms2);
00264
00265 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
00266 G4double PMinusMax=SqrtS-TargMassT;
00267
00268 PMinusNew=ChooseP(PMinusMin,PMinusMax);
00269 Qminus=PMinusNew-Pprojectile.minus();
00270
00271 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
00272 G4double TPlusMax=SqrtS-ProjMassT;
00273
00274 TPlusNew=ChooseP(TPlusMin, TPlusMax);
00275 Qplus=-(TPlusNew-Ptarget.plus());
00276
00277 Qmomentum.setPz( (Qplus-Qminus)/2 );
00278 Qmomentum.setE( (Qplus+Qminus)/2 );
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291 } while (( (Pprojectile+Qmomentum).mag2() < Mprojectile2 ||
00292 (Ptarget -Qmomentum).mag2() < Mtarget2 ) ||
00293 ( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 &&
00294 (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );
00295
00296 if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2)
00297 {
00298 G4double TMinusNew=SqrtS-PMinusNew;
00299 Qminus=Ptarget.minus()-TMinusNew;
00300 TPlusNew=TargMassT2/TMinusNew;
00301 Qplus=Ptarget.plus()-TPlusNew;
00302
00303 Qmomentum.setPz( (Qplus-Qminus)/2 );
00304 Qmomentum.setE( (Qplus+Qminus)/2 );
00305 }
00306 else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2)
00307 {
00308 G4double PPlusNew=SqrtS-TPlusNew;
00309 Qplus=PPlusNew-Pprojectile.plus();
00310 PMinusNew=ProjMassT2/PPlusNew;
00311 Qminus=PMinusNew-Pprojectile.minus();
00312
00313 Qmomentum.setPz( (Qplus-Qminus)/2 );
00314 Qmomentum.setE( (Qplus+Qminus)/2 );
00315 };
00316
00317 Pprojectile += Qmomentum;
00318 Ptarget -= Qmomentum;
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337 Pprojectile.transform(toLab);
00338 Ptarget.transform(toLab);
00339
00340
00341
00342
00343
00344
00345 target->Set4Momentum(Ptarget);
00346
00347
00348
00349 projectile->Set4Momentum(Pprojectile);
00350
00351 return true;
00352 }
00353
00354
00355 G4ExcitedString * G4QGSDiffractiveExcitation::
00356 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
00357 {
00358 hadron->SplitUp();
00359 G4Parton *start= hadron->GetNextParton();
00360 if ( start==NULL)
00361 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
00362 return NULL;
00363 }
00364 G4Parton *end = hadron->GetNextParton();
00365 if ( end==NULL)
00366 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
00367 return NULL;
00368 }
00369
00370 G4ExcitedString * string;
00371 if ( isProjectile )
00372 {
00373 string= new G4ExcitedString(end,start, +1);
00374 } else {
00375 string= new G4ExcitedString(start,end, -1);
00376 }
00377
00378 string->SetPosition(hadron->GetPosition());
00379
00380
00381 G4double ptSquared= hadron->Get4Momentum().perp2();
00382 G4double transverseMassSquared= hadron->Get4Momentum().plus()
00383 * hadron->Get4Momentum().minus();
00384
00385
00386 G4double maxAvailMomentumSquared=
00387 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
00388
00389 G4double widthOfPtSquare = 0.25;
00390 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
00391
00392 G4LorentzVector Pstart(G4LorentzVector(pt,0.));
00393 G4LorentzVector Pend;
00394 Pend.setPx(hadron->Get4Momentum().px() - pt.x());
00395 Pend.setPy(hadron->Get4Momentum().py() - pt.y());
00396
00397 G4double tm1=hadron->Get4Momentum().minus() +
00398 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
00399
00400 G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
00401 4. * Pend.perp2() * hadron->Get4Momentum().minus()
00402 / hadron->Get4Momentum().plus() ));
00403
00404 G4int Sign= isProjectile ? -1 : 1;
00405
00406 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
00407 G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
00408
00409 G4double startPlus= Pstart.perp2() / startMinus;
00410 G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
00411
00412 Pstart.setPz(0.5*(startPlus - startMinus));
00413 Pstart.setE(0.5*(startPlus + startMinus));
00414
00415 Pend.setPz(0.5*(endPlus - endMinus));
00416 Pend.setE(0.5*(endPlus + endMinus));
00417
00418 start->Set4Momentum(Pstart);
00419 end->Set4Momentum(Pend);
00420
00421 #ifdef G4_FTFDEBUG
00422 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
00423 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
00424 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
00425 G4cout << " sum of ends " << Pstart+Pend << G4endl;
00426 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
00427 #endif
00428
00429 return string;
00430 }
00431
00432
00433
00434
00435 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
00436 {
00437
00438
00439
00440
00441 G4double range=Pmax-Pmin;
00442
00443 if ( Pmin <= 0. || range <=0. )
00444 {
00445 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
00446 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
00447 }
00448
00449 G4double P;
00450
00451
00452
00453
00454
00455
00456 P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());
00457
00458
00459 return P;
00460 }
00461
00462 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
00463 {
00464
00465 G4double Pt2;
00466
00467
00468
00469
00470
00471
00472 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));
00473
00474 G4double Pt=std::sqrt(Pt2);
00475
00476 G4double phi=G4UniformRand() * twopi;
00477
00478 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
00479 }