| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // << 27 // $Id$
28 <<
29 // ------------------------------------------- 28 // ------------------------------------------------------------
30 // GEANT 4 class implemetation file 29 // GEANT 4 class implemetation file
31 // 30 //
32 // ---------------- G4DiffractiveExcitati 31 // ---------------- G4DiffractiveExcitation --------------
33 // by Gunter Folger, October 1998. 32 // by Gunter Folger, October 1998.
34 // diffractive Excitation used by strin 33 // diffractive Excitation used by strings models
35 // Take a projectile and a target 34 // Take a projectile and a target
36 // excite the projectile and targe 35 // excite the projectile and target
37 // Essential changed by V. Uzhinsky in Novemb 36 // Essential changed by V. Uzhinsky in November - December 2006
38 // in order to put it in a correspondence wit 37 // in order to put it in a correspondence with original FRITIOF
39 // model. Variant of FRITIOF with nucleon de- 38 // model. Variant of FRITIOF with nucleon de-excitation is implemented.
40 // Other changes by V.Uzhinsky in May 2007 we 39 // Other changes by V.Uzhinsky in May 2007 were introduced to fit
41 // meson-nucleon interactions. Additional cha 40 // meson-nucleon interactions. Additional changes by V. Uzhinsky
42 // were introduced in December 2006. They tre 41 // were introduced in December 2006. They treat diffraction dissociation
43 // processes more exactly. 42 // processes more exactly.
44 // Correct treatment of the diffraction disso <<
45 // Mass distributions for resonances and uu-d <<
46 // and dd-diquarks suppression in neutrons we <<
47 // ------------------------------------------- 43 // ---------------------------------------------------------------------
48 44
>> 45
49 #include "globals.hh" 46 #include "globals.hh"
50 #include "Randomize.hh" 47 #include "Randomize.hh"
51 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
53 50
54 #include "G4DiffractiveExcitation.hh" 51 #include "G4DiffractiveExcitation.hh"
55 #include "G4FTFParameters.hh" 52 #include "G4FTFParameters.hh"
56 #include "G4ElasticHNScattering.hh" 53 #include "G4ElasticHNScattering.hh"
57 54
>> 55 #include "G4LorentzRotation.hh"
58 #include "G4RotationMatrix.hh" 56 #include "G4RotationMatrix.hh"
>> 57 #include "G4ThreeVector.hh"
59 #include "G4ParticleDefinition.hh" 58 #include "G4ParticleDefinition.hh"
60 #include "G4ParticleTable.hh" <<
61 #include "G4SampleResonance.hh" <<
62 #include "G4VSplitableHadron.hh" 59 #include "G4VSplitableHadron.hh"
63 #include "G4ExcitedString.hh" 60 #include "G4ExcitedString.hh"
>> 61 #include "G4ParticleTable.hh"
64 #include "G4Neutron.hh" 62 #include "G4Neutron.hh"
65 << 63 #include "G4ParticleDefinition.hh"
66 #include "G4Exp.hh" <<
67 #include "G4Log.hh" <<
68 #include "G4Pow.hh" <<
69 64
70 //#include "G4ios.hh" 65 //#include "G4ios.hh"
>> 66 //#include "UZHI_diffraction.hh"
71 67
72 //============================================ << 68 G4DiffractiveExcitation::G4DiffractiveExcitation()
73 << 69 {
74 //#define debugFTFexcitation <<
75 //#define debug_heavyHadrons <<
76 <<
77 //============================================ <<
78 <<
79 G4DiffractiveExcitation::G4DiffractiveExcitati <<
80 <<
81 <<
82 //============================================ <<
83 <<
84 G4DiffractiveExcitation::~G4DiffractiveExcitat <<
85 <<
86 <<
87 //============================================ <<
88 <<
89 G4bool G4DiffractiveExcitation::ExciteParticip <<
90 <<
91 <<
92 <<
93 <<
94 #ifdef debugFTFexcitation <<
95 G4cout << G4endl << "FTF ExciteParticipants <<
96 #endif <<
97 <<
98 CommonVariables common; <<
99 <<
100 // Projectile parameters <<
101 common.Pprojectile = projectile->Get4Momentu <<
102 if ( common.Pprojectile.z() < 0.0 ) return f <<
103 common.ProjectilePDGcode = projectile->GetDe <<
104 common.absProjectilePDGcode = std::abs( comm <<
105 common.M0projectile = projectile->GetDefinit <<
106 G4double ProjectileRapidity = common.Pprojec <<
107 <<
108 // Target parameters <<
109 common.Ptarget = target->Get4Momentum(); <<
110 common.TargetPDGcode = target->GetDefinition <<
111 common.absTargetPDGcode = std::abs( common.T <<
112 common.M0target = target->GetDefinition()->G <<
113 G4double TargetRapidity = common.Ptarget.rap <<
114 <<
115 // Kinematical properties of the interaction <<
116 G4LorentzVector Psum = common.Pprojectile + <<
117 common.S = Psum.mag2(); <<
118 common.SqrtS = std::sqrt( common.S ); <<
119 <<
120 // Check off-shellness of the participants <<
121 G4bool toBePutOnMassShell = true; //Uzhi Au <<
122 common.MminProjectile = common.BrW.GetMinimu <<
123 /* Uzhi Aug.2019 <<
124 if ( common.M0projectile < common.MminProjec <<
125 toBePutOnMassShell = true; <<
126 common.M0projectile = common.BrW.SampleMas <<
127 <<
128 <<
129 } <<
130 */ <<
131 common.M0projectile2 = common.M0projectile * <<
132 common.ProjectileDiffStateMinMass = thePa <<
133 common.ProjectileNonDiffStateMinMass = thePa <<
134 if ( common.M0projectile > common.Projectile <<
135 common.ProjectileDiffStateMinMass = com <<
136 common.ProjectileNonDiffStateMinMass = com <<
137 if ( common.absProjectilePDGcode > 3000 ) <<
138 common.ProjectileDiffStateMinMass += <<
139 common.ProjectileNonDiffStateMinMass += <<
140 } <<
141 } <<
142 common.MminTarget = common.BrW.GetMinimumMas <<
143 /* Uzhi Aug.2019 <<
144 if ( common.M0target < common.MminTarget ) { <<
145 toBePutOnMassShell = true; <<
146 common.M0target = common.BrW.SampleMass( t <<
147 t <<
148 + <<
149 } <<
150 */ <<
151 common.M0target2 = common.M0target * common. <<
152 common.TargetDiffStateMinMass = theParame <<
153 common.TargetNonDiffStateMinMass = theParame <<
154 if ( common.M0target > common.TargetDiffStat <<
155 common.TargetDiffStateMinMass = common. <<
156 common.TargetNonDiffStateMinMass = common. <<
157 if ( common.absTargetPDGcode > 3000 ) { / <<
158 common.TargetDiffStateMinMass += 140. <<
159 common.TargetNonDiffStateMinMass += 140. <<
160 } <<
161 }; <<
162 #ifdef debugFTFexcitation <<
163 G4cout << "Proj Targ PDGcodes " << common.Pr <<
164 << "Mprojectile Y " << common.Pproje <<
165 << "M0projectile Y " << common.M0proj <<
166 G4cout << "Mtarget Y " << common.Ptarge <<
167 << "M0target Y " << common.M0targ <<
168 G4cout << "Pproj " << common.Pprojectile << <<
169 #endif <<
170 <<
171 // Transform momenta to cms and then rotate <<
172 common.toCms = G4LorentzRotation( -1 * Psum. <<
173 G4LorentzVector Ptmp = common.toCms * common <<
174 if ( Ptmp.pz() <= 0.0 ) return false; // "S <<
175 common.toCms.rotateZ( -1*Ptmp.phi() ); <<
176 common.toCms.rotateY( -1*Ptmp.theta() ); <<
177 common.toLab = common.toCms.inverse(); <<
178 common.Pprojectile.transform( common.toCms ) <<
179 common.Ptarget.transform( common.toCms ); <<
180 <<
181 G4double SumMasses = common.M0projectile + c <<
182 #ifdef debugFTFexcitation <<
183 G4cout << "SqrtS " << common.SqrtS << G4 <<
184 << " " << common.M0target << " " << S <<
185 #endif <<
186 if ( common.SqrtS < SumMasses ) return false <<
187 <<
188 common.PZcms2 = ( sqr( common.S ) + sqr( com <<
189 - 2.0 * ( common.S * ( com <<
190 + common.M0proje <<
191 #ifdef debugFTFexcitation <<
192 G4cout << "PZcms2 after toBePutOnMassShell " <<
193 #endif <<
194 if ( common.PZcms2 < 0.0 ) return false; // <<
195 <<
196 common.PZcms = std::sqrt( common.PZcms2 ); <<
197 if ( toBePutOnMassShell ) { <<
198 if ( common.Pprojectile.z() > 0.0 ) { <<
199 common.Pprojectile.setPz( common.PZcms <<
200 common.Ptarget.setPz( -common.PZcms <<
201 } else { <<
202 common.Pprojectile.setPz( -common.PZcms <<
203 common.Ptarget.setPz( common.PZcms <<
204 } <<
205 common.Pprojectile.setE( std::sqrt( common <<
206 + comm <<
207 + comm <<
208 + comm <<
209 common.Ptarget.setE( std::sqrt( common.M0t <<
210 + common.P <<
211 + common.P <<
212 + common.P <<
213 } <<
214 #ifdef debugFTFexcitation <<
215 G4cout << "Start --------------------" << G4 <<
216 << " " << common.ProjectileDiffStateM <<
217 << G4endl <<
218 << "Targ M0 Mdif Mndif " << common.M0 <<
219 << " " << common.TargetNonDiffStateMi <<
220 << "Proj CMS " << common.Pprojectile <<
221 #endif <<
222 <<
223 // Check for possible quark exchange <<
224 ProjectileRapidity = common.Pprojectile.rapi <<
225 TargetRapidity = common.Ptarget.rapidity(); <<
226 G4double QeNoExc = theParameters->GetProcPro <<
227 G4double QeExc = theParameters->GetProcPro <<
228 theParameters->GetProcPro <<
229 common.ProbProjectileDiffraction = <<
230 theParameters->GetProcProb( 2, ProjectileR <<
231 common.ProbTargetDiffraction = <<
232 theParameters->GetProcProb( 3, ProjectileR <<
233 common.ProbOfDiffraction = common.ProbProjec <<
234 #ifdef debugFTFexcitation <<
235 G4cout << "Proc Probs " << QeNoExc << " " << <<
236 << common.ProbProjectileDiffraction < <<
237 << "ProjectileRapidity " << Projectil <<
238 #endif <<
239 <<
240 if ( QeNoExc + QeExc + common.ProbProjectile <<
241 QeNoExc = 1.0 - QeExc - common.ProbProject <<
242 } <<
243 if ( QeExc + QeNoExc != 0.0 ) { <<
244 common.ProbExc = QeExc / ( QeExc + QeNoExc <<
245 } <<
246 if ( 1.0 - QeExc - QeNoExc > 0.0 ) { <<
247 common.ProbProjectileDiffraction /= ( 1.0 <<
248 common.ProbTargetDiffraction /= ( 1.0 <<
249 } <<
250 #ifdef debugFTFexcitation <<
251 G4cout << "Proc Probs " << QeNoExc << " " << <<
252 << common.ProbProjectileDiffraction < <<
253 << "ProjectileRapidity " << Projectil <<
254 #endif <<
255 <<
256 // Try out charge exchange <<
257 G4int returnCode = 1; <<
258 if ( G4UniformRand() < QeExc + QeNoExc ) { <<
259 returnCode = ExciteParticipants_doChargeEx <<
260 <<
261 } <<
262 <<
263 G4bool returnResult = false; <<
264 if ( returnCode == 0 ) { <<
265 returnResult = true; // Successfully ende <<
266 } else if ( returnCode == 1 ) { <<
267 <<
268 common.ProbOfDiffraction = common.ProbProj <<
269 #ifdef debugFTFexcitation <<
270 G4cout << "Excitation -------------------- <<
271 << "Proj M0 MdMin MndMin " << commo <<
272 << common.ProjectileDiffStateMinMas <<
273 << G4endl <<
274 << "Targ M0 MdMin MndMin " << commo <<
275 << " " << common.TargetNonDiffState <<
276 << G4endl <<
277 << "Prob: ProjDiff TargDiff + Sum " <<
278 << common.ProbTargetDiffraction << <<
279 #endif <<
280 if ( common.ProbOfDiffraction != 0.0 ) { <<
281 common.ProbProjectileDiffraction /= comm <<
282 } else { <<
283 common.ProbProjectileDiffraction = 0.0; <<
284 } <<
285 #ifdef debugFTFexcitation <<
286 G4cout << "Prob: ProjDiff TargDiff + Sum " <<
287 << common.ProbTargetDiffraction << <<
288 #endif <<
289 common.ProjectileDiffStateMinMass2 = sq <<
290 common.ProjectileNonDiffStateMinMass2 = sq <<
291 common.TargetDiffStateMinMass2 = sq <<
292 common.TargetNonDiffStateMinMass2 = sq <<
293 // Choose between diffraction and non-diff <<
294 if ( G4UniformRand() < common.ProbOfDiffra <<
295 <<
296 returnResult = ExciteParticipants_doDiff <<
297 <<
298 } else { <<
299 <<
300 returnResult = ExciteParticipants_doNonD <<
301 <<
302 } <<
303 if ( returnResult ) { <<
304 common.Pprojectile += common.Qmomentum; <<
305 common.Ptarget -= common.Qmomentum; <<
306 // Transform back and update SplitableHa <<
307 common.Pprojectile.transform( common.toL <<
308 common.Ptarget.transform( common.toLab ) <<
309 #ifdef debugFTFexcitation <<
310 G4cout << "Mproj " << common.Pprojectile <<
311 << G4endl; <<
312 #endif <<
313 projectile->Set4Momentum( common.Pprojec <<
314 target->Set4Momentum( common.Ptarget ); <<
315 projectile->IncrementCollisionCount( 1 ) <<
316 target->IncrementCollisionCount( 1 ); <<
317 } <<
318 } <<
319 <<
320 return returnResult; <<
321 } 70 }
322 71
323 //-------------------------------------------- << 72 // ---------------------------------------------------------------------
>> 73 G4bool G4DiffractiveExcitation::
>> 74 ExciteParticipants(G4VSplitableHadron *projectile,
>> 75 G4VSplitableHadron *target,
>> 76 G4FTFParameters *theParameters,
>> 77 G4ElasticHNScattering *theElastic) const
>> 78 {
>> 79 //G4cout<<G4endl<<"ExciteParticipants --------------"<<G4endl;
>> 80 // -------------------- Projectile parameters -----------------------
>> 81 G4LorentzVector Pprojectile=projectile->Get4Momentum();
>> 82
>> 83 if(Pprojectile.z() < 0.)
>> 84 {
>> 85 target->SetStatus(2);
>> 86 return false;
>> 87 }
>> 88
>> 89 G4double ProjectileRapidity = Pprojectile.rapidity();
>> 90
>> 91 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
>> 92 G4int absProjectilePDGcode=std::abs(ProjectilePDGcode);
>> 93
>> 94 G4bool PutOnMassShell(false);
>> 95 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
>> 96 G4double M0projectile = Pprojectile.mag(); // Without de-excitation
>> 97 //G4cout<<"M0projectile "<<M0projectile<<" "<<ProjectileRapidity<<G4endl;
>> 98
>> 99 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
>> 100 {
>> 101 PutOnMassShell=true;
>> 102 M0projectile=projectile->GetDefinition()->GetPDGMass();
>> 103 }
>> 104
>> 105 G4double M0projectile2 = M0projectile * M0projectile;
>> 106
>> 107 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
>> 108 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
>> 109 G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
>> 110
>> 111 // -------------------- Target parameters -------------------------
>> 112 G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
>> 113 G4int absTargetPDGcode=std::abs(TargetPDGcode);
>> 114 //G4cout<<"Entry to QE or Excit "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;
>> 115
>> 116 G4LorentzVector Ptarget=target->Get4Momentum();
>> 117 //G4cout<<"Pproj "<<Pprojectile<<G4endl;
>> 118 //G4cout<<"Ptarget "<<Ptarget<<G4endl;
>> 119 G4double M0target = Ptarget.mag();
>> 120
>> 121 // G4double TargetRapidity = Ptarget.rapidity();
>> 122 //G4cout<<"M0target "<<M0target<<" "<<TargetRapidity<<G4endl;
>> 123 if(M0target < target->GetDefinition()->GetPDGMass())
>> 124 {
>> 125 PutOnMassShell=true;
>> 126 M0target=target->GetDefinition()->GetPDGMass();
>> 127 }
324 128
325 G4int G4DiffractiveExcitation:: << 129 G4double M0target2 = M0target * M0target;
326 ExciteParticipants_doChargeExchange( G4VSplita << 130
327 G4VSplita << 131 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();
328 G4FTFPara << 132 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();
329 G4Elastic << 133 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
330 G4Diffrac << 134
331 // First of the three utility methods used o << 135 G4double AveragePt2=theParameters->GetAveragePt2();
332 // it does the sampling for the charge-excha << 136 G4double ProbLogDistr=theParameters->GetProbLogDistr(); // Uzhi 21.05.2012
333 // This method returns an integer code - ins << 137
334 // "0" : successfully ended and nothing el << 138 // G4double ProbOfDiffraction=ProbProjectileDiffraction +
335 // "1" : successfully completed, but the w << 139 // ProbTargetDiffraction;
336 // "99" : unsuccessfully ended, nothing els << 140
337 G4int returnCode = 99; << 141 // G4double SumMasses=M0projectile+M0target+220.*MeV; // 200->220 7 June 2011
338 << 142 G4double SumMasses=M0projectile+M0target+220.*MeV; // 200->220 7 June 2011
339 G4double DeltaProbAtQuarkExchange = theParam << 143
340 G4ParticleDefinition* TestParticle = 0; << 144 // Kinematical properties of the interactions --------------
341 G4double MtestPr = 0.0, MtestTr = 0.0; << 145 G4LorentzVector Psum; // 4-momentum in CMS
342 << 146 Psum=Pprojectile+Ptarget;
343 #ifdef debugFTFexcitation << 147 G4double S=Psum.mag2();
344 G4cout << "Q exchange ---------------------- << 148
345 #endif << 149 //Uzhi_SqrtS=std::sqrt(S);
346 << 150
347 // The target hadron is always a nucleon (i. << 151 // Transform momenta to cms and then rotate parallel to z axis;
348 // never an antinucleon), therefore only a q << 152 G4LorentzRotation toCms(-1*Psum.boostVector());
349 // exchanged between the projectile hadron a << 153
350 // we could get a quark-quark-antiquark syst << 154 G4LorentzVector Ptmp=toCms*Pprojectile;
351 // This implies that any projectile meson or << 155
352 // a constituent quark in all cases - can ha << 156 if ( Ptmp.pz() <= 0. )
353 // hadron. Instead, any projectile anti-bary << 157 {
354 // with a target hadron (because it has only << 158 target->SetStatus(2);
355 // projectile baryons, instead can have char << 159 // "String" moving backwards in CMS, abort collision !!
356 << 160 return false;
357 G4int NewProjCode = 0, NewTargCode = 0, Proj << 161 }
358 // Projectile unpacking << 162
359 if ( common.absProjectilePDGcode < 1000 ) { << 163 toCms.rotateZ(-1*Ptmp.phi());
360 UnpackMeson( common.ProjectilePDGcode, Pr << 164 toCms.rotateY(-1*Ptmp.theta());
361 } else { << 165
362 UnpackBaryon( common.ProjectilePDGcode, Pr << 166 G4LorentzRotation toLab(toCms.inverse());
363 } << 167
364 // Target unpacking << 168 Pprojectile.transform(toCms);
365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0; << 169 Ptarget.transform(toCms);
366 UnpackBaryon( common.TargetPDGcode, TargQ1, << 170
367 #ifdef debugFTFexcitation << 171 G4double PZcms2, PZcms;
368 G4cout << "Proj Quarks " << ProjQ1 << " " << << 172
369 << "Targ Quarks " << TargQ1 << " " << << 173 G4double SqrtS=std::sqrt(S);
370 #endif << 174
371 << 175 /*
372 // Sampling of exchanged quarks << 176 G4cout<<"Proj "<<Pprojectile<<G4endl;
373 G4int ProjExchangeQ = 0, TargExchangeQ = 0; << 177 G4cout<<"Targ "<<Ptarget<<G4endl;
374 if ( common.absProjectilePDGcode < 1000 ) { << 178 G4cout<<"SqrtS "<<SqrtS<<G4endl;
375 << 179 G4cout<<"M0pr M0tr "<<M0projectile<<" "<<M0target<<" "<<SumMasses<<G4endl;
376 G4bool isProjQ1Quark = false; << 180 */
377 ProjExchangeQ = ProjQ2; << 181 //G4cout<<"SqrtS < 2300*MeV Bary "<<SqrtS<<G4endl;
378 if ( ProjQ1 > 0 ) { // ProjQ1 is a quark << 182 // if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses)) // 7.06.11
379 isProjQ1Quark = true; << 183 if(absProjectilePDGcode > 1000 && (SqrtS < SumMasses))
380 ProjExchangeQ = ProjQ1; << 184 {target->SetStatus(2); return false;} // The model cannot work for
381 } << 185 // p+p-interactions
382 // Exchange of non-identical quarks is all << 186 // at Plab < 1.62 GeV/c.
383 G4int NpossibleStates = 0; << 187
384 if ( ProjExchangeQ != TargQ1 ) NpossibleSt << 188 //G4cout<<"SqrtS < 1600*MeV Pion "<<SqrtS<<G4endl;
385 if ( ProjExchangeQ != TargQ2 ) NpossibleSt << 189
386 if ( ProjExchangeQ != TargQ3 ) NpossibleSt << 190 // if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) && // 7.06.11
387 G4int Nsampled = (G4int)G4RandFlat::shootI << 191 // ((SqrtS < 1600*MeV) || (SqrtS < SumMasses)))
388 NpossibleStates = 0; << 192 if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) &&
389 if ( ProjExchangeQ != TargQ1 ) { << 193 (SqrtS < SumMasses))
390 if ( ++NpossibleStates == Nsampled ) { << 194 {target->SetStatus(2); return false;} // The model cannot work for
391 TargExchangeQ = TargQ1; TargQ1 = ProjE << 195 // Pi+p-interactions
392 isProjQ1Quark ? ProjQ1 = TargExchangeQ << 196 // at Plab < 1. GeV/c.
393 } << 197 //G4cout<<"SqrtS < 1600*MeV "<<SqrtS<<G4endl;
394 } << 198 //SumMasses=M0projectile+M0target+20.*MeV;
395 if ( ProjExchangeQ != TargQ2 ) { << 199 if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) ||
396 if ( ++NpossibleStates == Nsampled ) { << 200 (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
397 TargExchangeQ = TargQ2; TargQ2 = ProjE << 201 (absProjectilePDGcode == 310)) && (SqrtS < SumMasses))
398 isProjQ1Quark ? ProjQ1 = TargExchangeQ << 202 // (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses)))
399 } << 203 {target->SetStatus(2); return false;} // The model cannot work for
400 } << 204 // K+p-interactions
401 if ( ProjExchangeQ != TargQ3 ) { << 205 // at Plab < ??? GeV/c. ???
402 if ( ++NpossibleStates == Nsampled ) { << 206
403 TargExchangeQ = TargQ3; TargQ3 = ProjE << 207 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
404 isProjQ1Quark ? ProjQ1 = TargExchangeQ << 208 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
405 } << 209 /4./S;
406 } << 210 //G4cout<<"PZcms2 "<<PZcms2<<G4endl;
407 #ifdef debugFTFexcitation << 211 if(PZcms2 < 0)
408 G4cout << "Exchanged Qs in Pr Tr " << Proj << 212 {target->SetStatus(2); return false;} // It can be in an interaction with
409 #endif << 213 // off-shell nuclear nucleon
410 << 214
411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ << 215 PZcms = std::sqrt(PZcms2);
412 G4bool ProjExcited = false; << 216
413 const G4int maxNumberOfAttempts = 50; << 217 if(PutOnMassShell)
414 G4int attempts = 0; << 218 {
415 while ( attempts++ < maxNumberOfAttempts ) << 219 if(Pprojectile.z() > 0.)
416 << 220 {
417 // Determination of a new projectile ID << 221 Pprojectile.setPz( PZcms);
418 G4double ProbSpin0 = 0.5; << 222 Ptarget.setPz( -PZcms);
419 G4double Ksi = G4UniformRand(); << 223 }
420 if ( aProjQ1 == aProjQ2 ) { << 224 else
421 if ( G4UniformRand() < ProbSpin0 ) { << 225 {
422 if ( aProjQ1 < 3 ) { << 226 Pprojectile.setPz(-PZcms);
423 NewProjCode = 111; << 227 Ptarget.setPz( PZcms);
424 if ( Ksi < 0.5 ) { << 228 };
425 NewProjCode = 221; << 229
426 if ( Ksi < 0.25 ) { << 230 Pprojectile.setE(std::sqrt(M0projectile2 +
427 NewProjCode = 331; << 231 Pprojectile.x()*Pprojectile.x()+
428 } << 232 Pprojectile.y()*Pprojectile.y()+
429 } << 233 PZcms2));
430 } else if ( aProjQ1 == 3 ) { << 234 Ptarget.setE(std::sqrt(M0target2 +
431 NewProjCode = 221; << 235 Ptarget.x()*Ptarget.x()+
432 if ( Ksi < 0.5 ) { << 236 Ptarget.y()*Ptarget.y()+
433 NewProjCode = 331; << 237 PZcms2));
434 } << 238 }
435 } else if ( aProjQ1 == 4 ) { << 239
436 NewProjCode = 441; // eta << 240
437 } else if ( aProjQ1 == 5 ) { << 241 //G4cout<<"Proj "<<Pprojectile<<G4endl;
438 NewProjCode = 551; // eta << 242 //G4cout<<"Targ "<<Ptarget<<G4endl;
439 } << 243 G4double maxPtSquare; // = PZcms2;
440 } else { << 244 /*
441 if ( aProjQ1 < 3 ) { << 245 Uzhi_targetdiffraction=0;
442 NewProjCode = 113; << 246 Uzhi_projectilediffraction=0;
443 if ( Ksi < 0.5 ) { << 247 Uzhi_Mx2=1.0;
444 NewProjCode = 223; << 248 Uzhi_modT=0.;
445 } << 249 G4int Uzhi_QE=0;
446 } else if ( aProjQ1 == 3 ) { << 250 */
447 NewProjCode = 333; << 251 /*
448 } else if ( aProjQ1 == 4 ) { << 252 G4cout<<"Start --------------------"<<G4endl;
449 NewProjCode = 443; // J/p << 253 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<G4endl;
450 } else if ( aProjQ1 == 5 ) { << 254 G4cout<<"Targ "<<M0target <<" "<<TargetDiffStateMinMass <<" "<<TargetNonDiffStateMinMass<<G4endl;
451 NewProjCode = 553; // Ups << 255 G4cout<<"SqrtS "<<SqrtS<<G4endl;
452 } << 256 G4cout<<"Rapid "<<ProjectileRapidity<<G4endl; //" "<<TargetRapidity<<G4endl;
453 } << 257 */
454 } else { << 258 // Charge exchange can be possible for baryons -----------------
455 if ( aProjQ1 > aProjQ2 ) { << 259
456 NewProjCode = aProjQ1*100 + aProjQ2* << 260 // Getting the values needed for exchange ----------------------
457 } else { << 261 G4double MagQuarkExchange =theParameters->GetMagQuarkExchange();//0.045; //
458 NewProjCode = aProjQ2*100 + aProjQ1* << 262 G4double SlopeQuarkExchange =theParameters->GetSlopeQuarkExchange();//0.; //
>> 263 // 777777777777777777777777777777777777777777777777777777777777777777777777777777
>> 264 G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
>> 265
>> 266 // G4double NucleonMass=
>> 267 // (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();
>> 268 G4double DeltaMass=
>> 269 (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
>> 270
>> 271 //G4double TargetRapidity(0.);
>> 272 //G4cout<<"Prob Q Exch "<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*ProjectileRapidity)<<G4endl;
>> 273
>> 274 //G4cout<<"Q exc Mag Slop Wdelta"<<MagQuarkExchange<<" "<<SlopeQuarkExchange<<" "<<DeltaProbAtQuarkExchange<<G4endl;
>> 275 //G4cout<<"ProjectileRapidity "<<ProjectileRapidity<<G4endl;
>> 276 //G4cout<<"Prob Exc "<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity))<<G4endl;
>> 277 //G4int Uzhi; G4cin>>Uzhi;
>> 278 // Check for possible quark exchange -----------------------------------
>> 279
>> 280 if(G4UniformRand() < MagQuarkExchange*
>> 281 std::exp(-SlopeQuarkExchange*ProjectileRapidity)) //TargetRapidity))) 1.45
>> 282 {
>> 283 // std::exp(-SlopeQuarkExchange*(ProjectileRapidity - 1.36))) //TargetRapidity))) 1.45
>> 284 //G4cout<<"Q exchange"<<G4endl;
>> 285 //Uzhi_QE=1;
>> 286 G4int NewProjCode(0), NewTargCode(0);
>> 287
>> 288 G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
>> 289
>> 290 // Projectile unpacking --------------------------
>> 291 if(absProjectilePDGcode < 1000 )
>> 292 { // projectile is meson -----------------
>> 293 UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2);
>> 294 } else
>> 295 { // projectile is baryon ----------------
>> 296 UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3);
>> 297 } // End of the hadron's unpacking ----------
>> 298
>> 299 // Target unpacking ------------------------------
>> 300 G4int TargQ1(0), TargQ2(0), TargQ3(0);
>> 301 UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3);
>> 302
>> 303 //G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl;
>> 304 //G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
>> 305 // Sampling of exchanged quarks -------------------
>> 306 G4int ProjExchangeQ(0);
>> 307 G4int TargExchangeQ(0);
>> 308
>> 309 if(absProjectilePDGcode < 1000 )
>> 310 { // projectile is meson -----------------
>> 311
>> 312 if(ProjQ1 > 0 ) // ProjQ1 is quark
>> 313 {
>> 314 G4int Navailable=0;
>> 315 ProjExchangeQ = ProjQ1;
>> 316 if(ProjExchangeQ != TargQ1) Navailable++;
>> 317 if(ProjExchangeQ != TargQ2) Navailable++;
>> 318 if(ProjExchangeQ != TargQ3) Navailable++;
>> 319
>> 320 G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1;
>> 321 //G4cout<<"Navailable Nsampled "<<Navailable<<" "<<Nsampled<<G4endl;
>> 322 Navailable=0;
>> 323 if(ProjExchangeQ != TargQ1)
>> 324 {
>> 325 Navailable++;
>> 326 if(Navailable == Nsampled)
>> 327 {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;}
459 } 328 }
460 } <<
461 #ifdef debugFTFexcitation <<
462 G4cout << "NewProjCode " << NewProjCode <<
463 #endif <<
464 // Decide (with 50% probability) whether <<
465 // but not in the case of charmed and bo <<
466 // there are no excited charmed and bott <<
467 ProjExcited = false; <<
468 if ( aProjQ1 <= 3 && aProjQ2 <= 3 && <<
469 NewProjCode += 2; // Excited meson (J <<
470 ProjExcited = true; <<
471 } <<
472 329
473 // So far we have used the absolute valu << 330 if(ProjExchangeQ != TargQ2)
474 // now look at their signed values to se << 331 {
475 G4int value = ProjQ1, absValue = aProjQ1 << 332 Navailable++;
476 for ( G4int iQ = 0; iQ < 2; ++iQ ) { // << 333 if(Navailable == Nsampled)
477 if ( iQ == 1 ) { << 334 {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;}
478 value = ProjQ2; absValue = aProjQ2; << 335 }
479 } << 336
480 if ( absValue == 2 || absValue == 4 << 337 if(ProjExchangeQ != TargQ3)
481 Qquarks += 2*value/absValue; // u, << 338 {
482 } else { << 339 Navailable++;
483 Qquarks -= value/absValue; // d, << 340 if(Navailable == Nsampled)
>> 341 {TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;}
484 } 342 }
485 } << 343 } else // ProjQ2 is quark
486 // If Qquarks is positive, the sign of N << 344 {
487 // If Qquarks is negative, then the sign << 345 G4int Navailable=0;
488 // If Qquarks is zero, then we need to d << 346 ProjExchangeQ = ProjQ2;
489 // 1. If aProjQ1 and aProjQ2 are the sam << 347 if(ProjExchangeQ != TargQ1) Navailable++;
490 // (because the antiparticle is the s << 348 if(ProjExchangeQ != TargQ2) Navailable++;
491 // 2. If aProjQ1 and aProjQ2 are not the << 349 if(ProjExchangeQ != TargQ3) Navailable++;
492 // we have only two possibilities: << 350
493 // a. aProjQ1 and aProjQ2 are two dif << 351 G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1;
494 // (s,d) or (b,d), or (b,s). In th << 352 //G4cout<<"Navailable Nsampled "<<Navailable<<" "<<Nsampled<<G4endl;
495 // is fine (because the heaviest o << 353 Navailable=0;
496 // to be anti-quark belonging to t << 354 if(ProjExchangeQ != TargQ1)
497 // implies a meson with positive P << 355 {
498 // b. aProjQ1 and aProjQ2 are two dif << 356 Navailable++;
499 // The heaviest of the two (c) has << 357 if(Navailable == Nsampled)
500 // in the projectile, therefore th << 358 {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;}
501 // reverse: 421 -> -421 anti_D0 (c <<
502 if ( Qquarks < 0 || ( Qquarks == 0 && <<
503 NewProjCode *= -1; <<
504 } <<
505 #ifdef debugFTFexcitation <<
506 G4cout << "NewProjCode +2 or 0 " << NewP <<
507 G4cout<<"+++++++++++++++++++++++++++++++ <<
508 G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquark <<
509 G4cout<<"+++++++++++++++++++++++++++++++ <<
510 #endif <<
511 <<
512 // Proj <<
513 TestParticle = G4ParticleTable::GetParti <<
514 if ( ! TestParticle ) continue; <<
515 common.MminProjectile = common.BrW.GetMi <<
516 if ( common.SqrtS - common.M0target < co <<
517 MtestPr = common.BrW.SampleMass( TestPar <<
518 <<
519 #ifdef debugFTFexcitation <<
520 G4cout << "TestParticle Name " << NewPro <<
521 << G4endl <<
522 << "MtestPart MtestPart0 "<<Mtest <<
523 << "M0projectile projectile PDGMa <<
524 << projectile->GetDefinition()->G <<
525 #endif <<
526 <<
527 // Targ <<
528 NewTargCode = NewNucleonId( TargQ1, Targ <<
529 #ifdef debugFTFexcitation <<
530 G4cout << "New TrQ " << TargQ1 << " " << <<
531 << "NewTargCode " << NewTargCode <<
532 #endif <<
533 <<
534 // If the target has not heavy (charm or <<
535 // see whether a Delta isobar can be cre <<
536 if ( TargQ1 <= 3 && TargQ2 <= 3 && T <<
537 if ( TargQ1 != TargQ2 && TargQ1 != T <<
538 if ( G4UniformRand() < 0.5 ) { <<
539 NewTargCode += 2; <<
540 } else if ( G4UniformRand() < 0.75 ) <<
541 NewTargCode = 3122; // Lambda <<
542 } <<
543 } else if ( TargQ1 == TargQ2 && Targ <<
544 NewTargCode += 2; ProjExcited = true <<
545 } else if ( target->GetDefinition()->G <<
546 if ( G4UniformRand() > DeltaProbAtQu <<
547 NewTargCode += 2; ProjExcited = tr <<
548 } <<
549 } else if ( ! ProjExcited && <<
550 G4UniformRand() < DeltaPro <<
551 common.SqrtS > common.M0pr <<
552 G4ParticleTable::GetP <<
553 NewTargCode += 2; // Create Delta i <<
554 } 359 }
555 } <<
556 360
557 // Protection against: << 361 if(ProjExchangeQ != TargQ2)
558 // - Lambda* (i.e. excited Lambda state) << 362 {
559 // - Sigma* (i.e. excited Sigma hyperon << 363 Navailable++;
560 // - Xi* (i.e. excited Xi hyperon st << 364 if(Navailable == Nsampled)
561 if ( NewTargCode == 3124 || // Lambda << 365 {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;}
562 NewTargCode == 3224 || // Sigma*+ NOT <<
563 NewTargCode == 3214 || // Sigma*0 NOT <<
564 NewTargCode == 3114 || // Sigma*- NOT <<
565 NewTargCode == 3324 || // Xi*0 NOT <<
566 NewTargCode == 3314 ) { // Xi*- NOT <<
567 //G4cout << "G4DiffractiveExcitation::Excite <<
568 // << NewTargCode << G4endl; <<
569 NewTargCode -= 2; // Corresponding ground-s <<
570 } <<
571 <<
572 // Special treatment for charmed and bot <<
573 // so we need to transform them by hand <<
574 #ifdef debug_heavyHadrons <<
575 G4int initialNewTargCode = NewTargCode; <<
576 #endif <<
577 if ( NewTargCode == 4322 ) NewTargC <<
578 else if ( NewTargCode == 4312 ) NewTargC <<
579 else if ( NewTargCode == 5312 ) NewTargC <<
580 else if ( NewTargCode == 5322 ) NewTargC <<
581 #ifdef debug_heavyHadrons <<
582 if ( NewTargCode != initialNewTargCode ) <<
583 G4cout << "G4DiffractiveExcitation::ExcitePa <<
584 << "\t target heavy baryon with pdgCo <<
585 << " into pdgCode=" << NewTargCode << <<
586 } <<
587 #endif <<
588 <<
589 TestParticle = G4ParticleTable::GetParti <<
590 if ( ! TestParticle ) continue; <<
591 #ifdef debugFTFexcitation <<
592 G4cout << "New targ " << NewTargCode << <<
593 #endif <<
594 common.MminTarget = common.BrW.GetMinimu <<
595 if ( common.SqrtS - MtestPr < common.Mmi <<
596 MtestTr = common.BrW.SampleMass( TestPar <<
597 <<
598 if ( common.SqrtS > MtestPr + MtestTr ) <<
599 <<
600 } // End of while loop <<
601 if ( attempts >= maxNumberOfAttempts ) ret <<
602 <<
603 if ( MtestPr >= common.Pprojectile.mag() <<
604 common.M0projectile = MtestPr; <<
605 } <<
606 #ifdef debugFTFexcitation <<
607 G4cout << "M0projectile After check " << c <<
608 #endif <<
609 common.M0projectile2 = common.M0projectile <<
610 common.ProjectileDiffStateMinMass = com <<
611 common.ProjectileNonDiffStateMinMass = com <<
612 if ( MtestTr >= common.Ptarget.mag() || <<
613 common.M0target = MtestTr; <<
614 } <<
615 common.M0target2 = common.M0target * commo <<
616 #ifdef debugFTFexcitation <<
617 G4cout << "New targ M0 M0^2 " << common.M0 <<
618 #endif <<
619 common.TargetDiffStateMinMass = common. <<
620 common.TargetNonDiffStateMinMass = common. <<
621 <<
622 } else { // of the if ( common.absProjectil <<
623 <<
624 // If it is a projectile anti-baryon, no q <<
625 // therefore returns immediately 1 (which <<
626 // needs to be continued"). <<
627 if ( common.ProjectilePDGcode < 0 ) return <<
628 <<
629 // Choose randomly, with equal probability <<
630 // projectile or target hadron for selecti <<
631 G4bool isProjectileExchangedQ = false; <<
632 G4int firstQ = TargQ1, secondQ = <<
633 G4int otherFirstQ = ProjQ1, otherSecondQ = <<
634 if ( G4UniformRand() < 0.5 ) { <<
635 isProjectileExchangedQ = true; <<
636 firstQ = ProjQ1; secondQ = Pro <<
637 otherFirstQ = TargQ1; otherSecondQ = Tar <<
638 } <<
639 // Choose randomly, with equal probability <<
640 // selected (projectile or target) hadron <<
641 G4int exchangedQ = 0; <<
642 G4double Ksi = G4UniformRand(); <<
643 if ( Ksi < 0.333333 ) { <<
644 exchangedQ = firstQ; <<
645 } else if ( 0.333333 <= Ksi && Ksi < 0.6 <<
646 exchangedQ = secondQ; <<
647 } else { <<
648 exchangedQ = thirdQ; <<
649 } <<
650 #ifdef debugFTFexcitation <<
651 G4cout << "Exchange Qs isProjectile Q " << <<
652 #endif <<
653 // The exchanged quarks (one of the projec <<
654 // are always accepted if they have differ <<
655 // same flavour) they are accepted only wi <<
656 G4double probSame = theParameters->GetProb <<
657 const G4int MaxCount = 100; <<
658 G4int count = 0, otherExchangedQ = 0; <<
659 do { <<
660 if ( exchangedQ != otherFirstQ || G4Un <<
661 otherExchangedQ = otherFirstQ; otherFi <<
662 } else { <<
663 if ( exchangedQ != otherSecondQ || G <<
664 otherExchangedQ = otherSecondQ; othe <<
665 } else { <<
666 if ( exchangedQ != otherThirdQ || <<
667 otherExchangedQ = otherThirdQ; oth <<
668 } <<
669 } 366 }
670 } <<
671 } while ( otherExchangedQ == 0 && ++coun <<
672 if ( count >= MaxCount ) return returnCode <<
673 // Swap (between projectile and target had <<
674 // as "exchanged" quarks. <<
675 if ( Ksi < 0.333333 ) { <<
676 firstQ = exchangedQ; <<
677 } else if ( 0.333333 <= Ksi && Ksi < 0.6 <<
678 secondQ = exchangedQ; <<
679 } else { <<
680 thirdQ = exchangedQ; <<
681 } <<
682 if ( isProjectileExchangedQ ) { <<
683 ProjQ1 = firstQ; ProjQ2 = secondQ; <<
684 TargQ1 = otherFirstQ; TargQ2 = otherSeco <<
685 } else { <<
686 TargQ1 = firstQ; TargQ2 = secondQ; <<
687 ProjQ1 = otherFirstQ; ProjQ2 = otherSeco <<
688 } <<
689 #ifdef debugFTFexcitation <<
690 G4cout << "Exchange Qs Pr Tr " << ( isPro <<
691 << " " << ( isProjectileExchangedQ <<
692 #endif <<
693 <<
694 NewProjCode = NewNucleonId( ProjQ1, ProjQ2 <<
695 NewTargCode = NewNucleonId( TargQ1, TargQ2 <<
696 // Decide whether the new projectile hadro <<
697 // then decide whether the new target hadr <<
698 // Notice that a Delta particle has the la <<
699 // whereas a nucleon has "2" (because its <<
700 // If the new projectile hadron or the new <<
701 // constituent quark, then skip this part <<
702 // excited charm and bottom hadrons). <<
703 for ( G4int iHadron = 0; iHadron < 2; iHad <<
704 // First projectile hadron, then target <<
705 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, <<
706 G4double massConstraint = common.M0targe <<
707 G4bool isHadronADelta = ( projectile->Ge <<
708 if ( iHadron == 1 ) { // Target hadron <<
709 codeQ1 = TargQ1, codeQ2 = TargQ2, code <<
710 massConstraint = common.M0projectile; <<
711 isHadronADelta = ( target->GetDefiniti <<
712 } <<
713 if ( codeQ1 > 3 || codeQ2 > 3 || cod <<
714 if ( codeQ1 == codeQ2 && codeQ1 == cod <<
715 newHadCode += 2; // Delta++ (uuu) or <<
716 } else if ( isHadronADelta ) { // Hadro <<
717 if ( G4UniformRand() > DeltaProbAtQuar <<
718 newHadCode += 2; // Delta+ (uud) or <<
719 } else { <<
720 newHadCode += 0; // No delta (so th <<
721 } <<
722 } else { // Hadron (projectile or targe <<
723 if ( G4UniformRand() < DeltaProbAtQuar <<
724 common.SqrtS > G4ParticleTable::G <<
725 + massConstraint ) <<
726 newHadCode += 2; // Delta+ (uud) or <<
727 } else { <<
728 newHadCode += 0; // No delta (so th <<
729 } <<
730 } <<
731 if ( iHadron == 0 ) { // Projectile had <<
732 NewProjCode = newHadCode; <<
733 } else { // Target hadron <<
734 NewTargCode = newHadCode; <<
735 } <<
736 } <<
737 #ifdef debugFTFexcitation <<
738 G4cout << "NewProjCode NewTargCode " << Ne <<
739 #endif <<
740 <<
741 // Protection against: <<
742 // - Lambda* (i.e. excited Lambda state) <<
743 // - Sigma* (i.e. excited Sigma hyperon s <<
744 // - Xi* (i.e. excited Xi hyperon stat <<
745 if ( NewProjCode == 3124 || // Lambda* <<
746 NewProjCode == 3224 || // Sigma*+ NOT ex <<
747 NewProjCode == 3214 || // Sigma*0 NOT ex <<
748 NewProjCode == 3114 || // Sigma*- NOT ex <<
749 NewProjCode == 3324 || // Xi*0 NOT ex <<
750 NewProjCode == 3314 ) { // Xi*- NOT ex <<
751 //G4cout << "G4DiffractiveExcitation::Ex <<
752 // << NewProjCode << G4endl; <<
753 NewProjCode -= 2; // Corresponding grou <<
754 } <<
755 if ( NewTargCode == 3124 || // Lambda* <<
756 NewTargCode == 3224 || // Sigma*+ NOT ex <<
757 NewTargCode == 3214 || // Sigma*0 NOT ex <<
758 NewTargCode == 3114 || // Sigma*- NOT ex <<
759 NewTargCode == 3324 || // Xi*0 NOT ex <<
760 NewTargCode == 3314 ) { // Xi*- NOT ex <<
761 //G4cout << "G4DiffractiveExcitation::Ex <<
762 // << NewTargCode << G4endl; <<
763 NewTargCode -= 2; // Corresponding grou <<
764 } <<
765 367
766 // Special treatment for charmed and botto << 368 if(ProjExchangeQ != TargQ3)
767 // so we need to transform them by hand to << 369 {
768 #ifdef debug_heavyHadrons << 370 Navailable++;
769 G4int initialNewProjCode = NewProjCode, in << 371 if(Navailable == Nsampled)
770 #endif << 372 {TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;}
771 if ( NewProjCode == 4322 ) NewProjCod <<
772 else if ( NewProjCode == 4312 ) NewProjCod <<
773 else if ( NewProjCode == 5312 ) NewProjCod <<
774 else if ( NewProjCode == 5322 ) NewProjCod <<
775 if ( NewTargCode == 4322 ) NewTargCod <<
776 else if ( NewTargCode == 4312 ) NewTargCod <<
777 else if ( NewTargCode == 5312 ) NewTargCod <<
778 else if ( NewTargCode == 5322 ) NewTargCod <<
779 #ifdef debug_heavyHadrons <<
780 if ( NewProjCode != initialNewProjCode || <<
781 G4cout << "G4DiffractiveExcitation::Exci <<
782 << "\t heavy baryon into an exist <<
783 if ( NewProjCode != initialNewProjCode ) <<
784 G4cout << "\t \t projectile baryon with pdgC <<
785 << " into pdgCode=" << NewProjCode << <<
786 } <<
787 if ( NewTargCode != initialNewTargCode ) <<
788 G4cout << "\t \t target baryon with pd <<
789 << " into pdgCode=" << NewTargCode << <<
790 } <<
791 } <<
792 #endif <<
793 <<
794 // Sampling of the masses of the projectil <<
795 // Because of energy conservation, the ord <<
796 // randomly, half of the time we sample fi <<
797 // then the projectile nucleon mass, and t <<
798 // sample first the projectile nucleon mas <<
799 G4VSplitableHadron* firstHadron = target; <<
800 G4VSplitableHadron* secondHadron = project <<
801 G4int firstHadronCode = NewTargCode, secon <<
802 G4double massConstraint = common.M0project <<
803 G4bool isFirstTarget = true; <<
804 if ( G4UniformRand() < 0.5 ) { <<
805 // Sample first the projectile nucleon m <<
806 firstHadron = projectile; secondHa <<
807 firstHadronCode = NewProjCode; secondHa <<
808 massConstraint = common.M0target; <<
809 isFirstTarget = false; <<
810 } <<
811 G4double Mtest1st = 0.0, Mtest2nd = 0.0, M <<
812 for ( int iSamplingCase = 0; iSamplingCase <<
813 G4VSplitableHadron* aHadron = firstHadro <<
814 G4int aHadronCode = firstHadronCode; <<
815 if ( iSamplingCase == 1 ) { // Second n <<
816 aHadron = secondHadron; aHadronCode = <<
817 } <<
818 G4double MtestHadron = 0.0, MminHadron = <<
819 if ( aHadron->GetStatus() == 1 || aHad <<
820 TestParticle = G4ParticleTable::GetPar <<
821 if ( ! TestParticle ) return returnCod <<
822 MminHadron = common.BrW.GetMinimumMass <<
823 if ( common.SqrtS - massConstraint < M <<
824 if ( TestParticle->GetPDGWidth() == 0. <<
825 MtestHadron = common.BrW.SampleMass( <<
826 } else { <<
827 const G4int maxNumberOfAttempts = 50 <<
828 G4int attempts = 0; <<
829 while ( attempts < maxNumberOfAttemp <<
830 attempts++; <<
831 MtestHadron = common.BrW.SampleMas <<
832 <<
833 if ( common.SqrtS < MtestHadron + <<
834 continue; // Kinematically unac <<
835 } else { <<
836 break; // Kinematically acce <<
837 } <<
838 } <<
839 if ( attempts >= maxNumberOfAttempts <<
840 } 373 }
841 } << 374 } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark
842 if ( iSamplingCase == 0 ) { <<
843 Mtest1st = MtestHadron; Mmin1st = Mmi <<
844 } else { <<
845 Mtest2nd = MtestHadron; Mmin2nd = Mmi <<
846 } <<
847 } // End for loop on the two sampling cas <<
848 if ( isFirstTarget ) { <<
849 MtestTr = Mtest1st; MtestPr = Mtest2n <<
850 common.MminTarget = Mmin1st; common.Mmi <<
851 } else { <<
852 MtestTr = Mtest2nd; MtestPr = Mtest1s <<
853 common.MminTarget = Mmin2nd; common.Mmi <<
854 } <<
855 <<
856 if ( MtestPr != 0.0 ) { <<
857 common.M0projectile = MtestPr; <<
858 common.M0projectile2 = common.M0projecti <<
859 common.ProjectileDiffStateMinMass = c <<
860 common.ProjectileNonDiffStateMinMass = c <<
861 } <<
862 if ( MtestTr != 0.0 ) { <<
863 common.M0target = MtestTr; <<
864 common.M0target2 = common.M0target * com <<
865 common.TargetDiffStateMinMass = commo <<
866 common.TargetNonDiffStateMinMass = commo <<
867 } <<
868 375
869 } // End of if ( common.absProjectilePDGcod << 376 //G4cout<<"Exch Pr Tr "<<ProjExchangeQ<<" "<<TargExchangeQ<<G4endl;
>> 377 //G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl;
>> 378 //G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
>> 379
>> 380 G4int aProjQ1=std::abs(ProjQ1);
>> 381 G4int aProjQ2=std::abs(ProjQ2);
>> 382 if(aProjQ1 == aProjQ2) {NewProjCode = 111;} // Pi0-meson
>> 383 else // |ProjQ1| # |ProjQ2|
>> 384 {
>> 385 if(aProjQ1 > aProjQ2) {NewProjCode = aProjQ1*100+aProjQ2*10+1;}
>> 386 else {NewProjCode = aProjQ2*100+aProjQ1*10+1;}
>> 387 // NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode);
>> 388 }
>> 389
>> 390 G4bool ProjExcited=false;
>> 391 //G4cout<<"NewProjCode "<<NewProjCode<<G4endl;
>> 392 if(G4UniformRand() < 0.5)
>> 393 {
>> 394 NewProjCode +=2; // Excited Pi0-meson
>> 395 ProjExcited=true;
>> 396 }
>> 397 if(aProjQ1 != aProjQ2) NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode); // Uzhi
>> 398 //G4cout<<"NewProjCode +2 or 0 "<<NewProjCode<<G4endl;
870 399
871 // If we assume that final state hadrons aft << 400 G4ParticleDefinition* TestParticle=0;
872 // in the ground states, we have to put << 401 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode);
873 if ( common.SqrtS < common.M0projectile + co << 402 //G4cout<<"TestParticle ? "<<TestParticle<<G4endl;
874 << 403
875 common.PZcms2 = ( sqr( common.S ) + sqr( com << 404 if(TestParticle)
876 - 2.0 * ( common.S * ( com << 405 {
877 + common.M0proje << 406 G4double MtestPart= // 31.05.2012
878 #ifdef debugFTFexcitation << 407 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
879 G4cout << "At the end// NewProjCode " << New << 408 /*
880 << "At the end// NewTargCode " << New << 409 G4cout<<"TestParticle Name "<<NewProjCode<<" "<<TestParticle->GetParticleName()<<G4endl;
881 << "M0pr M0tr SqS " << common.M0pro << 410 G4cout<<"MtestPart M0projectile projectile->GetDefinition()->GetPDGMass() "<<MtestPart<<" "<<M0projectile<<" "<<projectile->GetDefinition()->GetPDGMass()<<G4endl;
882 << common.SqrtS << G4endl << 411 G4bool Test =M0projectile <= projectile->GetDefinition()->GetPDGMass();
883 << "M0pr2 M0tr2 SqS " << common.M0pro << 412 G4cout<<"M0projectile <= projectile->GetDefinition()->GetPDGMass() "<<Test<<G4endl;
884 << common.SqrtS << G4endl << 413 */
885 << "PZcms2 after the change " << comm << 414
886 #endif << 415 if(MtestPart > M0projectile)
887 if ( common.PZcms2 < 0.0 ) return returnCode << 416 {M0projectile = MtestPart;}
888 << 417 else
889 projectile->SetDefinition( G4ParticleTable:: << 418 {
890 target->SetDefinition( G4ParticleTable::GetP << 419 if(std::abs(M0projectile - projectile->GetDefinition()->GetPDGMass()) < 140.*MeV)
891 common.PZcms = std::sqrt( common.PZcms2 ); << 420 {M0projectile = MtestPart;}
892 common.Pprojectile.setPz( common.PZcms ); << 421 }
893 common.Pprojectile.setE( std::sqrt( common.M << 422 //G4cout<<"M0projectile After check "<<M0projectile<<G4endl;
894 common.Ptarget.setPz( -common.PZcms ); << 423 M0projectile2 = M0projectile * M0projectile;
895 common.Ptarget.setE( std::sqrt( common.M0tar << 424
896 #ifdef debugFTFexcitation << 425 ProjectileDiffStateMinMass =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
897 G4cout << "Proj Targ and Proj+Targ in CMS" < << 426 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
898 << common.Ptarget << G4endl << common << 427 } else
899 #endif << 428 {return false;}
900 << 429
901 if ( projectile->GetStatus() != 0 ) projecti << 430 //G4cout<<"New TrQ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
902 if ( target->GetStatus() != 0 ) target-> << 431 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
903 << 432 //G4cout<<"NewTargCode "<<NewTargCode<<G4endl;
904 // Check for possible excitation of the part << 433
905 if ( common.SqrtS < common.M0projectile + co << 434 // if( (TargQ1 != TargQ2) && (TargQ1 != TargQ3) && (TargQ2 != TargQ3) // Lambda or Sigma0
906 common.SqrtS < common.ProjectileDiffSta << 435 // {if(G4UniformRand() < 0.5) NewTargCode=
907 common.ProbOfDiffraction == 0.0 ) commo << 436
908 << 437
909 if ( G4UniformRand() > common.ProbExc ) { / << 438 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
910 #ifdef debugFTFexcitation << 439 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2; //Create Delta isobar
911 G4cout << "Make elastic scattering of new << 440 ProjExcited=true;}
912 #endif << 441 else if( target->GetDefinition()->GetPDGiIsospin() == 3 ) //Delta was the target
913 common.Pprojectile.transform( common.toLab << 442 { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2; //Save Delta isobar
914 common.Ptarget.transform( common.toLab ); << 443 ProjExcited=true;}
915 projectile->Set4Momentum( common.Pprojecti << 444 else {} // De-excite initial Delta isobar
916 target->Set4Momentum( common.Ptarget ); << 445 }
917 G4bool Result = theElastic->ElasticScatter << 446
918 #ifdef debugFTFexcitation << 447 // else if((!CreateDelta) &&
919 G4cout << "Result of el. scatt " << Result << 448 else if((!ProjExcited) &&
920 << G4endl << projectile->Get4Moment << 449 (G4UniformRand() < DeltaProbAtQuarkExchange) && //Nucleon was the target
921 << projectile->Get4Momentum() + tar << 450 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar
922 << (projectile->Get4Momentum() + ta << 451 // else if( CreateDelta) {NewTargCode +=2;}
923 #endif << 452 else {} //Save initial nucleon
924 if ( Result ) returnCode = 0; // successf << 453
925 return returnCode; << 454 // target->SetDefinition( // Fix 15.12.09
926 } << 455 // G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09
>> 456
>> 457 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
>> 458 //G4cout<<"New targ "<<NewTargCode<<" "<<TestParticle->GetParticleName()<<G4endl;
>> 459 if(TestParticle)
>> 460 {
>> 461 G4double MtestPart= // 31.05.2012
>> 462 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
>> 463
>> 464 if(MtestPart > M0target)
>> 465 {M0target=MtestPart;}
>> 466 else
>> 467 {
>> 468 if(std::abs(M0target - target->GetDefinition()->GetPDGMass()) < 140.*MeV)
>> 469 {M0target=MtestPart;}
>> 470 }
>> 471
>> 472 TargetDiffStateMinMass =M0target+220.*MeV; //220 MeV=m_pi+80 MeV;
>> 473 TargetNonDiffStateMinMass=M0target+220.*MeV; //220 MeV=m_pi+80 MeV;
>> 474 } else
>> 475 {return false;}
>> 476 } else
>> 477 { // projectile is baryon ----------------
>> 478 //=========================================================================
>> 479 G4double Same=theParameters->GetProbOfSameQuarkExchange(); //0.3; //0.5; 0.
>> 480 G4bool ProjDeltaHasCreated(false);
>> 481 G4bool TargDeltaHasCreated(false);
>> 482
>> 483 G4double Ksi=G4UniformRand();
>> 484 if(G4UniformRand() < 0.5) // Sampling exchange quark from proj. or targ.
>> 485 { // Sampling exchanged quark from the projectile ---
>> 486 if( Ksi < 0.333333 )
>> 487 {ProjExchangeQ = ProjQ1;}
>> 488 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
>> 489 {ProjExchangeQ = ProjQ2;}
>> 490 else
>> 491 {ProjExchangeQ = ProjQ3;}
>> 492
>> 493 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
>> 494 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same))
>> 495 {
>> 496 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
>> 497 } else
>> 498 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same))
>> 499 {
>> 500 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
>> 501 } else
>> 502 {
>> 503 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
>> 504 }
927 505
928 #ifdef debugFTFexcitation << 506 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
929 G4cout << "Make excitation of new hadrons" < << 507 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
930 #endif << 508 if( Ksi < 0.333333 )
931 << 509 {ProjQ1=ProjExchangeQ;}
932 // Redefinition of ProbOfDiffraction because << 510 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
933 common.ProbOfDiffraction = common.ProbProjec << 511 {ProjQ2=ProjExchangeQ;}
934 if ( common.ProbOfDiffraction != 0.0 ) { << 512 else
935 common.ProbProjectileDiffraction /= common << 513 {ProjQ3=ProjExchangeQ;}
936 common.ProbTargetDiffraction /= common << 514
937 } << 515 } else
>> 516 { // Sampling exchanged quark from the target -------
>> 517 if( Ksi < 0.333333 )
>> 518 {TargExchangeQ = TargQ1;}
>> 519 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
>> 520 {TargExchangeQ = TargQ2;}
>> 521 else
>> 522 {TargExchangeQ = TargQ3;}
>> 523
>> 524 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same))
>> 525 {
>> 526 ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
>> 527 } else
>> 528 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same))
>> 529 {
>> 530 ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
>> 531 } else
>> 532 {
>> 533 ProjExchangeQ = ProjQ3; ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
>> 534 }
938 535
939 return returnCode = 1; // successfully comp << 536 if( Ksi < 0.333333 )
940 } << 537 {TargQ1=TargExchangeQ;}
>> 538 else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
>> 539 {TargQ2=TargExchangeQ;}
>> 540 else
>> 541 {TargQ3=TargExchangeQ;}
>> 542
>> 543 } // End of sampling baryon
>> 544
>> 545 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
>> 546
>> 547 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
>> 548 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
>> 549 { if(G4UniformRand() > DeltaProbAtQuarkExchange)
>> 550 {NewProjCode +=2; ProjDeltaHasCreated=true;}
>> 551 else {NewProjCode +=0; ProjDeltaHasCreated=false;}
>> 552 }
>> 553 else // Projectile was Nucleon
>> 554 {
>> 555 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target))
>> 556 {NewProjCode +=2; ProjDeltaHasCreated=true;}
>> 557 else {NewProjCode +=0; ProjDeltaHasCreated=false;}
>> 558 }
>> 559
>> 560 //G4ParticleDefinition* NewTestParticle=
>> 561 // G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode);
>> 562 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl;
>> 563 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewProjCode=TestParticleID;}
>> 564
>> 565 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
941 566
942 //-------------------------------------------- << 567 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
943 568
944 G4bool G4DiffractiveExcitation:: << 569 //G4cout<<"TargQ1, TargQ2, TargQ3 "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<" "<<NewTargCode<<G4endl;
945 ExciteParticipants_doDiffraction( G4VSplitable <<
946 G4FTFParamet <<
947 G4Diffractiv <<
948 // Second of the three utility methods used <<
949 // it does the sampling for the diffraction <<
950 <<
951 G4bool isProjectileDiffraction = false; <<
952 if ( G4UniformRand() < common.ProbProjectile <<
953 isProjectileDiffraction = true; <<
954 #ifdef debugFTFexcitation <<
955 G4cout << "projectile diffraction" << G4en <<
956 #endif <<
957 common.ProjMassT2 = common.ProjectileDiffS <<
958 common.ProjMassT = common.ProjectileDiffS <<
959 common.TargMassT2 = common.M0target2; <<
960 common.TargMassT = common.M0target; <<
961 } else { <<
962 #ifdef debugFTFexcitation <<
963 G4cout << "Target diffraction" << G4endl; <<
964 #endif <<
965 common.ProjMassT2 = common.M0projectile2; <<
966 common.ProjMassT = common.M0projectile; <<
967 common.TargMassT2 = common.TargetDiffState <<
968 common.TargMassT = common.TargetDiffState <<
969 } <<
970 570
971 // Check whether the interaction is possible << 571 //TestParticleID=NewTargCode;
972 if ( common.SqrtS < common.ProjMassT + commo << 572 //TestParticleMass=DBL_MAX;
973 <<
974 common.PZcms2 = ( sqr( common.S ) + sqr( com <<
975 - 2.0 * ( common.S * ( com <<
976 + common.ProjMassT2 <<
977 if ( common.PZcms2 < 0.0 ) return false; <<
978 common.maxPtSquare = common.PZcms2; <<
979 <<
980 G4double DiffrAveragePt2 = theParameters->Ge <<
981 G4bool loopCondition = true; <<
982 G4int whilecount = 0; <<
983 do { // Generate pt and mass of projectile <<
984 <<
985 whilecount++; <<
986 if ( whilecount > 1000 ) { <<
987 common.Qmomentum = G4LorentzVector( 0.0, <<
988 return false; // Ignore this interacti <<
989 }; <<
990 <<
991 common.Qmomentum = G4LorentzVector( Gaussi <<
992 common.Pt2 = G4ThreeVector( common.Qmoment <<
993 if ( isProjectileDiffraction ) { // proje <<
994 common.ProjMassT2 = common.ProjectileDif <<
995 common.TargMassT2 = common.M0target2 + c <<
996 } else { // targe <<
997 common.ProjMassT2 = common.M0projectile2 <<
998 common.TargMassT2 = common.TargetDiffSta <<
999 } <<
1000 common.ProjMassT = std::sqrt( common.Proj <<
1001 common.TargMassT = std::sqrt( common.Targ <<
1002 if ( common.SqrtS < common.ProjMassT + co <<
1003 <<
1004 common.PZcms2 = ( sqr( common.S ) + sqr( <<
1005 - 2.0 * ( common.S * ( <<
1006 + common.Proj <<
1007 if ( common.PZcms2 < 0.0 ) continue; <<
1008 <<
1009 common.PZcms = std::sqrt( common.PZcms2 ) <<
1010 if ( isProjectileDiffraction ) { // proj <<
1011 common.PMinusMin = std::sqrt( common.Pr <<
1012 common.PMinusMax = common.SqrtS - commo <<
1013 common.PMinusNew = ChooseP( common.PMin <<
1014 common.TMinusNew = common.SqrtS - commo <<
1015 common.Qminus = common.Ptarget.minus() <<
1016 common.TPlusNew = common.TargMassT2 / c <<
1017 common.Qplus = common.Ptarget.plus() - <<
1018 common.Qmomentum.setPz( (common.Qplus - <<
1019 common.Qmomentum.setE( (common.Qplus + <<
1020 loopCondition = ( ( common.Pprojectile <<
1021 common.ProjectileDiff <<
1022 } else { // targ <<
1023 common.TPlusMin = std::sqrt( common.Tar <<
1024 common.TPlusMax = common.SqrtS - common <<
1025 common.TPlusNew = ChooseP( common.TPlus <<
1026 common.PPlusNew = common.SqrtS - common <<
1027 common.Qplus = common.PPlusNew - common <<
1028 common.PMinusNew = common.ProjMassT2 / <<
1029 common.Qminus = common.PMinusNew - comm <<
1030 common.Qmomentum.setPz( (common.Qplus - <<
1031 common.Qmomentum.setE( (common.Qplus + <<
1032 loopCondition = ( ( common.Ptarget - c <<
1033 common.TargetDiffSta <<
1034 } <<
1035 573
1036 } while ( loopCondition ); /* Loop checkin << 574 //TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
1037 // Repeat the sampling because ther << 575 //if(TestParticle) TestParticleMass=TestParticle->GetPDGMass();
>> 576
>> 577 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;}
>> 578 else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta
>> 579 { if(G4UniformRand() > DeltaProbAtQuarkExchange)
>> 580 {NewTargCode +=2; TargDeltaHasCreated=true;}
>> 581 else {NewTargCode +=0; TargDeltaHasCreated=false;}
>> 582 }
>> 583 else // Target was Nucleon
>> 584 {
>> 585 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass))
>> 586 {NewTargCode +=2; TargDeltaHasCreated=true;}
>> 587 else {NewTargCode +=0; TargDeltaHasCreated=false;}
>> 588 }
>> 589
>> 590 //NewTestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
>> 591 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl;
>> 592 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewTargCode=TestParticleID;}
>> 593
>> 594 //G4cout<<"NewProjCode NewTargCode "<<NewProjCode<<" "<<NewTargCode<<G4endl;
>> 595 //G4int Uzhi; G4cin>>Uzhi;
>> 596
>> 597 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
>> 598 { // Nothing was changed! It is not right!?
>> 599 }
>> 600 // Forming baryons --------------------------------------------------
>> 601 if(ProjDeltaHasCreated) {ProbProjectileDiffraction=1.; ProbTargetDiffraction=0.;}
>> 602 if(TargDeltaHasCreated) {ProbProjectileDiffraction=0.; ProbTargetDiffraction=1.;}
>> 603 if(ProjDeltaHasCreated)
>> 604 {
>> 605 G4double MtestPart= // 31.05.2012
>> 606 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
>> 607
>> 608 if(MtestPart >= M0projectile) // 31.05.2012
>> 609 { // 31.05.2012
>> 610 M0projectile = MtestPart; // 31.05.2012
>> 611 M0projectile2 = M0projectile * M0projectile; // 31.05.2012
>> 612 } // 31.05.2012
>> 613
>> 614 ProjectileDiffStateMinMass =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
>> 615 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
>> 616 }
>> 617
>> 618 // if(M0target <
>> 619 // (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
>> 620 if(TargDeltaHasCreated)
>> 621 {
>> 622 G4double MtestPart= // 31.05.2012
>> 623 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
>> 624
>> 625 if(MtestPart >=M0target) // 31.05.2012
>> 626 { // 31.05.2012
>> 627 M0target=MtestPart; // 31.05.2012
>> 628 M0target2 = M0target * M0target; // 31.05.2012
>> 629 } // 31.05.2012
>> 630
>> 631 TargetDiffStateMinMass =M0target+210.*MeV; //210 MeV=m_pi+70 MeV;
>> 632 TargetNonDiffStateMinMass=M0target+210.*MeV; //210 MeV=m_pi+70 MeV;
>> 633 }
>> 634 } // End of if projectile is baryon ---------------------------
>> 635
>> 636 //G4cout<<"At end// NewProjCode "<<NewProjCode<<G4endl;
>> 637 //G4cout<<"At end// NewTargCode "<<NewTargCode<<G4endl;
>> 638
>> 639 // If we assume that final state hadrons after the charge exchange will be
>> 640 // in the ground states, we have to put ----------------------------------
>> 641 //G4cout<<"M0pr M0tr SqS "<<M0projectile<<" "<<M0target<<" "<<SqrtS<<G4endl;
>> 642
>> 643 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
>> 644 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
>> 645 /4./S;
>> 646 //G4cout<<"PZcms2 1 "<<PZcms2<<G4endl<<G4endl;
>> 647 if(PZcms2 < 0) {return false;} // It can be if energy is not sufficient for Delta
>> 648 //----------------------------------------------------------
>> 649 projectile->SetDefinition(
>> 650 G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode));
>> 651
>> 652 target->SetDefinition(
>> 653 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));
>> 654 //----------------------------------------------------------
>> 655 PZcms = std::sqrt(PZcms2);
>> 656
>> 657 Pprojectile.setPz( PZcms);
>> 658 Pprojectile.setE(std::sqrt(M0projectile2+PZcms2));
>> 659
>> 660 Ptarget.setPz( -PZcms);
>> 661 Ptarget.setE(std::sqrt(M0target2+PZcms2));
>> 662
>> 663 // ----------------------------------------------------------
>> 664
>> 665 if(absProjectilePDGcode < 1000)
>> 666 { // For projectile meson
>> 667 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity);
>> 668 Wexcit=0.;
>> 669 if(G4UniformRand() > Wexcit)
>> 670 { // Make elastic scattering
>> 671 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl;
>> 672 Pprojectile.transform(toLab);
>> 673 Ptarget.transform(toLab);
>> 674
>> 675 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
>> 676 projectile->SetPosition(target->GetPosition());
>> 677
>> 678 projectile->Set4Momentum(Pprojectile);
>> 679 target->Set4Momentum(Ptarget);
>> 680
>> 681 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
>> 682 //G4cout<<"Result of el. scatt "<<Result<<G4endl;
>> 683 return Result;
>> 684 } // end of if(Make elastic scattering for projectile meson?)
>> 685 } else
>> 686 { // For projectile baryon
>> 687 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity);
>> 688 //Wexcit=0.;
>> 689 if(G4UniformRand() > Wexcit)
>> 690 { // Make elastic scattering
>> 691 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl;
>> 692 Pprojectile.transform(toLab);
>> 693 Ptarget.transform(toLab);
>> 694
>> 695 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
>> 696 projectile->SetPosition(target->GetPosition());
>> 697
>> 698 projectile->Set4Momentum(Pprojectile);
>> 699 target->Set4Momentum(Ptarget);
>> 700
>> 701 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
>> 702 return Result;
>> 703 } // end of if(Make elastic scattering for projectile baryon?)
>> 704 }
>> 705 //G4cout<<"Make excitation of new hadrons"<<G4endl;
>> 706 } // End of charge exchange part ------------------------------
>> 707
>> 708 // ------------------------------------------------------------------
>> 709 G4double ProbOfDiffraction=ProbProjectileDiffraction + ProbTargetDiffraction;
>> 710 /*
>> 711 G4cout<<"Excite --------------------"<<G4endl;
>> 712 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<G4endl;
>> 713 G4cout<<"Targ "<<M0target <<" "<<TargetDiffStateMinMass <<" "<<TargetNonDiffStateMinMass<<G4endl;
>> 714 G4cout<<"SqrtS "<<SqrtS<<G4endl;
>> 715
>> 716 G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl;
>> 717 G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl;
>> 718 //G4int Uzhi; G4cin>>Uzhi;
>> 719 */
>> 720 /*
>> 721 if(ProjectileNonDiffStateMinMass + TargetNonDiffStateMinMass > SqrtS) // 24.07.10
>> 722 {
>> 723 if(ProbOfDiffraction!=0.)
>> 724 {
>> 725 ProbProjectileDiffraction/=ProbOfDiffraction;
>> 726 ProbOfDiffraction=1.;
>> 727 } else {return false;}
>> 728 }
>> 729
>> 730 */
>> 731
>> 732 if(ProbOfDiffraction!=0.)
>> 733 {
>> 734 ProbProjectileDiffraction/=ProbOfDiffraction;
>> 735 }
>> 736 else
>> 737 {
>> 738 ProbProjectileDiffraction=0.;
>> 739 }
>> 740
>> 741 //G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl;
>> 742
>> 743 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass *
>> 744 ProjectileDiffStateMinMass;
>> 745 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
>> 746 ProjectileNonDiffStateMinMass;
>> 747
>> 748 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass *
>> 749 TargetDiffStateMinMass;
>> 750 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass *
>> 751 TargetNonDiffStateMinMass;
>> 752
>> 753 G4double Pt2;
>> 754 G4double ProjMassT2, ProjMassT;
>> 755 G4double TargMassT2, TargMassT;
>> 756 G4double PMinusMin, PMinusMax;
>> 757 // G4double PPlusMin , PPlusMax;
>> 758 G4double TPlusMin , TPlusMax;
>> 759 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
>> 760
>> 761 G4LorentzVector Qmomentum;
>> 762 G4double Qminus, Qplus;
>> 763
>> 764 G4int whilecount=0;
>> 765
>> 766 // Choose a process ---------------------------
>> 767 //ProbOfDiffraction=1.; // Uzhi Difr
>> 768 //ProbProjectileDiffraction=1.; // Uzhi
>> 769 if(G4UniformRand() < ProbOfDiffraction)
>> 770 {
>> 771 if(G4UniformRand() < ProbProjectileDiffraction)
>> 772 { //-------- projectile diffraction ---------------
>> 773 //G4cout<<"projectile diffraction"<<G4endl;
>> 774
>> 775 do {
>> 776 /*
>> 777 Uzhi_projectilediffraction=1;
>> 778 Uzhi_targetdiffraction=0;
>> 779 Uzhi_Mx2=1.;
>> 780 */
>> 781 // Generate pt
>> 782 // if (whilecount++ >= 500 && (whilecount%100)==0)
>> 783 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
>> 784 // << ", loop count/ maxPtSquare : "
>> 785 // << whilecount << " / " << maxPtSquare << G4endl;
>> 786
>> 787 // whilecount++;
>> 788 if (whilecount > 1000 )
>> 789 {
>> 790 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
>> 791 target->SetStatus(2); return false; // Ignore this interaction
>> 792 };
>> 793
>> 794 // --------------- Check that the interaction is possible -----------
>> 795 ProjMassT2=ProjectileDiffStateMinMass2;
>> 796 ProjMassT =ProjectileDiffStateMinMass;
>> 797
>> 798 TargMassT2=M0target2;
>> 799 TargMassT =M0target;
>> 800 //G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl;
>> 801 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
>> 802 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
>> 803 /4./S;
>> 804
>> 805 //G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl;
>> 806 if(PZcms2 < 0 )
>> 807 {
>> 808 target->SetStatus(2);
>> 809 return false;
>> 810 }
>> 811 maxPtSquare=PZcms2;
>> 812
>> 813 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
>> 814 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
>> 815
>> 816 ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
>> 817 ProjMassT =std::sqrt(ProjMassT2);
>> 818
>> 819 TargMassT2=M0target2+Pt2;
>> 820 TargMassT =std::sqrt(TargMassT2);
>> 821
>> 822 //G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl;
>> 823
>> 824 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
>> 825 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
>> 826 /4./S;
>> 827 //G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl;
>> 828 if(PZcms2 < 0 ) continue;
>> 829 PZcms =std::sqrt(PZcms2);
>> 830
>> 831 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
>> 832 PMinusMax=SqrtS-TargMassT;
>> 833
>> 834 PMinusNew=ChooseP(PMinusMin, PMinusMax);
>> 835 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
>> 836 //PMinusNew=1./sqr(1./std::sqrt(PMinusMin)-G4UniformRand()*(1./std::sqrt(PMinusMin)-1./std::sqrt(PMinusMax)));
>> 837
>> 838 TMinusNew=SqrtS-PMinusNew;
>> 839 Qminus=Ptarget.minus()-TMinusNew;
>> 840 TPlusNew=TargMassT2/TMinusNew;
>> 841 Qplus=Ptarget.plus()-TPlusNew;
>> 842
>> 843 Qmomentum.setPz( (Qplus-Qminus)/2 );
>> 844 Qmomentum.setE( (Qplus+Qminus)/2 );
>> 845
>> 846 } while ((Pprojectile+Qmomentum).mag2() < ProjectileDiffStateMinMass2); //||
>> 847 //Repeat the sampling because there was not any excitation
>> 848 //((Ptarget -Qmomentum).mag2() < M0target2 )) );
>> 849 projectile->SetStatus(1*projectile->GetStatus()); // VU 10.04.2012
>> 850 }
>> 851 else
>> 852 { // -------------- Target diffraction ----------------
1038 853
1039 if ( isProjectileDiffraction ) { // projec << 854 //G4cout<<"Target diffraction"<<G4endl;
1040 projectile->SetStatus( 0 ); << 855 do {
1041 if ( projectile->GetStatus() == 2 ) proje << 856 /*
1042 if ( target->GetStatus() == 1 && target << 857 Uzhi_projectilediffraction=0;
1043 } else { // target << 858 Uzhi_targetdiffraction=1;
1044 target->SetStatus( 0 ); << 859 Uzhi_Mx2=1.;
1045 } << 860 */
>> 861 // Generate pt
>> 862 // if (whilecount++ >= 500 && (whilecount%100)==0)
>> 863 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
>> 864 // << ", loop count/ maxPtSquare : "
>> 865 // << whilecount << " / " << maxPtSquare << G4endl;
>> 866
>> 867 // whilecount++;
>> 868 if (whilecount > 1000 )
>> 869 {
>> 870 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
>> 871 target->SetStatus(2); return false; // Ignore this interaction
>> 872 };
>> 873 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl;
>> 874 // --------------- Check that the interaction is possible -----------
>> 875 ProjMassT2=M0projectile2;
>> 876 ProjMassT =M0projectile;
>> 877
>> 878 TargMassT2=TargetDiffStateMinMass2;
>> 879 TargMassT =TargetDiffStateMinMass;
>> 880
>> 881 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
>> 882 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
>> 883 /4./S;
>> 884
>> 885 //G4cout<<"PZcms2 TrD <0 "<<PZcms2<<" return"<<G4endl;
>> 886 if(PZcms2 < 0 )
>> 887 {
>> 888 target->SetStatus(2);
>> 889 return false;
>> 890 }
>> 891 maxPtSquare=PZcms2;
>> 892
>> 893 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
>> 894
>> 895 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl;
>> 896 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
>> 897
>> 898 ProjMassT2=M0projectile2+Pt2;
>> 899 ProjMassT =std::sqrt(ProjMassT2);
>> 900
>> 901 TargMassT2=TargetDiffStateMinMass2+Pt2;
>> 902 TargMassT =std::sqrt(TargMassT2);
>> 903
>> 904 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
>> 905 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
>> 906 /4./S;
>> 907
>> 908 //G4cout<<"PZcms2 <0 "<<PZcms2<<" continue"<<G4endl;
>> 909 if(PZcms2 < 0 ) continue;
>> 910 PZcms =std::sqrt(PZcms2);
>> 911
>> 912 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
>> 913 TPlusMax=SqrtS-ProjMassT;
>> 914
>> 915 TPlusNew=ChooseP(TPlusMin, TPlusMax);
>> 916 //TPlusNew=1./sqr(1./std::sqrt(TPlusMin)-G4UniformRand()*(1./std::sqrt(TPlusMin)-1./std::sqrt(TPlusMax)));
>> 917
>> 918 //TPlusNew=TPlusMax;
>> 919
>> 920 PPlusNew=SqrtS-TPlusNew;
>> 921 Qplus=PPlusNew-Pprojectile.plus();
>> 922 PMinusNew=ProjMassT2/PPlusNew;
>> 923 Qminus=PMinusNew-Pprojectile.minus();
>> 924
>> 925 Qmomentum.setPz( (Qplus-Qminus)/2 );
>> 926 Qmomentum.setE( (Qplus+Qminus)/2 );
>> 927
>> 928 /*
>> 929 G4cout<<(Pprojectile+Qmomentum).mag()<<" "<<M0projectile<<G4endl;
>> 930 G4bool First=(Pprojectile+Qmomentum).mag2() < M0projectile2;
>> 931 G4cout<<First<<G4endl;
>> 932
>> 933 G4cout<<(Ptarget -Qmomentum).mag()<<" "<<TargetDiffStateMinMass<<" "<<TargetDiffStateMinMass2<<G4endl;
>> 934 G4bool Seco=(Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2;
>> 935 G4cout<<Seco<<G4endl;
>> 936 */
>> 937
>> 938 } while ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2);
>> 939 // Repeat the sampling because there was not any excitation
>> 940 // (((Pprojectile+Qmomentum).mag2() < M0projectile2 ) || //No without excitation
>> 941 // ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2)) );
>> 942 //G4cout<<"Go out"<<G4endl;
>> 943 target->SetStatus(1*target->GetStatus()); // VU 10.04.2012
>> 944 } // End of if(G4UniformRand() < ProbProjectileDiffraction)
>> 945 }
>> 946 else //----------- Non-diffraction process ------------
>> 947 {
1046 948
1047 return true; << 949 //G4cout<<"Non-diffraction process"<<G4endl;
>> 950 do {
>> 951 /*
>> 952 Uzhi_projectilediffraction=0;
>> 953 Uzhi_targetdiffraction=0;
>> 954 Uzhi_Mx2=1.;
>> 955 */
>> 956 // Generate pt
>> 957 // if (whilecount++ >= 500 && (whilecount%100)==0)
>> 958 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
>> 959 // << ", loop count/ maxPtSquare : "
>> 960 // << whilecount << " / " << maxPtSquare << G4endl;
>> 961
>> 962 // whilecount++;
>> 963 if (whilecount > 1000 )
>> 964 {
>> 965 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
>> 966 target->SetStatus(2); return false; // Ignore this interaction
>> 967 };
>> 968 // --------------- Check that the interaction is possible -----------
>> 969 ProjMassT2=ProjectileNonDiffStateMinMass2;
>> 970 ProjMassT =ProjectileNonDiffStateMinMass;
>> 971
>> 972 TargMassT2=TargetNonDiffStateMinMass2;
>> 973 TargMassT =TargetNonDiffStateMinMass;
>> 974
>> 975 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
>> 976 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
>> 977 /4./S;
>> 978
>> 979 if(PZcms2 < 0 )
>> 980 {
>> 981 target->SetStatus(2);
>> 982 return false;
>> 983 }
>> 984 maxPtSquare=PZcms2;
>> 985 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
>> 986 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
>> 987
>> 988 ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
>> 989 ProjMassT =std::sqrt(ProjMassT2);
>> 990
>> 991 TargMassT2=TargetNonDiffStateMinMass2+Pt2;
>> 992 TargMassT =std::sqrt(TargMassT2);
>> 993
>> 994 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
>> 995 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
>> 996 /4./S;
>> 997 //G4cout<<"PZcms2 ND"<<PZcms2<<G4endl;
>> 998
>> 999 if(PZcms2 < 0 ) continue;
>> 1000 PZcms =std::sqrt(PZcms2);
>> 1001
>> 1002 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
>> 1003 PMinusMax=SqrtS-TargMassT;
>> 1004
>> 1005 if(G4UniformRand() < ProbLogDistr) // Uzhi 25.04.2012
>> 1006 { PMinusNew=ChooseP(PMinusMin, PMinusMax);} // 12.06.11
>> 1007 else {PMinusNew=(PMinusMax-PMinusMin)*G4UniformRand() + PMinusMin;}
>> 1008 Qminus=PMinusNew-Pprojectile.minus();
>> 1009
>> 1010 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
>> 1011 TPlusMax=SqrtS-PMinusNew; // Why is it closed???
>> 1012 // TPlusMax=SqrtS-ProjMassT;
>> 1013
>> 1014 if(G4UniformRand() < 0.5) //ProbLogDistr) // Uzhi 29.05.2012 0.5)
>> 1015 { TPlusNew=ChooseP(TPlusMin, TPlusMax);} // 12.06.11
>> 1016 else {TPlusNew=(TPlusMax-TPlusMin)*G4UniformRand() +TPlusMin;}
>> 1017
>> 1018 Qplus=-(TPlusNew-Ptarget.plus());
>> 1019
>> 1020 Qmomentum.setPz( (Qplus-Qminus)/2 );
>> 1021 Qmomentum.setE( (Qplus+Qminus)/2 );
>> 1022 /*
>> 1023 G4cout<<(Pprojectile+Qmomentum).mag2()<<" "<<ProjectileNonDiffStateMinMass2<<G4endl;
>> 1024 G4cout<<(Ptarget -Qmomentum).mag2()<<" "<<TargetNonDiffStateMinMass2<<G4endl;
>> 1025 G4int Uzhi; G4cin>>Uzhi;
>> 1026 */
>> 1027 } while (
>> 1028 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction
>> 1029 ((Ptarget -Qmomentum).mag2() < TargetNonDiffStateMinMass2 ));
>> 1030
>> 1031 projectile->SetStatus(0*projectile->GetStatus()); // VU 10.04.2012
>> 1032 target->SetStatus(0*target->GetStatus()); // VU 10.04.2012
>> 1033 }
>> 1034
>> 1035 Pprojectile += Qmomentum;
>> 1036 Ptarget -= Qmomentum;
>> 1037
>> 1038 //G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl;
>> 1039
>> 1040 // Transform back and update SplitableHadron Participant.
>> 1041 Pprojectile.transform(toLab);
>> 1042 Ptarget.transform(toLab);
>> 1043
>> 1044 // Calculation of the creation time ---------------------
>> 1045 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
>> 1046 projectile->SetPosition(target->GetPosition());
>> 1047 // Creation time and position of target nucleon were determined at
>> 1048 // ReggeonCascade() of G4FTFModel
>> 1049 // ------------------------------------------------------
>> 1050 /*
>> 1051 if(Uzhi_projectilediffraction != 0)
>> 1052 {Uzhi_Mx2=Pprojectile.mag2(); Uzhi_modT=(target->Get4Momentum()-Ptarget).mag2();}
>> 1053
>> 1054 if(Uzhi_targetdiffraction != 0)
>> 1055 {Uzhi_Mx2=Ptarget.mag2(); Uzhi_modT=(projectile->Get4Momentum()-Pprojectile).mag2();}
>> 1056
>> 1057 if(Uzhi_QE!= 0)
>> 1058 {
>> 1059 Uzhi_projectilediffraction=0;
>> 1060 Uzhi_targetdiffraction =0;
>> 1061 Uzhi_Mx2 =1.;
1048 } 1062 }
1049 << 1063 */
1050 //------------------------------------------- << 1064 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl;
>> 1065 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl;
>> 1066 projectile->Set4Momentum(Pprojectile);
>> 1067 target->Set4Momentum(Ptarget);
1051 1068
1052 G4bool G4DiffractiveExcitation:: << 1069 projectile->IncrementCollisionCount(1);
1053 ExciteParticipants_doNonDiffraction( G4VSplit << 1070 target->IncrementCollisionCount(1);
1054 G4VSplit <<
1055 G4FTFPar <<
1056 G4Diffra <<
1057 // Third of the three utility methods used <<
1058 // it does the sampling for the non-diffrac <<
1059 <<
1060 #ifdef debugFTFexcitation <<
1061 G4cout << "Non-diffraction process" << G4en <<
1062 #endif <<
1063 <<
1064 // Check whether the interaction is possibl <<
1065 common.ProjMassT2 = common.ProjectileNonDif <<
1066 common.ProjMassT = common.ProjectileNonDif <<
1067 common.TargMassT2 = common.TargetNonDiffSta <<
1068 common.TargMassT = common.TargetNonDiffSta <<
1069 if ( common.SqrtS < common.ProjMassT + comm <<
1070 <<
1071 common.PZcms2 = ( sqr( common.S ) + sqr( co <<
1072 - 2.0 * ( common.S * ( comm <<
1073 + common.ProjMassT2 * <<
1074 common.maxPtSquare = common.PZcms2; <<
1075 <<
1076 G4int whilecount = 0; <<
1077 do { // Generate pt and masses <<
1078 <<
1079 whilecount++; <<
1080 if ( whilecount > 1000 ) { <<
1081 common.Qmomentum = G4LorentzVector( 0.0 <<
1082 return false; // Ignore this interacti <<
1083 }; <<
1084 <<
1085 common.Qmomentum = G4LorentzVector( Gauss <<
1086 <<
1087 common.Pt2 = G4ThreeVector( common.Qmomen <<
1088 common.ProjMassT2 = common.ProjectileNonD <<
1089 common.ProjMassT = std::sqrt( common.Pro <<
1090 common.TargMassT2 = common.TargetNonDiffS <<
1091 common.TargMassT = std::sqrt( common.Tar <<
1092 if ( common.SqrtS < common.ProjMassT + co <<
1093 <<
1094 common.PZcms2 =( sqr( common.S ) + sqr( c <<
1095 - 2.0 * ( common.S * ( c <<
1096 + common.ProjM <<
1097 if ( common.PZcms2 < 0.0 ) continue; <<
1098 <<
1099 common.PZcms = std::sqrt( common.PZcms2 ) <<
1100 common.PMinusMin = std::sqrt( common.Proj <<
1101 common.PMinusMax = common.SqrtS - common. <<
1102 common.TPlusMin = std::sqrt( common.TargM <<
1103 common.TPlusMax = common.SqrtS - common.P <<
1104 <<
1105 if ( G4UniformRand() <= 0.5 ) { // Rando <<
1106 if ( G4UniformRand() < theParameters->G <<
1107 common.PMinusNew = ChooseP( common.PM <<
1108 } else { <<
1109 common.PMinusNew = ( common.PMinusMax <<
1110 } <<
1111 if ( G4UniformRand() < theParameters->G <<
1112 common.TPlusNew = ChooseP( common.TPl <<
1113 } else { <<
1114 common.TPlusNew = ( common.TPlusMax - <<
1115 } <<
1116 } else { <<
1117 if ( G4UniformRand() < theParameters->G <<
1118 common.TPlusNew = ChooseP( common.TPl <<
1119 } else { <<
1120 common.TPlusNew = ( common.TPlusMax - <<
1121 } <<
1122 if ( G4UniformRand() < theParameters->G <<
1123 common.PMinusNew = ChooseP( common.PM <<
1124 } else { <<
1125 common.PMinusNew = ( common.PMinusMax <<
1126 } <<
1127 } <<
1128 1071
1129 common.Qminus = common.PMinusNew - common << 1072 return true;
1130 common.Qplus = -( common.TPlusNew - commo <<
1131 common.Qmomentum.setPz( (common.Qplus - c <<
1132 common.Qmomentum.setE( (common.Qplus + c <<
1133 #ifdef debugFTFexcitation <<
1134 G4cout <<"Sampled: Mpr, MdifPr, Mtr, Mdif <<
1135 << ( common.Pprojectile + common.Q <<
1136 << common.ProjectileNonDiffStateMi <<
1137 << ( common.Ptarget - common.Qmome <<
1138 << common.TargetNonDiffStateMinMas <<
1139 #endif <<
1140 <<
1141 } while ( ( common.Pprojectile + common.Qmo <<
1142 common.ProjectileNonDiffStateMinM <<
1143 ( common.Ptarget - common.Qmoment <<
1144 common.TargetNonDiffStateMinMass2 <<
1145 <<
1146 projectile->SetStatus( 0 ); <<
1147 target->SetStatus( 0 ); <<
1148 return true; <<
1149 } 1073 }
1150 1074
>> 1075 // ---------------------------------------------------------------------
>> 1076 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron,
>> 1077 G4bool isProjectile,
>> 1078 G4ExcitedString * &FirstString,
>> 1079 G4ExcitedString * &SecondString,
>> 1080 G4FTFParameters *theParameters) const
>> 1081 {
>> 1082 /*
>> 1083 G4cout<<"Create Strings SplitUp "<<hadron<<G4endl;
>> 1084 G4cout<<"Defin "<<hadron->GetDefinition()<<G4endl;
>> 1085 G4cout<<"Defin "<<hadron->GetDefinition()->GetPDGEncoding()<<G4endl;
>> 1086 */
>> 1087 hadron->SplitUp();
>> 1088 G4Parton *start= hadron->GetNextParton();
>> 1089
>> 1090 if ( start==NULL)
>> 1091 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
>> 1092 FirstString=0; SecondString=0;
>> 1093 return;
>> 1094 }
>> 1095 G4Parton *end = hadron->GetNextParton();
>> 1096 if ( end==NULL)
>> 1097 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
>> 1098 FirstString=0; SecondString=0;
>> 1099 return;
>> 1100 }
>> 1101 //G4cout<<start<<" "<<start->GetPDGcode()<<" "<<end<<" "<<end->GetPDGcode()<<G4endl;
>> 1102 //G4cout<<"Create string "<<start->GetPDGcode()<<" "<<end->GetPDGcode()<<G4endl;
>> 1103 G4LorentzVector Phadron=hadron->Get4Momentum();
>> 1104 //G4cout<<"String mom "<<Phadron<<G4endl;
>> 1105 G4LorentzVector Pstart(0.,0.,0.,0.);
>> 1106 G4LorentzVector Pend(0.,0.,0.,0.);
>> 1107 G4LorentzVector Pkink(0.,0.,0.,0.);
>> 1108 G4LorentzVector PkinkQ1(0.,0.,0.,0.);
>> 1109 G4LorentzVector PkinkQ2(0.,0.,0.,0.);
>> 1110
>> 1111 G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding());
>> 1112 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding());
>> 1113 //G4cout<<"PDGcode_startQ "<<PDGcode_startQ<<" PDGcode_endQ "<<PDGcode_endQ <<G4endl;
>> 1114
>> 1115 //--------------------------------------------------------------------------------
>> 1116 G4double Wmin(0.);
>> 1117 if(isProjectile)
>> 1118 {
>> 1119 Wmin=theParameters->GetProjMinDiffMass();
>> 1120 } else
>> 1121 {
>> 1122 Wmin=theParameters->GetTarMinDiffMass();
>> 1123 } // end of if(isProjectile)
>> 1124
>> 1125 G4double W = hadron->Get4Momentum().mag();
>> 1126 //G4cout<<"Wmin W "<<Wmin<<" "<<W<<G4endl;
>> 1127 G4double W2=W*W;
>> 1128
>> 1129 G4double Pt(0.), x1(0.), x3(0.); // x2(0.),
>> 1130 G4bool Kink=false;
>> 1131
>> 1132 if(!(((start->GetDefinition()->GetParticleSubType() == "di_quark") &&
>> 1133 ( end->GetDefinition()->GetParticleSubType() == "di_quark") ) ||
>> 1134 ((start->GetDefinition()->GetParticleSubType() == "quark") &&
>> 1135 ( end->GetDefinition()->GetParticleSubType() == "quark") )))
>> 1136 { // Kinky strings are allowed only for qq-q strings
>> 1137 // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
>> 1138 // according to the analysis of Pbar P interactions
>> 1139 //G4cout<<G4endl<<"Check for Kink!##############"<<G4endl<<G4endl;
>> 1140 if(W > Wmin)
>> 1141 { // Kink is possible
>> 1142 if(hadron->GetStatus() == 0) // VU 10.04.2012
>> 1143 {
>> 1144 G4double Pt2kink=theParameters->GetPt2Kink(); // For non-diffractive
>> 1145 Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
>> 1146 }
>> 1147 else
>> 1148 {Pt=0.;} // For diffractive
1151 1149
1152 //=========================================== << 1150 if(Pt > 500.*MeV)
1153 << 1151 {
1154 void G4DiffractiveExcitation::CreateStrings( << 1152 G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
1155 << 1153 G4double Y=Ymax*(1.- 2.*G4UniformRand());
1156 << 1154
1157 << 1155 x1=1.-Pt/W*std::exp( Y);
1158 << 1156 x3=1.-Pt/W*std::exp(-Y);
1159 << 1157 // x2=2.-x1-x3;
1160 //G4cout << "Create Strings SplitUp " << ha << 1158
1161 // << "Defin " << hadron->GetDefiniti << 1159 G4double Mass_startQ = 650.*MeV;
1162 // << "Defin " << hadron->GetDefiniti << 1160 if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV;
1163 << 1161 if(PDGcode_startQ == 3) Mass_startQ = 500.*MeV;
1164 G4bool HadronIsString = hadron->IsSplit(); << 1162 if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV;
1165 if( ! HadronIsString ) hadron->SplitUp(); << 1163
1166 << 1164 G4double Mass_endQ = 650.*MeV;
1167 G4Parton* start = hadron->GetNextParton(); << 1165 if(PDGcode_endQ < 3) Mass_endQ = 325.*MeV;
1168 if ( start == nullptr ) { << 1166 if(PDGcode_endQ == 3) Mass_endQ = 500.*MeV;
1169 G4cout << " G4FTFModel::String() Error: N << 1167 if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV;
1170 FirstString = 0; SecondString = 0; << 1168
1171 return; << 1169 G4double P2_1=W2*x1*x1/4.-Mass_endQ *Mass_endQ;
1172 } << 1170 G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ;
1173 << 1171
1174 G4Parton* end = hadron->GetNextParton(); << 1172 G4double P2_2=sqr((2.-x1-x3)*W/2.);
1175 if ( end == nullptr ) { << 1173
1176 G4cout << " G4FTFModel::String() Error: N << 1174 if((P2_1 <= 0.) || (P2_3 <= 0.))
1177 FirstString = 0; SecondString = 0; << 1175 { Kink=false;}
1178 return; << 1176 else
1179 } << 1177 {
1180 << 1178 G4double P_1=std::sqrt(P2_1);
1181 //G4cout << start << " " << start->GetPDGco << 1179 G4double P_2=std::sqrt(P2_2);
1182 // << G4endl << 1180 G4double P_3=std::sqrt(P2_3);
1183 // << "Create string " << start->GetP << 1181
1184 << 1182 G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
1185 if ( HadronIsString ) { << 1183 G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
1186 if ( isProjectile ) { << 1184 // Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09
1187 FirstString = new G4ExcitedString( end, << 1185
1188 } else { << 1186 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.))
1189 FirstString = new G4ExcitedString( end, << 1187 { Kink=false;}
1190 } << 1188 else
1191 FirstString->SetTimeOfCreation( hadron->G << 1189 {
1192 FirstString->SetPosition( hadron->GetPosi << 1190 Kink=true;
1193 SecondString = 0; << 1191 Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09
1194 return; << 1192 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13);
1195 } << 1193 Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1);
>> 1194 Pkink.setPx( Pt); Pkink.setPy( 0.); Pkink.setPz( P_2*CosT12);
>> 1195 Pstart.setE(x3*W/2.);
>> 1196 Pkink.setE(Pkink.vect().mag());
>> 1197 Pend.setE(x1*W/2.);
>> 1198
>> 1199 G4double XkQ=GetQuarkFractionOfKink(0.,1.);
>> 1200 if(Pkink.getZ() > 0.)
>> 1201 {
>> 1202 if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;}
>> 1203 } else {
>> 1204 if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;}
>> 1205 }
>> 1206
>> 1207 PkinkQ2=Pkink - PkinkQ1;
>> 1208 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------
>> 1209
>> 1210 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
>> 1211 std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
>> 1212 G4double Psi=std::acos(Cos2Psi);
>> 1213
>> 1214 G4LorentzRotation Rotate;
>> 1215 if(isProjectile) {Rotate.rotateY(Psi);}
>> 1216 else {Rotate.rotateY(pi-Psi);}
>> 1217 Rotate.rotateZ(twopi*G4UniformRand());
>> 1218
>> 1219 Pstart*=Rotate;
>> 1220 Pkink*=Rotate;
>> 1221 PkinkQ1*=Rotate;
>> 1222 PkinkQ2*=Rotate;
>> 1223 Pend*=Rotate;
>> 1224
>> 1225 }
>> 1226 } // end of if((P2_1 < 0.) || (P2_3 <0.))
>> 1227 } // end of if(Pt > 500.*MeV)
>> 1228 } // end of if(W > Wmin) Check for a kink
>> 1229 } // end of qq-q string selection
>> 1230 //--------------------------------------------------------------------------------
>> 1231 /*
>> 1232 G4cout<<"Kink "<<Kink<<" "
>> 1233 <<start->GetDefinition()->GetParticleSubType()<<" "
>> 1234 << end->GetDefinition()->GetParticleSubType() <<G4endl;
>> 1235
>> 1236 G4cout<<"Kink "<<Kink<<" "
>> 1237 <<start->GetDefinition()->GetPDGEncoding()<<" "
>> 1238 << end->GetDefinition()->GetPDGEncoding() <<G4endl;
>> 1239 G4int Uzhi; G4cin>>Uzhi;
>> 1240 */
>> 1241
>> 1242 if(Kink)
>> 1243 { // Kink is possible
>> 1244 //G4cout<<"Kink is sampled!"<<G4endl;
>> 1245 std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp =
>> 1246 theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
>> 1247
>> 1248 G4int QuarkInGluon(1); G4double Ksi=G4UniformRand();
>> 1249 for(unsigned int Iq=0; Iq <3; Iq++)
>> 1250 {
>> 1251 //G4cout<<"Iq "<<Iq<<G4endl;
>> 1252
>> 1253 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
>> 1254
>> 1255 //G4cout<<"Last Iq "<<QuarkInGluon<<G4endl;
>> 1256 G4Parton * Gquark = new G4Parton(QuarkInGluon);
>> 1257 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
>> 1258 //G4cout<<"Lorentz "<<G4endl;
>> 1259
>> 1260 //-------------------------------------------------------------------------------
>> 1261 G4LorentzRotation toCMS(-1*Phadron.boostVector());
>> 1262
>> 1263 G4LorentzRotation toLab(toCMS.inverse());
>> 1264 //G4cout<<"Pstart "<<Pstart<<G4endl;
>> 1265 //G4cout<<"Pend "<<Pend<<G4endl;
>> 1266 Pstart.transform(toLab); start->Set4Momentum(Pstart);
>> 1267 PkinkQ1.transform(toLab);
>> 1268 PkinkQ2.transform(toLab);
>> 1269 Pend.transform(toLab); end->Set4Momentum(Pend);
>> 1270 //G4cout<<"Pstart "<<Pstart<<G4endl;
>> 1271 //G4cout<<"Pend "<<Pend<<G4endl;
>> 1272 //G4cout<<"Defin "<<hadron->GetDefinition()<<G4endl;
>> 1273 //G4cout<<"Defin "<<hadron->GetDefinition()->GetPDGEncoding()<<G4endl;
>> 1274
>> 1275 // G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding());
>> 1276 G4int absPDGcode=1500; // 23 Dec
>> 1277 if((start->GetDefinition()->GetParticleSubType() == "quark") &&
>> 1278 ( end->GetDefinition()->GetParticleSubType() == "quark") )
>> 1279 absPDGcode=110;
>> 1280
>> 1281 //G4cout<<"absPDGcode "<<absPDGcode<<G4endl;
>> 1282
>> 1283 if(absPDGcode < 1000)
>> 1284 { // meson
>> 1285 if ( isProjectile )
>> 1286 { // Projectile
>> 1287 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end
>> 1288 { // Quark on the end
>> 1289 FirstString = new G4ExcitedString(end ,Ganti_quark, +1);
>> 1290 SecondString= new G4ExcitedString(Gquark,start ,+1);
>> 1291 Ganti_quark->Set4Momentum(PkinkQ1);
>> 1292 Gquark->Set4Momentum(PkinkQ2);
>> 1293
>> 1294 } else
>> 1295 { // Anti_Quark on the end
>> 1296 FirstString = new G4ExcitedString(end ,Gquark, +1);
>> 1297 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
>> 1298 Gquark->Set4Momentum(PkinkQ1);
>> 1299 Ganti_quark->Set4Momentum(PkinkQ2);
>> 1300
>> 1301 } // end of if(end->GetPDGcode() > 0)
>> 1302 } else { // Target
>> 1303 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end
>> 1304 { // Quark on the end
>> 1305 FirstString = new G4ExcitedString(Ganti_quark,end ,-1);
>> 1306 SecondString= new G4ExcitedString(start ,Gquark,-1);
>> 1307 Ganti_quark->Set4Momentum(PkinkQ2);
>> 1308 Gquark->Set4Momentum(PkinkQ1);
>> 1309
>> 1310 } else
>> 1311 { // Anti_Quark on the end
>> 1312 FirstString = new G4ExcitedString(Gquark,end ,-1);
>> 1313 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
>> 1314 Gquark->Set4Momentum(PkinkQ2);
>> 1315 Ganti_quark->Set4Momentum(PkinkQ1);
>> 1316
>> 1317 } // end of if(end->GetPDGcode() > 0)
>> 1318 } // end of if ( isProjectile )
>> 1319 } else // if(absPDGCode < 1000)
>> 1320 { // Baryon/AntiBaryon
>> 1321 if ( isProjectile )
>> 1322 { // Projectile
>> 1323 if((end->GetDefinition()->GetParticleType() == "diquarks") &&
>> 1324 (end->GetDefinition()->GetPDGEncoding() > 0 ) )
>> 1325 { // DiQuark on the end
>> 1326 FirstString = new G4ExcitedString(end ,Gquark, +1);
>> 1327 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
>> 1328 Gquark->Set4Momentum(PkinkQ1);
>> 1329 Ganti_quark->Set4Momentum(PkinkQ2);
>> 1330
>> 1331 } else
>> 1332 { // Anti_DiQuark on the end or quark
>> 1333 FirstString = new G4ExcitedString(end ,Ganti_quark, +1);
>> 1334 SecondString= new G4ExcitedString(Gquark,start ,+1);
>> 1335 Ganti_quark->Set4Momentum(PkinkQ1);
>> 1336 Gquark->Set4Momentum(PkinkQ2);
>> 1337
>> 1338 } // end of if(end->GetPDGcode() > 0)
>> 1339 } else { // Target
>> 1340
>> 1341 if((end->GetDefinition()->GetParticleType() == "diquarks") &&
>> 1342 (end->GetDefinition()->GetPDGEncoding() > 0 ) )
>> 1343 { // DiQuark on the end
>> 1344 FirstString = new G4ExcitedString(Gquark,end ,-1);
>> 1345
>> 1346 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
>> 1347 Gquark->Set4Momentum(PkinkQ1);
>> 1348 Ganti_quark->Set4Momentum(PkinkQ2);
>> 1349
>> 1350 } else
>> 1351 { // Anti_DiQuark on the end or Q
>> 1352 FirstString = new G4ExcitedString(Ganti_quark,end ,-1);
>> 1353 SecondString= new G4ExcitedString(start ,Gquark,-1);
>> 1354 Gquark->Set4Momentum(PkinkQ2);
>> 1355 Ganti_quark->Set4Momentum(PkinkQ1);
>> 1356
>> 1357 } // end of if(end->GetPDGcode() > 0)
>> 1358 } // end of if ( isProjectile )
>> 1359 } // end of if(absPDGcode < 1000)
>> 1360
>> 1361 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
>> 1362 FirstString->SetPosition(hadron->GetPosition());
>> 1363
>> 1364 SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation());
>> 1365 SecondString->SetPosition(hadron->GetPosition());
>> 1366
>> 1367 // -------------------------------------------------------------------------
>> 1368 } else // End of kink is possible
>> 1369 { // Kink is impossible
>> 1370 //G4cout<<start<<" "<<start->GetPDGcode()<<" "<<end<<" "<<end->GetPDGcode()<<G4endl;
>> 1371 if ( isProjectile )
>> 1372 {
>> 1373 FirstString= new G4ExcitedString(end,start, +1);
>> 1374 } else {
>> 1375 FirstString= new G4ExcitedString(start,end, -1);
>> 1376 }
1196 1377
1197 G4LorentzVector Phadron = hadron->Get4Momen << 1378 SecondString=0;
1198 //G4cout << "String mom " << Phadron << G4e <<
1199 G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 <<
1200 G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 ); <<
1201 G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 ) <<
1202 G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 <<
1203 G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 <<
1204 <<
1205 G4int PDGcode_startQ = std::abs( start->Get <<
1206 G4int PDGcode_endQ = std::abs( end->Get <<
1207 //G4cout << "PDGcode_startQ " << PDGcode_st <<
1208 <<
1209 G4double Wmin( 0.0 ); <<
1210 if ( isProjectile ) { <<
1211 Wmin = theParameters->GetProjMinDiffMass( <<
1212 } else { <<
1213 Wmin = theParameters->GetTarMinDiffMass() <<
1214 } <<
1215 1379
1216 G4double W = hadron->Get4Momentum().mag(); << 1380 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
1217 G4double W2 = W*W; << 1381 FirstString->SetPosition(hadron->GetPosition());
1218 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); <<
1219 G4bool Kink = false; <<
1220 <<
1221 if ( ! ( ( start->GetDefinition()->GetParti <<
1222 end->GetDefinition()->GetParti <<
1223 ( start->GetDefinition()->GetParti <<
1224 end->GetDefinition()->GetParti <<
1225 // Kinky strings are allowed only for qq- <<
1226 // Kinky strings are impossible for other <<
1227 // according to the analysis of Pbar P in <<
1228 <<
1229 if ( W > Wmin ) { // Kink is possible <<
1230 if ( hadron->GetStatus() == 0 ) { <<
1231 G4double Pt2kink = theParameters->Get <<
1232 if ( Pt2kink ) { <<
1233 Pt = std::sqrt( Pt2kink * ( G4Pow:: <<
1234 } else { <<
1235 Pt = 0.0; <<
1236 } <<
1237 } else { <<
1238 Pt = 0.0; <<
1239 } <<
1240 1382
1241 if ( Pt > 500.0*MeV ) { << 1383 // momenta of string ends
1242 G4double Ymax = G4Log( W/2.0/Pt + std << 1384 //
1243 G4double Y = Ymax*( 1.0 - 2.0*G4Unifo << 1385 G4double Momentum=hadron->Get4Momentum().vect().mag();
1244 x1 = 1.0 - Pt/W * G4Exp( Y ); << 1386 G4double Plus=hadron->Get4Momentum().e() + Momentum;
1245 x3 = 1.0 - Pt/W * G4Exp(-Y ); << 1387 G4double Minus=hadron->Get4Momentum().e() - Momentum;
1246 //x2 = 2.0 - x1 - x3; << 1388
1247 << 1389 G4ThreeVector tmp;
1248 G4double Mass_startQ = 650.0*MeV; << 1390 if(Momentum > 0.)
1249 if ( PDGcode_startQ < 3 ) Mass_start << 1391 {
1250 if ( PDGcode_startQ == 3 ) Mass_start << 1392 tmp.set(hadron->Get4Momentum().px(),
1251 if ( PDGcode_startQ == 4 ) Mass_start << 1393 hadron->Get4Momentum().py(),
1252 if ( PDGcode_startQ == 5 ) Mass_start << 1394 hadron->Get4Momentum().pz());
1253 G4double Mass_endQ = 650.0*MeV; << 1395 tmp/=Momentum;
1254 if ( PDGcode_endQ < 3 ) Mass_endQ = << 1396 }
1255 if ( PDGcode_endQ == 3 ) Mass_endQ = << 1397 else
1256 if ( PDGcode_endQ == 4 ) Mass_endQ = << 1398 {
1257 if ( PDGcode_endQ == 5 ) Mass_endQ = << 1399 tmp.set(0.,0.,1.);
1258 <<
1259 G4double P2_1 = W2*x1*x1/4.0 - Mass_e <<
1260 G4double P2_3 = W2*x3*x3/4.0 - Mass_s <<
1261 G4double P2_2 = sqr( (2.0 - x1 - x3)* <<
1262 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) { <<
1263 Kink = false; <<
1264 } else { <<
1265 G4double P_1 = std::sqrt( P2_1 ); <<
1266 G4double P_2 = std::sqrt( P2_2 ); <<
1267 G4double P_3 = std::sqrt( P2_3 ); <<
1268 G4double CosT12 = ( P2_3 - P2_1 - P <<
1269 G4double CosT13 = ( P2_2 - P2_1 - P <<
1270 <<
1271 if ( std::abs( CosT12 ) > 1.0 || <<
1272 Kink = false; <<
1273 } else { <<
1274 Kink = true; <<
1275 Pt = P_2 * std::sqrt( 1.0 - CosT1 <<
1276 Pstart.setPx( -Pt ); Pstart.setPy <<
1277 Pend.setPx( 0.0 ); Pend.setPy( <<
1278 Pkink.setPx( Pt ); Pkink.setPy( <<
1279 Pstart.setE( x3*W/2.0 ); <<
1280 Pkink.setE( Pkink.vect().mag() ); <<
1281 Pend.setE( x1*W/2.0 ); <<
1282 <<
1283 G4double XkQ = GetQuarkFractionOf <<
1284 if ( Pkink.getZ() > 0.0 ) { <<
1285 if ( XkQ > 0.5 ) { <<
1286 PkinkQ1 = XkQ*Pkink; <<
1287 } else { <<
1288 PkinkQ1 = (1.0 - XkQ)*Pkink; <<
1289 } <<
1290 } else { <<
1291 if ( XkQ > 0.5 ) { <<
1292 PkinkQ1 = (1.0 - XkQ)*Pkink; <<
1293 } else { <<
1294 PkinkQ1 = XkQ*Pkink; <<
1295 } <<
1296 } <<
1297 <<
1298 PkinkQ2 = Pkink - PkinkQ1; <<
1299 // Minimizing Pt1^2+Pt3^2 <<
1300 G4double Cos2Psi = ( sqr(x1) - sq <<
1301 std::sqrt( sqr <<
1302 G4double Psi = std::acos( Cos2Psi <<
1303 <<
1304 G4LorentzRotation Rotate; <<
1305 if ( isProjectile ) { <<
1306 Rotate.rotateY( Psi ); <<
1307 } else { <<
1308 Rotate.rotateY( pi + Psi ); <<
1309 } <<
1310 Rotate.rotateZ( twopi * G4Uniform <<
1311 Pstart *= Rotate; <<
1312 Pkink *= Rotate; <<
1313 PkinkQ1 *= Rotate; <<
1314 PkinkQ2 *= Rotate; <<
1315 Pend *= Rotate; <<
1316 } 1400 }
1317 } // End of if ( P2_1 <= 0.0 || P2 <<
1318 } // End of if ( Pt > 500.0*MeV ) <<
1319 } // End of if ( W > Wmin ) : check for a <<
1320 } // end of qq-q string selection <<
1321 <<
1322 if ( Kink ) { // Kink is possible <<
1323 <<
1324 //G4cout << "Kink is sampled!" << G4endl; <<
1325 std::vector< G4double > QuarkProbabilitie <<
1326 theParameters->GetQuarkProbabilitiesA <<
1327 <<
1328 G4int QuarkInGluon( 1 ); G4double Ksi = G <<
1329 for ( unsigned int Iq = 0; Iq < 3; Iq++ ) <<
1330 //G4cout << "Iq " << Iq << G4endl; <<
1331 if ( Ksi > QuarkProbabilitiesAtGluonSpl <<
1332 } <<
1333 //G4cout << "Last Iq " << QuarkInGluon << <<
1334 G4Parton* Gquark = new G4Parton( QuarkInG <<
1335 G4Parton* Ganti_quark = new G4Parton( -Qu <<
1336 //G4cout << "Lorentz " << G4endl; <<
1337 <<
1338 G4LorentzRotation toCMS( -1 * Phadron.boo <<
1339 G4LorentzRotation toLab( toCMS.inverse() <<
1340 //G4cout << "Pstart " << Pstart << G4endl <<
1341 //G4cout << "Pend " << Pend << G4endl; <<
1342 //G4cout << "Kink1 " <<PkinkQ1<<G4endl; <<
1343 //G4cout << "Kink2 " <<PkinkQ2<<G4endl; <<
1344 //G4cout << "Pstart " << Pstart << G4endl <<
1345 <<
1346 Pstart.transform( toLab ); start->Set4Mo <<
1347 PkinkQ1.transform( toLab ); <<
1348 PkinkQ2.transform( toLab ); <<
1349 Pend.transform( toLab ); end->Set4Mome <<
1350 //G4cout << "Pstart " << Pstart << G4endl <<
1351 //G4cout << "Pend " << Pend << G4endl; <<
1352 //G4cout << "Defin " << hadron->GetDefini <<
1353 //G4cout << "Defin " << hadron->GetDefini <<
1354 <<
1355 //G4int absPDGcode = std::abs( hadron->Ge <<
1356 G4int absPDGcode = 1500; <<
1357 if ( start->GetDefinition()->GetParticleS <<
1358 end->GetDefinition()->GetParticleSub <<
1359 absPDGcode = 110; <<
1360 } <<
1361 //G4cout << "absPDGcode " << absPDGcode < <<
1362 1401
1363 if ( absPDGcode < 1000 ) { // meson << 1402 G4LorentzVector Pstart1(tmp,0.);
1364 if ( isProjectile ) { // Projectile << 1403 G4LorentzVector Pend1(tmp,0.);
1365 if ( end->GetDefinition()->GetPDGEnco <<
1366 FirstString = new G4ExcitedString( <<
1367 SecondString = new G4ExcitedString( <<
1368 Ganti_quark->Set4Momentum( PkinkQ1 <<
1369 Gquark->Set4Momentum( PkinkQ2 ); <<
1370 } else { // Anti_Quark on the end <<
1371 FirstString = new G4ExcitedString( <<
1372 SecondString = new G4ExcitedString( <<
1373 Gquark->Set4Momentum( PkinkQ1 ); <<
1374 Ganti_quark->Set4Momentum( PkinkQ2 <<
1375 } <<
1376 } else { // Target <<
1377 if ( end->GetDefinition()->GetPDGEnco <<
1378 FirstString = new G4ExcitedString( <<
1379 SecondString = new G4ExcitedString( <<
1380 Ganti_quark->Set4Momentum( PkinkQ2 <<
1381 Gquark->Set4Momentum( PkinkQ1 ); <<
1382 } else { // Anti_Quark on the end <<
1383 FirstString = new G4ExcitedString( <<
1384 SecondString = new G4ExcitedString( <<
1385 Gquark->Set4Momentum( PkinkQ2 ); <<
1386 Ganti_quark->Set4Momentum( PkinkQ1 <<
1387 } <<
1388 } <<
1389 } else { // Baryon/AntiBaryon <<
1390 if ( isProjectile ) { // Projectile <<
1391 if ( end->GetDefinition()->GetParticl <<
1392 end->GetDefinition()->GetPDGEnco <<
1393 FirstString = new G4ExcitedString( <<
1394 SecondString = new G4ExcitedString( <<
1395 Gquark->Set4Momentum( PkinkQ1 ); <<
1396 Ganti_quark->Set4Momentum( PkinkQ2 <<
1397 } else { / <<
1398 FirstString = new G4ExcitedString( <<
1399 SecondString = new G4ExcitedString( <<
1400 Ganti_quark->Set4Momentum( PkinkQ1 <<
1401 Gquark->Set4Momentum( PkinkQ2 ); <<
1402 } <<
1403 } else { // Target <<
1404 if ( end->GetDefinition()->GetParticl <<
1405 end->GetDefinition()->GetPDGEnco <<
1406 Gquark->Set4Momentum( PkinkQ1 ); <<
1407 Ganti_quark->Set4Momentum( PkinkQ2 <<
1408 FirstString = new G4ExcitedString( <<
1409 SecondString = new G4ExcitedString( <<
1410 } else { // Anti_DiQuark on the end <<
1411 FirstString = new G4ExcitedString( <<
1412 SecondString = new G4ExcitedString( <<
1413 Gquark->Set4Momentum( PkinkQ2 ); <<
1414 Ganti_quark->Set4Momentum( PkinkQ1 <<
1415 } <<
1416 } <<
1417 } <<
1418 1404
1419 FirstString->SetTimeOfCreation( hadron->G << 1405 if(isProjectile)
1420 FirstString->SetPosition( hadron->GetPosi << 1406 {
1421 SecondString->SetTimeOfCreation( hadron-> << 1407 Pstart1*=(-1.)*Minus/2.;
1422 SecondString->SetPosition( hadron->GetPos << 1408 Pend1 *=(+1.)*Plus /2.;
1423 << 1409 }
1424 } else { // Kink is impossible << 1410 else
1425 << 1411 {
1426 if ( isProjectile ) { << 1412 Pstart1*=(+1.)*Plus/2.;
1427 FirstString = new G4ExcitedString( end, << 1413 Pend1 *=(-1.)*Minus/2.;
1428 } else { << 1414 }
1429 FirstString = new G4ExcitedString( end, <<
1430 } <<
1431 FirstString->SetTimeOfCreation( hadron->G <<
1432 FirstString->SetPosition( hadron->GetPosi <<
1433 SecondString = 0; <<
1434 // momenta of string ends <<
1435 G4LorentzVector HadronMom = hadron->Get4M <<
1436 G4LorentzVector Pstart1 = G4LorentzVector <<
1437 G4LorentzVector Pend1 = G4LorentzVector <<
1438 G4double Pz = HadronMom.pz(); <<
1439 G4double Eh = HadronMom.e(); <<
1440 G4double Pt2 = sqr( HadronMom.px() ) + sq <<
1441 G4double Mt2 = HadronMom.mt2(); <<
1442 G4double Exp = std::sqrt( sqr(Pz) + ( sqr <<
1443 G4double Pzq = Pz/2.0 - Exp; <<
1444 G4double Eq = std::sqrt( sqr(Pzq) + Pt2 <<
1445 G4double Pzqq = Pz/2.0 + Exp; <<
1446 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt <<
1447 start->Set4Momentum( Pstart1 ); <<
1448 end->Set4Momentum( Pend1 ); <<
1449 Pstart = Pstart1; Pend = Pend1; <<
1450 <<
1451 } // End of "if (Kink)" <<
1452 <<
1453 //G4cout << "Quarks in the string at creati <<
1454 // << " " << FirstString->GetLeftPart <<
1455 // << FirstString << " " << SecondStr <<
1456 <<
1457 #ifdef G4_FTFDEBUG <<
1458 G4cout << " generated string flavors <<
1459 << end->GetPDGcode() << G4endl << " <<
1460 << start->Get4Momentum() << "mass : <<
1461 << " generated string momenta: Diqua <<
1462 << end->Get4Momentum().mag() << G4en <<
1463 << Pstart + Pend << G4endl << " Orig <<
1464 << hadron->Get4Momentum() << " "<<ha <<
1465 #endif <<
1466 1415
1467 return; << 1416 Momentum=-Pstart1.mag();
1468 } << 1417 Pstart1.setT(Momentum); // It is assumed that quark has m=0.
1469 1418
>> 1419 Momentum=-Pend1.mag();
>> 1420 Pend1.setT(Momentum); // It is assumed that di-quark has m=0.
1470 1421
1471 //=========================================== << 1422 start->Set4Momentum(Pstart1);
>> 1423 end->Set4Momentum(Pend1);
>> 1424 SecondString=0;
>> 1425 } // End of kink is impossible
>> 1426
>> 1427 //G4cout<<"Quarks in the string at creation"<<FirstString->GetRightParton()->GetPDGcode()<<" "<<FirstString->GetLeftParton()->GetPDGcode()<<G4endl;
>> 1428 //G4cout<<FirstString<<" "<<SecondString<<G4endl;
>> 1429
>> 1430 #ifdef G4_FTFDEBUG
>> 1431 G4cout << " generated string flavors "
>> 1432 << start->GetPDGcode() << " / "
>> 1433 << end->GetPDGcode() << G4endl;
>> 1434 G4cout << " generated string momenta: quark "
>> 1435 << start->Get4Momentum() << "mass : "
>> 1436 <<start->Get4Momentum().mag() << G4endl;
>> 1437 G4cout << " generated string momenta: Diquark "
>> 1438 << end ->Get4Momentum()
>> 1439 << "mass : " <<end->Get4Momentum().mag()<< G4endl;
>> 1440 G4cout << " sum of ends " << Pstart+Pend << G4endl;
>> 1441 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
>> 1442 #endif
1472 1443
1473 G4double G4DiffractiveExcitation::ChooseP( G4 << 1444 return;
1474 // Choose an x between Xmin and Xmax with P << 1445
1475 // To be improved... <<
1476 G4double range = Pmax - Pmin; <<
1477 if ( Pmin <= 0.0 || range <= 0.0 ) { <<
1478 G4cout << " Pmin, range : " << Pmin << " <<
1479 throw G4HadronicException( __FILE__, __LI <<
1480 "G4Diffractive <<
1481 } <<
1482 G4double P = Pmin * G4Pow::GetInstance()->p <<
1483 //G4double P = (Pmax - Pmin) * G4UniformRan <<
1484 return P; <<
1485 } 1446 }
1486 1447
1487 1448
1488 //=========================================== << 1449 // --------- private methods ----------------------
1489 1450
1490 G4ThreeVector G4DiffractiveExcitation::Gaussi << 1451 // ---------------------------------------------------------------------
1491 // @@ this method is used in FTFModel as w << 1452 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
1492 G4double Pt2( 0.0 ); << 1453 {
1493 if ( AveragePt2 <= 0.0 ) { << 1454 // choose an x between Xmin and Xmax with P(x) ~ 1/x
1494 Pt2 = 0.0; << 1455 // to be improved...
1495 } else { << 1456
1496 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor << 1457 G4double range=Pmax-Pmin;
1497 ( G4Ex << 1458
1498 } << 1459 if ( Pmin <= 0. || range <=0. )
1499 G4double Pt = std::sqrt( Pt2 ); << 1460 {
1500 G4double phi = G4UniformRand() * twopi; << 1461 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1501 return G4ThreeVector( Pt * std::cos( phi ), << 1462 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
>> 1463 }
>> 1464
>> 1465 G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());
>> 1466 // G4double P=(Pmax-Pmin)*G4UniformRand()+Pmin;
>> 1467 return P;
1502 } 1468 }
1503 1469
1504 << 1470 // ---------------------------------------------------------------------
1505 //=========================================== << 1471 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2,
1506 << 1472 G4double maxPtSquare) const
1507 G4double G4DiffractiveExcitation::GetQuarkFra << 1473 { // @@ this method is used in FTFModel as well. Should go somewhere common!
1508 G4double z, yf; << 1474
1509 const G4int maxNumberOfLoops = 10000; << 1475 G4double Pt2(0.);
1510 G4int loopCounter = 0; << 1476 if(AveragePt2 <= 0.) {Pt2=0.;}
1511 do { << 1477 else
1512 z = zmin + G4UniformRand() * (zmax - zmin << 1478 {
1513 yf = z*z + sqr(1.0 - z); << 1479 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
1514 } while ( ( G4UniformRand() > yf ) && << 1480 (std::exp(-maxPtSquare/AveragePt2)-1.));
1515 ++loopCounter < maxNumberOfLoops << 1481 }
1516 if ( loopCounter >= maxNumberOfLoops ) { << 1482 G4double Pt=std::sqrt(Pt2);
1517 z = 0.5*(zmin + zmax); // Just something << 1483 G4double phi=G4UniformRand() * twopi;
1518 } << 1484 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
1519 return z; <<
1520 } 1485 }
1521 1486
1522 << 1487 // ---------------------------------------------------------------------
1523 //=========================================== << 1488 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
1524 << 1489 {
1525 void G4DiffractiveExcitation::UnpackMeson( co << 1490 G4double z, yf;
1526 G4int absIdPDG = std::abs( IdPDG ); << 1491 do {
1527 if ( ! ( absIdPDG == 111 || absIdPDG == 221 << 1492 z = zmin + G4UniformRand()*(zmax-zmin);
1528 absIdPDG == 441 || absIdPDG == 443 << 1493 yf = z*z +sqr(1 - z);
1529 // All other projectile mesons (including << 1494 }
1530 Q1 = absIdPDG / 100; << 1495 while (G4UniformRand() > yf);
1531 Q2 = (absIdPDG % 100) / 10; << 1496 return z;
1532 G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) <<
1533 if ( IdPDG < 0 ) anti *= -1; <<
1534 Q1 *= anti; <<
1535 Q2 *= -1 * anti; <<
1536 } else { <<
1537 if ( absIdPDG == 441 || absIdPDG == 443 ) <<
1538 Q1 = 4; Q2 = -4; <<
1539 } else if ( absIdPDG == 553 ) { <<
1540 Q1 = 5; Q2 = -5; <<
1541 } else { <<
1542 if ( G4UniformRand() < 0.5 ) { <<
1543 Q1 = 1; Q2 = -1; <<
1544 } else { <<
1545 Q1 = 2; Q2 = -2; <<
1546 } <<
1547 } 1497 }
1548 } << 1498 // ---------------------------------------------------------------------
1549 return; << 1499 void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09
1550 } << 1500 {
1551 << 1501 G4int absIdPDG = std::abs(IdPDG);
1552 << 1502 Q1= absIdPDG/ 100;
1553 //=========================================== << 1503 Q2= (absIdPDG %100)/10;
1554 << 1504
1555 void G4DiffractiveExcitation::UnpackBaryon( G << 1505 G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 );
1556 G << 1506
1557 Q1 = IdPDG / 1000; << 1507 if (IdPDG < 0 ) anti *=-1;
1558 Q2 = (IdPDG % 1000) / 100; << 1508 Q1 *= anti;
1559 Q3 = (IdPDG % 100) / 10; << 1509 Q2 *= -1 * anti;
1560 return; << 1510 return;
1561 } << 1511 }
1562 << 1512 // ---------------------------------------------------------------------
1563 << 1513 void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG,
1564 //=========================================== << 1514 G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09
>> 1515 {
>> 1516 Q1 = IdPDG / 1000;
>> 1517 Q2 = (IdPDG % 1000) / 100;
>> 1518 Q3 = (IdPDG % 100) / 10;
>> 1519 return;
>> 1520 }
>> 1521 // ---------------------------------------------------------------------
>> 1522 G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09
>> 1523 {
>> 1524 G4int TmpQ(0);
>> 1525 if( Q3 > Q2 )
>> 1526 {
>> 1527 TmpQ = Q2;
>> 1528 Q2 = Q3;
>> 1529 Q3 = TmpQ;
>> 1530 } else if( Q3 > Q1 )
>> 1531 {
>> 1532 TmpQ = Q1;
>> 1533 Q1 = Q3;
>> 1534 Q3 = TmpQ;
>> 1535 }
>> 1536
>> 1537 if( Q2 > Q1 )
>> 1538 {
>> 1539 TmpQ = Q1;
>> 1540 Q1 = Q2;
>> 1541 Q2 = TmpQ;
>> 1542 }
1565 1543
1566 G4int G4DiffractiveExcitation::NewNucleonId( << 1544 G4int NewCode = Q1*1000 + Q2* 100 + Q3* 10 + 2;
1567 // Order the three integers in such a way t << 1545 return NewCode;
1568 G4int TmpQ( 0 ); << 1546 }
1569 if ( Q3 > Q2 ) { << 1547 // ---------------------------------------------------------------------
1570 TmpQ = Q2; << 1548 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
1571 Q2 = Q3; << 1549 {
1572 Q3 = TmpQ; << 1550 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
1573 } else if ( Q3 > Q1 ) { <<
1574 TmpQ = Q1; <<
1575 Q1 = Q3; <<
1576 Q3 = TmpQ; <<
1577 } <<
1578 if ( Q2 > Q1 ) { <<
1579 TmpQ = Q1; <<
1580 Q1 = Q2; <<
1581 Q2 = TmpQ; <<
1582 } <<
1583 // By now Q1 >= Q2 >= Q3 <<
1584 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + <<
1585 return NewCode; <<
1586 } 1551 }
1587 1552
1588 1553
1589 //=========================================== << 1554 G4DiffractiveExcitation::~G4DiffractiveExcitation()
1590 << 1555 {
1591 G4DiffractiveExcitation::G4DiffractiveExcitat <<
1592 throw G4HadronicException( __FILE__, __LINE <<
1593 "G4DiffractiveEx <<
1594 } 1556 }
1595 1557
1596 1558
1597 //=========================================== << 1559 const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
1598 << 1560 {
1599 const G4DiffractiveExcitation & G4Diffractive << 1561 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
1600 throw G4HadronicException( __FILE__, __LINE << 1562 return *this;
1601 "G4DiffractiveEx <<
1602 return *this; <<
1603 } 1563 }
1604 1564
1605 1565
1606 //=========================================== << 1566 int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
1607 << 1567 {
1608 G4bool G4DiffractiveExcitation::operator==( c << 1568 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
1609 throw G4HadronicException( __FILE__, __LINE << 1569 return false;
1610 "G4DiffractiveEx <<
1611 } 1570 }
1612 1571
1613 << 1572 int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
1614 //=========================================== << 1573 {
1615 << 1574 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
1616 G4bool G4DiffractiveExcitation::operator!= ( << 1575 return true;
1617 throw G4HadronicException( __FILE__, __LINE <<
1618 "G4DiffractiveEx <<
1619 } 1576 }
1620 <<
1621 1577