source: trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc@ 1315

Last change on this file since 1315 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 44.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4DiffractiveExcitation.cc,v 1.21 2009/12/15 19:14:31 vuzhinsk Exp $
28// ------------------------------------------------------------
29// GEANT 4 class implemetation file
30//
31// ---------------- G4DiffractiveExcitation --------------
32// by Gunter Folger, October 1998.
33// diffractive Excitation used by strings models
34// Take a projectile and a target
35// excite the projectile and target
36// Essential changed by V. Uzhinsky in November - December 2006
37// in order to put it in a correspondence with original FRITIOF
38// model. Variant of FRITIOF with nucleon de-excitation is implemented.
39// Other changes by V.Uzhinsky in May 2007 were introduced to fit
40// meson-nucleon interactions. Additional changes by V. Uzhinsky
41// were introduced in December 2006. They treat diffraction dissociation
42// processes more exactly.
43// ---------------------------------------------------------------------
44
45
46#include "globals.hh"
47#include "Randomize.hh"
48
49#include "G4DiffractiveExcitation.hh"
50#include "G4FTFParameters.hh"
51#include "G4ElasticHNScattering.hh"
52
53#include "G4LorentzRotation.hh"
54#include "G4RotationMatrix.hh"
55#include "G4ThreeVector.hh"
56#include "G4ParticleDefinition.hh"
57#include "G4VSplitableHadron.hh"
58#include "G4ExcitedString.hh"
59#include "G4ParticleTable.hh"
60#include "G4Neutron.hh"
61#include "G4ParticleDefinition.hh"
62
63//#include "G4ios.hh"
64//#include "UZHI_diffraction.hh"
65
66G4DiffractiveExcitation::G4DiffractiveExcitation()
67{
68}
69
70// ---------------------------------------------------------------------
71G4bool G4DiffractiveExcitation::
72 ExciteParticipants(G4VSplitableHadron *projectile,
73 G4VSplitableHadron *target,
74 G4FTFParameters *theParameters,
75 G4ElasticHNScattering *theElastic) const
76{
77// -------------------- Projectile parameters -----------------------
78 G4LorentzVector Pprojectile=projectile->Get4Momentum();
79
80 if(Pprojectile.z() < 0.)
81 {
82 target->SetStatus(2);
83 return false;
84 }
85
86 G4double ProjectileRapidity = Pprojectile.rapidity();
87
88 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
89 G4int absProjectilePDGcode=std::abs(ProjectilePDGcode);
90
91 G4bool PutOnMassShell(false);
92// G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
93 G4double M0projectile = Pprojectile.mag(); // Without de-excitation
94
95 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
96 {
97 PutOnMassShell=true;
98 M0projectile=projectile->GetDefinition()->GetPDGMass();
99 }
100
101 G4double M0projectile2 = M0projectile * M0projectile;
102
103 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
104 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
105 G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
106
107// -------------------- Target parameters -------------------------
108 G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
109 G4int absTargetPDGcode=std::abs(TargetPDGcode);
110
111 G4LorentzVector Ptarget=target->Get4Momentum();
112
113 G4double M0target = Ptarget.mag();
114
115 G4double TargetRapidity = Ptarget.rapidity();
116
117 if(M0target < target->GetDefinition()->GetPDGMass())
118 {
119 PutOnMassShell=true;
120 M0target=target->GetDefinition()->GetPDGMass();
121 }
122
123 G4double M0target2 = M0target * M0target;
124
125 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();
126 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();
127 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
128
129 G4double AveragePt2=theParameters->GetAveragePt2();
130
131 G4double ProbOfDiffraction=ProbProjectileDiffraction +
132 ProbTargetDiffraction;
133
134 G4double SumMasses=M0projectile+M0target+200.*MeV;
135
136// Kinematical properties of the interactions --------------
137 G4LorentzVector Psum; // 4-momentum in CMS
138 Psum=Pprojectile+Ptarget;
139 G4double S=Psum.mag2();
140
141// Transform momenta to cms and then rotate parallel to z axis;
142 G4LorentzRotation toCms(-1*Psum.boostVector());
143
144 G4LorentzVector Ptmp=toCms*Pprojectile;
145
146 if ( Ptmp.pz() <= 0. )
147 {
148 target->SetStatus(2);
149 // "String" moving backwards in CMS, abort collision !!
150 return false;
151 }
152
153 toCms.rotateZ(-1*Ptmp.phi());
154 toCms.rotateY(-1*Ptmp.theta());
155
156 G4LorentzRotation toLab(toCms.inverse());
157
158 Pprojectile.transform(toCms);
159 Ptarget.transform(toCms);
160
161 G4double PZcms2, PZcms;
162
163 G4double SqrtS=std::sqrt(S);
164
165 if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses))
166 {target->SetStatus(2); return false;} // The model cannot work for
167 // p+p-interactions
168 // at Plab < 1.62 GeV/c.
169
170 if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) &&
171 ((SqrtS < 1600*MeV) || (SqrtS < SumMasses)))
172 {target->SetStatus(2); return false;} // The model cannot work for
173 // Pi+p-interactions
174 // at Plab < 1. GeV/c.
175
176 if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) ||
177 (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
178 (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses)))
179 {target->SetStatus(2); return false;} // The model cannot work for
180 // K+p-interactions
181 // at Plab < ??? GeV/c. ???
182
183 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
184 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
185 /4./S;
186
187 if(PZcms2 < 0)
188 {target->SetStatus(2); return false;} // It can be in an interaction with
189 // off-shell nuclear nucleon
190
191 PZcms = std::sqrt(PZcms2);
192
193 if(PutOnMassShell)
194 {
195 if(Pprojectile.z() > 0.)
196 {
197 Pprojectile.setPz( PZcms);
198 Ptarget.setPz( -PZcms);
199 }
200 else
201 {
202 Pprojectile.setPz(-PZcms);
203 Ptarget.setPz( PZcms);
204 };
205
206 Pprojectile.setE(std::sqrt(M0projectile2 +
207 Pprojectile.x()*Pprojectile.x()+
208 Pprojectile.y()*Pprojectile.y()+
209 PZcms2));
210 Ptarget.setE(std::sqrt(M0target2 +
211 Ptarget.x()*Ptarget.x()+
212 Ptarget.y()*Ptarget.y()+
213 PZcms2));
214 }
215
216 G4double maxPtSquare; // = PZcms2;
217
218// Charge exchange can be possible for baryons -----------------
219
220// Getting the values needed for exchange ----------------------
221 G4double MagQuarkExchange =theParameters->GetMagQuarkExchange();
222 G4double SlopeQuarkExchange =theParameters->GetSlopeQuarkExchange();
223 G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
224
225// G4double NucleonMass=
226// (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();
227 G4double DeltaMass=
228 (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
229
230// Check for possible quark excjane -----------------------------------
231 if(G4UniformRand() < MagQuarkExchange*
232 std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity)))
233 {
234 G4int NewProjCode(0), NewTargCode(0);
235
236 G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
237
238// Projectile unpacking --------------------------
239 if(absProjectilePDGcode < 1000 )
240 { // projectile is meson -----------------
241 UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2);
242 } else
243 { // projectile is baryon ----------------
244 UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3);
245 } // End of the hadron's unpacking ----------
246
247// Target unpacking ------------------------------
248 G4int TargQ1(0), TargQ2(0), TargQ3(0);
249 UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3);
250
251// Sampling of exchanged quarks -------------------
252 G4int ProjExchangeQ(0);
253 G4int TargExchangeQ(0);
254
255 if(absProjectilePDGcode < 1000 )
256 { // projectile is meson -----------------
257 if(ProjQ1 > 0 ) // ProjQ1 is quark
258 {
259 ProjExchangeQ = ProjQ1;
260 if(ProjExchangeQ != TargQ1)
261 {
262 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;
263 } else
264 if(ProjExchangeQ != TargQ2)
265 {
266 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;
267 } else
268 {
269 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;
270 }
271 } else // ProjQ2 is quark
272 {
273 ProjExchangeQ = ProjQ2;
274 if(ProjExchangeQ != TargQ1)
275 {
276 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;
277 } else
278 if(ProjExchangeQ != TargQ2)
279 {
280 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;
281 } else
282 {
283 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;
284 }
285
286 } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark
287
288 G4int aProjQ1=std::abs(ProjQ1);
289 G4int aProjQ2=std::abs(ProjQ2);
290 if(aProjQ1 == aProjQ2) {NewProjCode = 111;} // Pi0-meson
291 else // |ProjQ1| # |ProjQ2|
292 {
293 if(aProjQ1 > aProjQ2) {NewProjCode = aProjQ1*100+aProjQ2*10+1;}
294 else {NewProjCode = aProjQ2*100+aProjQ1*10+1;}
295 }
296
297 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
298
299 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
300 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar
301
302 else if( target->GetDefinition()->GetPDGiIsospin() == 3 ) //Delta was the target
303 { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;} //Save Delta isobar
304 else {} // De-excite initial Delta isobar
305 }
306
307 else if((G4UniformRand() < DeltaProbAtQuarkExchange) && //Nucleon was the target
308 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar
309 else {} //Save initial nucleon
310// target->SetDefinition( // Fix 15.12.09
311// G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09
312 } else
313 { // projectile is baryon ----------------
314 G4double Same=0.; //0.3; //0.5;
315 G4bool ProjDeltaHasCreated(false);
316 G4bool TargDeltaHasCreated(false);
317
318 G4double Ksi=G4UniformRand();
319 if(G4UniformRand() < 0.5) // Sampling exchange quark from proj. or targ.
320 { // Sampling exchanged quark from the projectile ---
321 if( Ksi < 0.333333 )
322 {ProjExchangeQ = ProjQ1;}
323 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
324 {ProjExchangeQ = ProjQ2;}
325 else
326 {ProjExchangeQ = ProjQ3;}
327
328 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same))
329 {
330 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
331 } else
332 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same))
333 {
334 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
335 } else
336 {
337 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
338 }
339
340 if( Ksi < 0.333333 )
341 {ProjQ1=ProjExchangeQ;}
342 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
343 {ProjQ2=ProjExchangeQ;}
344 else
345 {ProjQ3=ProjExchangeQ;}
346
347 } else
348 { // Sampling exchanged quark from the target -------
349 if( Ksi < 0.333333 )
350 {TargExchangeQ = TargQ1;}
351 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
352 {TargExchangeQ = TargQ2;}
353 else
354 {TargExchangeQ = TargQ3;}
355
356 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same))
357 {
358 ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
359 } else
360 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same))
361 {
362 ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
363 } else
364 {
365 ProjExchangeQ = ProjQ3; ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
366 }
367
368 if( Ksi < 0.333333 )
369 {TargQ1=TargExchangeQ;}
370 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
371 {TargQ2=TargExchangeQ;}
372 else
373 {TargQ3=TargExchangeQ;}
374
375 } // End of sampling baryon
376
377 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
378
379 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
380 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
381 { if(G4UniformRand() > DeltaProbAtQuarkExchange)
382 {NewProjCode +=2; ProjDeltaHasCreated=true;}
383 else {NewProjCode +=0; ProjDeltaHasCreated=false;}
384 }
385 else // Projectile was Nucleon
386 {
387 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target))
388 {NewProjCode +=2; ProjDeltaHasCreated=true;}
389 else {NewProjCode +=0; ProjDeltaHasCreated=false;}
390 }
391
392 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
393
394 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;}
395 else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta
396 { if(G4UniformRand() > DeltaProbAtQuarkExchange)
397 {NewTargCode +=2; TargDeltaHasCreated=true;}
398 else {NewTargCode +=0; TargDeltaHasCreated=false;}
399 }
400 else // Target was Nucleon
401 {
402 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass))
403 {NewTargCode +=2; TargDeltaHasCreated=true;}
404 else {NewTargCode +=0; TargDeltaHasCreated=false;}
405 }
406
407 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
408 { // Nothing was changed! It is not right!?
409 }
410// Forming baryons --------------------------------------------------
411
412 } // End of if projectile is baryon ---------------------------
413
414
415// If we assume that final state hadrons after the charge exchange will be
416// in the ground states, we have to put ----------------------------------
417
418 if(M0projectile <
419 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass())
420 {
421 M0projectile=
422 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
423 M0projectile2 = M0projectile * M0projectile;
424 }
425
426 if(M0target <
427 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
428 {
429 M0target=
430 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
431 M0target2 = M0target * M0target;
432 }
433
434 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
435 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
436 /4./S;
437
438 if(PZcms2 < 0) {return false;} // It can be if energy is not sufficient for Delta
439//----------------------------------------------------------
440 projectile->SetDefinition(
441 G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode));
442
443 target->SetDefinition(
444 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));
445//----------------------------------------------------------
446 PZcms = std::sqrt(PZcms2);
447
448 Pprojectile.setPz( PZcms);
449 Pprojectile.setE(std::sqrt(M0projectile2+PZcms2));
450
451 Ptarget.setPz( -PZcms);
452 Ptarget.setE(std::sqrt(M0target2+PZcms2));
453
454 {
455 Pprojectile.transform(toLab);
456 Ptarget.transform(toLab);
457
458 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
459 projectile->SetPosition(target->GetPosition());
460
461 projectile->Set4Momentum(Pprojectile);
462 target->Set4Momentum(Ptarget);
463
464 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
465
466 return Result;
467 }
468 } // End of charge exchange part ------------------------------
469
470// ------------------------------------------------------------------
471 if(ProbOfDiffraction!=0.)
472 {
473 ProbProjectileDiffraction/=ProbOfDiffraction;
474 }
475 else
476 {
477 ProbProjectileDiffraction=0.;
478 }
479
480 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass *
481 ProjectileDiffStateMinMass;
482 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
483 ProjectileNonDiffStateMinMass;
484
485 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass *
486 TargetDiffStateMinMass;
487 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass *
488 TargetNonDiffStateMinMass;
489
490 G4double Pt2;
491 G4double ProjMassT2, ProjMassT;
492 G4double TargMassT2, TargMassT;
493 G4double PMinusMin, PMinusMax;
494// G4double PPlusMin , PPlusMax;
495 G4double TPlusMin , TPlusMax;
496 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
497
498 G4LorentzVector Qmomentum;
499 G4double Qminus, Qplus;
500
501 G4int whilecount=0;
502// Choose a process ---------------------------
503
504 if(G4UniformRand() < ProbOfDiffraction)
505 {
506 if(G4UniformRand() < ProbProjectileDiffraction)
507 { //-------- projectile diffraction ---------------
508 do {
509// Generate pt
510// if (whilecount++ >= 500 && (whilecount%100)==0)
511// G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
512// << ", loop count/ maxPtSquare : "
513// << whilecount << " / " << maxPtSquare << G4endl;
514 if (whilecount > 1000 )
515 {
516 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
517 target->SetStatus(2); return false; // Ignore this interaction
518 };
519// --------------- Check that the interaction is possible -----------
520 ProjMassT2=ProjectileDiffStateMinMass2;
521 ProjMassT =ProjectileDiffStateMinMass;
522
523 TargMassT2=M0target2;
524 TargMassT =M0target;
525
526 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
527 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
528 /4./S;
529
530 if(PZcms2 < 0 )
531 {
532 target->SetStatus(2);
533 return false;
534 }
535 maxPtSquare=PZcms2;
536
537 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
538 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
539
540 ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
541 ProjMassT =std::sqrt(ProjMassT2);
542
543 TargMassT2=M0target2+Pt2;
544 TargMassT =std::sqrt(TargMassT2);
545
546 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
547 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
548 /4./S;
549
550 if(PZcms2 < 0 ) continue;
551 PZcms =std::sqrt(PZcms2);
552
553 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
554 PMinusMax=SqrtS-TargMassT;
555
556 PMinusNew=ChooseP(PMinusMin, PMinusMax);
557// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
558
559 TMinusNew=SqrtS-PMinusNew;
560 Qminus=Ptarget.minus()-TMinusNew;
561 TPlusNew=TargMassT2/TMinusNew;
562 Qplus=Ptarget.plus()-TPlusNew;
563
564 Qmomentum.setPz( (Qplus-Qminus)/2 );
565 Qmomentum.setE( (Qplus+Qminus)/2 );
566 } while (
567((Pprojectile+Qmomentum).mag2() < ProjectileDiffStateMinMass2) || //No without excitation
568((Ptarget -Qmomentum).mag2() < M0target2 ));
569 }
570 else
571 { // -------------- Target diffraction ----------------
572 do {
573// Generate pt
574// if (whilecount++ >= 500 && (whilecount%100)==0)
575// G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
576// << ", loop count/ maxPtSquare : "
577// << whilecount << " / " << maxPtSquare << G4endl;
578 if (whilecount > 1000 )
579 {
580 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
581 target->SetStatus(2); return false; // Ignore this interaction
582 };
583// --------------- Check that the interaction is possible -----------
584 ProjMassT2=M0projectile2;
585 ProjMassT =M0projectile;
586
587 TargMassT2=TargetDiffStateMinMass2;
588 TargMassT =TargetDiffStateMinMass;
589
590 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
591 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
592 /4./S;
593
594 if(PZcms2 < 0 )
595 {
596 target->SetStatus(2);
597 return false;
598 }
599 maxPtSquare=PZcms2;
600
601 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
602 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
603
604 ProjMassT2=M0projectile2+Pt2;
605 ProjMassT =std::sqrt(ProjMassT2);
606
607 TargMassT2=TargetDiffStateMinMass2+Pt2;
608 TargMassT =std::sqrt(TargMassT2);
609
610 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
611 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
612 /4./S;
613
614 if(PZcms2 < 0 ) continue;
615 PZcms =std::sqrt(PZcms2);
616
617 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
618 TPlusMax=SqrtS-ProjMassT;
619
620 TPlusNew=ChooseP(TPlusMin, TPlusMax);
621
622//TPlusNew=TPlusMax;
623
624 PPlusNew=SqrtS-TPlusNew;
625 Qplus=PPlusNew-Pprojectile.plus();
626 PMinusNew=ProjMassT2/PPlusNew;
627 Qminus=PMinusNew-Pprojectile.minus();
628
629 Qmomentum.setPz( (Qplus-Qminus)/2 );
630 Qmomentum.setE( (Qplus+Qminus)/2 );
631
632 } while (
633 ((Pprojectile+Qmomentum).mag2() < M0projectile2 ) || //No without excitation
634 ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2));
635 }
636 }
637 else //----------- Non-diffraction process ------------
638 {
639 do {
640// Generate pt
641// if (whilecount++ >= 500 && (whilecount%100)==0)
642// G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
643// << ", loop count/ maxPtSquare : "
644// << whilecount << " / " << maxPtSquare << G4endl;
645 if (whilecount > 1000 )
646 {
647 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
648 target->SetStatus(2); return false; // Ignore this interaction
649 };
650// --------------- Check that the interaction is possible -----------
651 ProjMassT2=ProjectileNonDiffStateMinMass2;
652 ProjMassT =ProjectileNonDiffStateMinMass;
653
654 TargMassT2=TargetNonDiffStateMinMass2;
655 TargMassT =TargetNonDiffStateMinMass;
656
657 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
658 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
659 /4./S;
660
661 if(PZcms2 < 0 )
662 {
663 target->SetStatus(2);
664 return false;
665 }
666 maxPtSquare=PZcms2;
667 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
668 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
669
670 ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
671 ProjMassT =std::sqrt(ProjMassT2);
672
673 TargMassT2=TargetNonDiffStateMinMass2+Pt2;
674 TargMassT =std::sqrt(TargMassT2);
675
676 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
677 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
678 /4./S;
679
680 if(PZcms2 < 0 ) continue;
681 PZcms =std::sqrt(PZcms2);
682
683 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
684 PMinusMax=SqrtS-TargMassT;
685
686 PMinusNew=ChooseP(PMinusMin, PMinusMax);
687
688 Qminus=PMinusNew-Pprojectile.minus();
689
690 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
691// TPlusMax=SqrtS-PMinusNew;
692 TPlusMax=SqrtS-ProjMassT;
693
694 TPlusNew=ChooseP(TPlusMin, TPlusMax);
695
696 Qplus=-(TPlusNew-Ptarget.plus());
697
698 Qmomentum.setPz( (Qplus-Qminus)/2 );
699 Qmomentum.setE( (Qplus+Qminus)/2 );
700
701 } while (
702 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction
703 ((Ptarget -Qmomentum).mag2() < TargetNonDiffStateMinMass2 ));
704 }
705
706 Pprojectile += Qmomentum;
707 Ptarget -= Qmomentum;
708
709// Transform back and update SplitableHadron Participant.
710 Pprojectile.transform(toLab);
711 Ptarget.transform(toLab);
712
713// Calculation of the creation time ---------------------
714 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
715 projectile->SetPosition(target->GetPosition());
716// Creation time and position of target nucleon were determined at
717// ReggeonCascade() of G4FTFModel
718// ------------------------------------------------------
719
720 projectile->Set4Momentum(Pprojectile);
721 target->Set4Momentum(Ptarget);
722
723 projectile->IncrementCollisionCount(1);
724 target->IncrementCollisionCount(1);
725
726 return true;
727}
728
729// ---------------------------------------------------------------------
730void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron,
731 G4bool isProjectile,
732 G4ExcitedString * &FirstString,
733 G4ExcitedString * &SecondString,
734 G4FTFParameters *theParameters) const
735{
736 hadron->SplitUp();
737 G4Parton *start= hadron->GetNextParton();
738 if ( start==NULL)
739 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
740 FirstString=0; SecondString=0;
741 return;
742 }
743 G4Parton *end = hadron->GetNextParton();
744 if ( end==NULL)
745 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
746 FirstString=0; SecondString=0;
747 return;
748 }
749
750 G4LorentzVector Phadron=hadron->Get4Momentum();
751
752 G4LorentzVector Pstart(0.,0.,0.,0.);
753 G4LorentzVector Pend(0.,0.,0.,0.);
754 G4LorentzVector Pkink(0.,0.,0.,0.);
755 G4LorentzVector PkinkQ1(0.,0.,0.,0.);
756 G4LorentzVector PkinkQ2(0.,0.,0.,0.);
757
758 G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding());
759 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding());
760
761//--------------------------------------------------------------------------------
762 G4double Wmin(0.);
763 if(isProjectile)
764 {
765 Wmin=theParameters->GetProjMinDiffMass();
766 } else
767 {
768 Wmin=theParameters->GetTarMinDiffMass();
769 } // end of if(isProjectile)
770
771 G4double W = hadron->Get4Momentum().mag();
772 G4double W2=W*W;
773
774 G4double Pt(0.), x1(0.), x2(0.), x3(0.);
775 G4bool Kink=false;
776
777 if(W > Wmin)
778 { // Kink is possible
779 G4double Pt2kink=theParameters->GetPt2Kink();
780 Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
781
782 if(Pt > 500.*MeV)
783 {
784 G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
785 G4double Y=Ymax*(1.- 2.*G4UniformRand());
786
787 x1=1.-Pt/W*std::exp( Y);
788 x3=1.-Pt/W*std::exp(-Y);
789 x2=2.-x1-x3;
790
791 G4double Mass_startQ = 650.*MeV;
792 if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV;
793 if(PDGcode_startQ == 3) Mass_startQ = 500.*MeV;
794 if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV;
795
796 G4double Mass_endQ = 650.*MeV;
797 if(PDGcode_endQ < 3) Mass_endQ = 325.*MeV;
798 if(PDGcode_endQ == 3) Mass_endQ = 500.*MeV;
799 if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV;
800
801 G4double P2_1=W2*x1*x1/4.-Mass_endQ *Mass_endQ;
802 G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ;
803
804 G4double P2_2=sqr((2.-x1-x3)*W/2.);
805
806 if((P2_1 <= 0.) || (P2_3 <= 0.))
807 { Kink=false;}
808 else
809 {
810 G4double P_1=std::sqrt(P2_1);
811 G4double P_2=std::sqrt(P2_2);
812 G4double P_3=std::sqrt(P2_3);
813
814 G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
815 G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
816// Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09
817
818 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.))
819 { Kink=false;}
820 else
821 {
822 Kink=true;
823 Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09
824 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13);
825 Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1);
826 Pkink.setPx( Pt); Pkink.setPy( 0.); Pkink.setPz( P_2*CosT12);
827 Pstart.setE(x3*W/2.);
828 Pkink.setE(Pkink.vect().mag());
829 Pend.setE(x1*W/2.);
830
831 G4double XkQ=GetQuarkFractionOfKink(0.,1.);
832 if(Pkink.getZ() > 0.)
833 {
834 if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;}
835 } else {
836 if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;}
837 }
838
839 PkinkQ2=Pkink - PkinkQ1;
840//------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------
841
842 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
843 std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
844 G4double Psi=std::acos(Cos2Psi);
845
846G4LorentzRotation Rotate;
847if(isProjectile) {Rotate.rotateY(Psi);}
848else {Rotate.rotateY(pi-Psi);}
849Rotate.rotateZ(twopi*G4UniformRand());
850
851Pstart*=Rotate;
852Pkink*=Rotate;
853PkinkQ1*=Rotate;
854PkinkQ2*=Rotate;
855Pend*=Rotate;
856
857 }
858 } // end of if((P2_1 < 0.) || (P2_3 <0.))
859 } // end of if(Pt > 500.*MeV)
860 } // end of if(W > Wmin) Check for a kink
861
862//--------------------------------------------------------------------------------
863
864 if(Kink)
865 { // Kink is possible
866 std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp =
867 theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
868
869 G4int QuarkInGluon(1); G4double Ksi=G4UniformRand();
870 for(unsigned int Iq=0; Iq <3; Iq++)
871 {
872
873if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
874
875 G4Parton * Gquark = new G4Parton(QuarkInGluon);
876 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
877
878//-------------------------------------------------------------------------------
879 G4LorentzRotation toCMS(-1*Phadron.boostVector());
880
881 G4LorentzRotation toLab(toCMS.inverse());
882
883 Pstart.transform(toLab); start->Set4Momentum(Pstart);
884 PkinkQ1.transform(toLab);
885 PkinkQ2.transform(toLab);
886 Pend.transform(toLab); end->Set4Momentum(Pend);
887
888 G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding());
889
890 if(absPDGcode < 1000)
891 { // meson
892 if ( isProjectile )
893 { // Projectile
894 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end
895 { // Quark on the end
896 FirstString = new G4ExcitedString(end ,Ganti_quark, +1);
897 SecondString= new G4ExcitedString(Gquark,start ,+1);
898 Ganti_quark->Set4Momentum(PkinkQ1);
899 Gquark->Set4Momentum(PkinkQ2);
900
901 } else
902 { // Anti_Quark on the end
903 FirstString = new G4ExcitedString(end ,Gquark, +1);
904 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
905 Gquark->Set4Momentum(PkinkQ1);
906 Ganti_quark->Set4Momentum(PkinkQ2);
907
908 } // end of if(end->GetPDGcode() > 0)
909 } else { // Target
910 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end
911 { // Quark on the end
912 FirstString = new G4ExcitedString(Ganti_quark,end ,-1);
913 SecondString= new G4ExcitedString(start ,Gquark,-1);
914 Ganti_quark->Set4Momentum(PkinkQ2);
915 Gquark->Set4Momentum(PkinkQ1);
916
917 } else
918 { // Anti_Quark on the end
919 FirstString = new G4ExcitedString(Gquark,end ,-1);
920 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
921 Gquark->Set4Momentum(PkinkQ2);
922 Ganti_quark->Set4Momentum(PkinkQ1);
923
924 } // end of if(end->GetPDGcode() > 0)
925 } // end of if ( isProjectile )
926 } else // if(absPDGCode < 1000)
927 { // Baryon/AntiBaryon
928 if ( isProjectile )
929 { // Projectile
930 if((end->GetDefinition()->GetParticleType() == "diquarks") &&
931 (end->GetDefinition()->GetPDGEncoding() > 0 ) )
932 { // DiQuark on the end
933 FirstString = new G4ExcitedString(end ,Gquark, +1);
934 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
935 Gquark->Set4Momentum(PkinkQ1);
936 Ganti_quark->Set4Momentum(PkinkQ2);
937
938 } else
939 { // Anti_DiQuark on the end or quark
940 FirstString = new G4ExcitedString(end ,Ganti_quark, +1);
941 SecondString= new G4ExcitedString(Gquark,start ,+1);
942 Ganti_quark->Set4Momentum(PkinkQ1);
943 Gquark->Set4Momentum(PkinkQ2);
944
945 } // end of if(end->GetPDGcode() > 0)
946 } else { // Target
947
948 if((end->GetDefinition()->GetParticleType() == "diquarks") &&
949 (end->GetDefinition()->GetPDGEncoding() > 0 ) )
950 { // DiQuark on the end
951 FirstString = new G4ExcitedString(Gquark,end ,-1);
952
953 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
954 Gquark->Set4Momentum(PkinkQ1);
955 Ganti_quark->Set4Momentum(PkinkQ2);
956
957 } else
958 { // Anti_DiQuark on the end or Q
959 FirstString = new G4ExcitedString(Ganti_quark,end ,-1);
960 SecondString= new G4ExcitedString(start ,Gquark,-1);
961 Gquark->Set4Momentum(PkinkQ2);
962 Ganti_quark->Set4Momentum(PkinkQ1);
963
964 } // end of if(end->GetPDGcode() > 0)
965 } // end of if ( isProjectile )
966 } // end of if(absPDGcode < 1000)
967
968 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
969 FirstString->SetPosition(hadron->GetPosition());
970
971 SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation());
972 SecondString->SetPosition(hadron->GetPosition());
973
974// -------------------------------------------------------------------------
975 } else // End of kink is possible
976 { // Kink is impossible
977 if ( isProjectile )
978 {
979 FirstString= new G4ExcitedString(end,start, +1);
980 } else {
981 FirstString= new G4ExcitedString(start,end, -1);
982 }
983
984 SecondString=0;
985
986 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
987 FirstString->SetPosition(hadron->GetPosition());
988
989// momenta of string ends
990//
991 G4double Momentum=hadron->Get4Momentum().vect().mag();
992 G4double Plus=hadron->Get4Momentum().e() + Momentum;
993 G4double Minus=hadron->Get4Momentum().e() - Momentum;
994
995 G4ThreeVector tmp;
996 if(Momentum > 0.)
997 {
998 tmp.set(hadron->Get4Momentum().px(),
999 hadron->Get4Momentum().py(),
1000 hadron->Get4Momentum().pz());
1001 tmp/=Momentum;
1002 }
1003 else
1004 {
1005 tmp.set(0.,0.,1.);
1006 }
1007
1008 G4LorentzVector Pstart(tmp,0.);
1009 G4LorentzVector Pend(tmp,0.);
1010
1011 if(isProjectile)
1012 {
1013 Pstart*=(-1.)*Minus/2.;
1014 Pend *=(+1.)*Plus /2.;
1015 }
1016 else
1017 {
1018 Pstart*=(+1.)*Plus/2.;
1019 Pend *=(-1.)*Minus/2.;
1020 }
1021
1022 Momentum=-Pstart.mag();
1023 Pstart.setT(Momentum); // It is assumed that quark has m=0.
1024
1025 Momentum=-Pend.mag();
1026 Pend.setT(Momentum); // It is assumed that di-quark has m=0.
1027
1028 start->Set4Momentum(Pstart);
1029 end->Set4Momentum(Pend);
1030 SecondString=0;
1031 } // End of kink is impossible
1032
1033#ifdef G4_FTFDEBUG
1034 G4cout << " generated string flavors "
1035 << start->GetPDGcode() << " / "
1036 << end->GetPDGcode() << G4endl;
1037 G4cout << " generated string momenta: quark "
1038 << start->Get4Momentum() << "mass : "
1039 <<start->Get4Momentum().mag() << G4endl;
1040 G4cout << " generated string momenta: Diquark "
1041 << end ->Get4Momentum()
1042 << "mass : " <<end->Get4Momentum().mag()<< G4endl;
1043 G4cout << " sum of ends " << Pstart+Pend << G4endl;
1044 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
1045#endif
1046
1047 return;
1048
1049}
1050
1051
1052// --------- private methods ----------------------
1053
1054// ---------------------------------------------------------------------
1055G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
1056{
1057// choose an x between Xmin and Xmax with P(x) ~ 1/x
1058// to be improved...
1059
1060 G4double range=Pmax-Pmin;
1061
1062 if ( Pmin <= 0. || range <=0. )
1063 {
1064 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1065 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
1066 }
1067
1068 G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());
1069 return P;
1070}
1071
1072// ---------------------------------------------------------------------
1073G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2,
1074 G4double maxPtSquare) const
1075{ // @@ this method is used in FTFModel as well. Should go somewhere common!
1076
1077 G4double Pt2(0.);
1078 if(AveragePt2 <= 0.) {Pt2=0.;}
1079 else
1080 {
1081 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
1082 (std::exp(-maxPtSquare/AveragePt2)-1.));
1083 }
1084 G4double Pt=std::sqrt(Pt2);
1085 G4double phi=G4UniformRand() * twopi;
1086 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
1087}
1088
1089// ---------------------------------------------------------------------
1090G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
1091 {
1092 G4double z, yf;
1093 do {
1094 z = zmin + G4UniformRand()*(zmax-zmin);
1095 yf = z*z +sqr(1 - z);
1096 }
1097 while (G4UniformRand() > yf);
1098 return z;
1099 }
1100// ---------------------------------------------------------------------
1101void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09
1102 {
1103 G4int absIdPDG = std::abs(IdPDG);
1104 Q1= absIdPDG/ 100;
1105 Q2= (absIdPDG %100)/10;
1106
1107 G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 );
1108
1109 if (IdPDG < 0 ) anti *=-1;
1110 Q1 *= anti;
1111 Q2 *= -1 * anti;
1112 return;
1113 }
1114// ---------------------------------------------------------------------
1115void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG,
1116 G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09
1117 {
1118 Q1 = IdPDG / 1000;
1119 Q2 = (IdPDG % 1000) / 100;
1120 Q3 = (IdPDG % 100) / 10;
1121 return;
1122 }
1123// ---------------------------------------------------------------------
1124G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09
1125 {
1126 G4int TmpQ(0);
1127 if( Q3 > Q2 )
1128 {
1129 TmpQ = Q2;
1130 Q2 = Q3;
1131 Q3 = TmpQ;
1132 } else if( Q3 > Q1 )
1133 {
1134 TmpQ = Q1;
1135 Q1 = Q3;
1136 Q3 = TmpQ;
1137 }
1138
1139 if( Q2 > Q1 )
1140 {
1141 TmpQ = Q1;
1142 Q1 = Q2;
1143 Q2 = TmpQ;
1144 }
1145
1146 G4int NewCode = Q1*1000 + Q2* 100 + Q3* 10 + 2;
1147 return NewCode;
1148 }
1149// ---------------------------------------------------------------------
1150G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
1151{
1152 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
1153}
1154
1155
1156G4DiffractiveExcitation::~G4DiffractiveExcitation()
1157{
1158}
1159
1160
1161const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
1162{
1163 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
1164 return *this;
1165}
1166
1167
1168int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
1169{
1170 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
1171 return false;
1172}
1173
1174int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
1175{
1176 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
1177 return true;
1178}
Note: See TracBrowser for help on using the repository browser.