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

Last change on this file since 1337 was 1228, checked in by garnier, 15 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.