| 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 // $Id: G4QGSDiffractiveExcitation.cc 69570 2013-05-08 13:23:39Z gcosmo $
27 // ------------------------------------------- 28 // ------------------------------------------------------------
28 // GEANT 4 class implemetation file 29 // GEANT 4 class implemetation file
29 // 30 //
30 // ---------------- G4QGSDiffractiveExcit 31 // ---------------- G4QGSDiffractiveExcitation --------------
31 // by Gunter Folger, October 1998. 32 // by Gunter Folger, October 1998.
32 // diffractive Excitation used by stri 33 // diffractive Excitation used by strings models
33 // Take a projectile and a target 34 // Take a projectile and a target
34 // excite the projectile and target 35 // excite the projectile and target
35 // Essential changed by V. Uzhinsky in Novemb 36 // Essential changed by V. Uzhinsky in November - December 2006
36 // in order to put it in a correspondence wit 37 // in order to put it in a correspondence with original FRITIOF
37 // model. Variant of FRITIOF with nucleon de- 38 // model. Variant of FRITIOF with nucleon de-excitation is implemented.
38 // ------------------------------------------- 39 // ---------------------------------------------------------------------
39 40
40 // Modified: 41 // Modified:
41 // 25-05-07 : G.Folger 42 // 25-05-07 : G.Folger
42 // move from management/G4DiffractiveExc 43 // move from management/G4DiffractiveExcitation to to qgsm/G4QGSDiffractiveExcitation
43 // 44 //
44 45
45 #include "globals.hh" 46 #include "globals.hh"
46 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh" 49 #include "Randomize.hh"
49 50
50 #include "G4QGSDiffractiveExcitation.hh" 51 #include "G4QGSDiffractiveExcitation.hh"
51 #include "G4LorentzRotation.hh" 52 #include "G4LorentzRotation.hh"
52 #include "G4ThreeVector.hh" 53 #include "G4ThreeVector.hh"
53 #include "G4ParticleDefinition.hh" 54 #include "G4ParticleDefinition.hh"
54 #include "G4VSplitableHadron.hh" 55 #include "G4VSplitableHadron.hh"
55 #include "G4ExcitedString.hh" 56 #include "G4ExcitedString.hh"
56 //#include "G4ios.hh" 57 //#include "G4ios.hh"
57 58
58 #include "G4Exp.hh" << 59 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation() // Uzhi
59 #include "G4Log.hh" <<
60 #include "G4Pow.hh" <<
61 <<
62 //============================================ <<
63 <<
64 //#define debugDoubleDiffraction <<
65 <<
66 //============================================ <<
67 <<
68 G4QGSDiffractiveExcitation::G4QGSDiffractiveEx <<
69 { 60 {
70 } 61 }
71 62
72 G4QGSDiffractiveExcitation::~G4QGSDiffractiveE 63 G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation()
73 { 64 {
74 } 65 }
75 66
76 67
77 G4bool G4QGSDiffractiveExcitation:: 68 G4bool G4QGSDiffractiveExcitation::
78 ExciteParticipants(G4VSplitableHadron *project << 69 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
79 { 70 {
80 #ifdef debugDoubleDiffraction <<
81 G4cout<<G4endl<<"G4QGSDiffractiveExcitatio <<
82 G4cout<<"Proj Targ "<<projectile->GetDefin <<
83 G4cout<<"Proj 4 Mom "<<projectile->Get4Mom <<
84 G4cout<<"Targ 4 Mom "<<target->Get4Momentu <<
85 #endif <<
86 <<
87 G4LorentzVector Pprojectile=projectile->Get4 <<
88 <<
89 // -------------------- Projectile parameter <<
90 G4bool PutOnMassShell=0; <<
91 <<
92 G4double M0projectile = Pprojectile.mag(); <<
93 <<
94 if(M0projectile < projectile->GetDefinition( <<
95 { <<
96 PutOnMassShell=1; <<
97 M0projectile=projectile->GetDefinition()-> <<
98 } <<
99 <<
100 // -------------------- Target parameters -- <<
101 G4LorentzVector Ptarget=target->Get4Momentum <<
102 <<
103 G4double M0target = Ptarget.mag(); <<
104 <<
105 if(M0target < target->GetDefinition()->GetPD <<
106 { <<
107 PutOnMassShell=1; <<
108 M0target=target->GetDefinition()->GetPDGMa <<
109 } <<
110 <<
111 G4LorentzVector Psum=Pprojectile+Ptarget; <<
112 G4double S=Psum.mag2(); <<
113 G4double SqrtS=std::sqrt(S); <<
114 <<
115 if (SqrtS < M0projectile + M0target) {return <<
116 <<
117 <<
118 G4double Mprojectile2 = M0projectile * M0pro <<
119 G4double Mtarget2 = M0target * M0tar <<
120 <<
121 // Transform momenta to cms and then rotate <<
122 <<
123 G4LorentzRotation toCms(-1*Psum.boostVector( <<
124 <<
125 G4LorentzVector Ptmp=toCms*Pprojectile; <<
126 <<
127 if ( Ptmp.pz() <= 0. ) <<
128 { <<
129 // "String" moving backwards in CMS, abor <<
130 //G4cout << " abort Collision!! " << G4end <<
131 return false; <<
132 } <<
133 <<
134 toCms.rotateZ(-1*Ptmp.phi()); <<
135 toCms.rotateY(-1*Ptmp.theta()); <<
136 <<
137 G4LorentzRotation toLab(toCms.inverse()); <<
138 <<
139 Pprojectile.transform(toCms); <<
140 Ptarget.transform(toCms); <<
141 <<
142 G4double PZcms2=(S*S+Mprojectile2*Mprojectil <<
143 2*S*Mprojectile2-2*S*Mtarget2-2*Mproj <<
144 <<
145 if (PZcms2 < 0) {return false;} // It can <<
146 <<
147 G4double PZcms = std::sqrt(PZcms2); <<
148 <<
149 if (PutOnMassShell) <<
150 { <<
151 if (Pprojectile.z() > 0.) <<
152 { <<
153 Pprojectile.setPz( PZcms); <<
154 Ptarget.setPz( -PZcms); <<
155 } <<
156 else <<
157 { <<
158 Pprojectile.setPz(-PZcms); <<
159 Ptarget.setPz( PZcms); <<
160 }; <<
161 <<
162 Pprojectile.setE(std::sqrt(Mprojectile2+sq <<
163 Ptarget.setE( std::sqrt( Mtarget2+sq <<
164 } <<
165 <<
166 G4double maxPtSquare = PZcms2; <<
167 <<
168 #ifdef debugDoubleDiffraction <<
169 G4cout << "Pprojectile after boost to CMS: <<
170 G4cout << "Ptarget after boost to CMS: <<
171 #endif <<
172 <<
173 G4int PrPDGcode=projectile->GetDefinition <<
174 G4int absPrPDGcode=std::abs(PrPDGcode); <<
175 G4double MinPrDiffMass(0.); <<
176 G4double AveragePt2(0.); <<
177 <<
178 if (M0projectile <= projectile->GetDefinitio <<
179 { // Normal projectile <<
180 if( absPrPDGcode > 1000 ) <<
181 { <<
182 if ( absPrPDGcode > 4000 && absPrPDGcode <<
183 { <<
184 MinPrDiffMass = projectile->GetDefinit <<
185 AveragePt2 = 0.3; <<
186 } <<
187 else <<
188 { <<
189 MinPrDiffMass = 1.16; // <<
190 AveragePt2 = 0.3; <<
191 } <<
192 } <<
193 else if( absPrPDGcode == 211 || PrPDGcode <<
194 { <<
195 MinPrDiffMass = 1.0; <<
196 AveragePt2 = 0.3; <<
197 } <<
198 else if( absPrPDGcode == 321 || absPrPDGco <<
199 { <<
200 MinPrDiffMass = 1.1; <<
201 AveragePt2 = 0.3; <<
202 } <<
203 else if( absPrPDGcode > 400 && absPrPDGcod <<
204 { <<
205 MinPrDiffMass = projectile->GetDefinitio <<
206 AveragePt2 = 0.3; <<
207 } <<
208 else <<
209 { <<
210 MinPrDiffMass = 1.16; <<
211 AveragePt2 = 0.3; <<
212 } <<
213 } <<
214 else <<
215 { // Excited projectile <<
216 MinPrDiffMass = M0projectile + 220.0*MeV; <<
217 AveragePt2 = 0.3; <<
218 } <<
219 <<
220 MinPrDiffMass = MinPrDiffMass * GeV; <<
221 AveragePt2 = AveragePt2 * GeV*GeV; <<
222 //------------------------------------------ <<
223 G4double MinTrDiffMass = 1.16*GeV; <<
224 <<
225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) { <<
226 <<
227 G4double MinPrDiffMass2 = MinPrDiffMass * Mi <<
228 G4double MinTrDiffMass2 = MinTrDiffMass * Mi <<
229 <<
230 G4double Pt2; <<
231 G4double ProjMassT2, ProjMassT; <<
232 G4double TargMassT2, TargMassT; <<
233 G4double PMinusNew, TPlusNew; <<
234 <<
235 G4LorentzVector Qmomentum; <<
236 G4double Qminus, Qplus; <<
237 <<
238 G4int whilecount=0; <<
239 do { <<
240 if (whilecount++ >= 500 && (whilecount%100 <<
241 if (whilecount > 1000 ) { <<
242 Qmomentum=G4LorentzVector(0.,0.,0.,0.); <<
243 return false; // Ignore this interactio <<
244 } <<
245 <<
246 // Generate pt <<
247 Qmomentum=G4LorentzVector(GaussianPt(Avera <<
248 <<
249 Pt2=G4ThreeVector(Qmomentum.vect()).mag2() <<
250 ProjMassT2=MinPrDiffMass2+Pt2; <<
251 ProjMassT =std::sqrt(ProjMassT2); <<
252 <<
253 TargMassT2=MinTrDiffMass2+Pt2; <<
254 TargMassT =std::sqrt(TargMassT2); <<
255 <<
256 if (SqrtS < ProjMassT + TargMassT) continu <<
257 <<
258 PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMass <<
259 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjM <<
260 if (PZcms2 < 0 ) {PZcms2=0;}; <<
261 PZcms =std::sqrt(PZcms2); <<
262 <<
263 G4double PMinusMin=std::sqrt(ProjMassT2+PZ <<
264 G4double PMinusMax=SqrtS-TargMassT; <<
265 <<
266 PMinusNew=ChooseP(PMinusMin,PMinusMax); <<
267 Qminus=PMinusNew-Pprojectile.minus(); <<
268 <<
269 G4double TPlusMin=std::sqrt(TargMassT2+PZc <<
270 G4double TPlusMax=SqrtS-ProjMassT; <<
271 <<
272 TPlusNew=ChooseP(TPlusMin, TPlusMax); <<
273 Qplus=-(TPlusNew-Ptarget.plus()); <<
274 <<
275 Qmomentum.setPz( (Qplus-Qminus)/2 ); <<
276 Qmomentum.setE( (Qplus+Qminus)/2 ); <<
277 <<
278 } while ( (Pprojectile+Qmomentum).mag2() < <<
279 (Ptarget -Qmomentum).mag2() < MinTrD <<
280 <<
281 Pprojectile += Qmomentum; <<
282 Ptarget -= Qmomentum; <<
283 <<
284 // Transform back and update SplitableHadron <<
285 Pprojectile.transform(toLab); <<
286 Ptarget.transform(toLab); <<
287 <<
288 #ifdef debugDoubleDiffraction <<
289 G4cout << "Pprojectile after boost to Lab: <<
290 G4cout << "Ptarget after boost to Lab: <<
291 #endif <<
292 71
293 target->Set4Momentum(Ptarget); << 72 G4LorentzVector Pprojectile=projectile->Get4Momentum();
294 projectile->Set4Momentum(Pprojectile); << 73
>> 74 // -------------------- Projectile parameters -----------------------------------
>> 75 G4bool PutOnMassShell=0;
>> 76
>> 77 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
>> 78 G4double M0projectile = Pprojectile.mag(); // Without de-excitation
>> 79
>> 80 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
>> 81 {
>> 82 PutOnMassShell=1;
>> 83 M0projectile=projectile->GetDefinition()->GetPDGMass();
>> 84 }
>> 85
>> 86 G4double Mprojectile2 = M0projectile * M0projectile;
>> 87
>> 88 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();
>> 89 G4int absPDGcode=std::abs(PDGcode);
>> 90 G4double ProjectileDiffCut;
>> 91 G4double AveragePt2;
>> 92
>> 93 if( absPDGcode > 1000 ) //------Projectile is baryon --------
>> 94 {
>> 95 ProjectileDiffCut = 1.1; // GeV
>> 96 AveragePt2 = 0.3; // GeV^2
>> 97 }
>> 98 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion -----------
>> 99 {
>> 100 ProjectileDiffCut = 1.0; // GeV
>> 101 AveragePt2 = 0.3; // GeV^2
>> 102 }
>> 103 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
>> 104 {
>> 105 ProjectileDiffCut = 1.1; // GeV
>> 106 AveragePt2 = 0.3; // GeV^2
>> 107 }
>> 108 else //------Projectile is undefined, Nucleon assumed
>> 109 {
>> 110 ProjectileDiffCut = 1.1; // GeV
>> 111 AveragePt2 = 0.3; // GeV^2
>> 112 };
>> 113
>> 114 ProjectileDiffCut = ProjectileDiffCut * GeV;
>> 115 AveragePt2 = AveragePt2 * GeV*GeV;
>> 116
>> 117 // -------------------- Target parameters ----------------------------------------------
>> 118 G4LorentzVector Ptarget=target->Get4Momentum();
>> 119
>> 120 G4double M0target = Ptarget.mag();
>> 121
>> 122 if(M0target < target->GetDefinition()->GetPDGMass())
>> 123 {
>> 124 PutOnMassShell=1;
>> 125 M0target=target->GetDefinition()->GetPDGMass();
>> 126 }
>> 127
>> 128 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
>> 129
>> 130 G4double NuclearNucleonDiffCut = 1.1*GeV;
>> 131
>> 132 G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
>> 133 G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
>> 134
>> 135 // Transform momenta to cms and then rotate parallel to z axis;
>> 136
>> 137 G4LorentzVector Psum;
>> 138 Psum=Pprojectile+Ptarget;
>> 139
>> 140 G4LorentzRotation toCms(-1*Psum.boostVector());
>> 141
>> 142 G4LorentzVector Ptmp=toCms*Pprojectile;
>> 143
>> 144 if ( Ptmp.pz() <= 0. )
>> 145 {
>> 146 // "String" moving backwards in CMS, abort collision !!
>> 147 //G4cout << " abort Collision!! " << G4endl;
>> 148 return false;
>> 149 }
>> 150
>> 151 toCms.rotateZ(-1*Ptmp.phi());
>> 152 toCms.rotateY(-1*Ptmp.theta());
>> 153
>> 154 G4LorentzRotation toLab(toCms.inverse());
>> 155
>> 156 Pprojectile.transform(toCms);
>> 157 Ptarget.transform(toCms);
>> 158
>> 159 G4double Pt2;
>> 160 G4double ProjMassT2, ProjMassT;
>> 161 G4double TargMassT2, TargMassT;
>> 162 G4double PZcms2, PZcms;
>> 163 G4double PMinusNew, TPlusNew;
>> 164
>> 165 G4double S=Psum.mag2();
>> 166 G4double SqrtS=std::sqrt(S);
>> 167
>> 168 if(SqrtS < 2200*MeV) {return false;} // The model cannot work for pp-interactions
>> 169 // at Plab < 1.3 GeV/c. Uzhi
>> 170
>> 171 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
>> 172 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
>> 173 if(PZcms2 < 0)
>> 174 {return false;} // It can be in an interaction with off-shell nuclear nucleon
>> 175
>> 176 PZcms = std::sqrt(PZcms2);
>> 177
>> 178 if(PutOnMassShell)
>> 179 {
>> 180 if(Pprojectile.z() > 0.)
>> 181 {
>> 182 Pprojectile.setPz( PZcms);
>> 183 Ptarget.setPz( -PZcms);
>> 184 }
>> 185 else
>> 186 {
>> 187 Pprojectile.setPz(-PZcms);
>> 188 Ptarget.setPz( PZcms);
>> 189 };
>> 190
>> 191 Pprojectile.setE(std::sqrt(Mprojectile2+
>> 192 Pprojectile.x()*Pprojectile.x()+
>> 193 Pprojectile.y()*Pprojectile.y()+
>> 194 PZcms2));
>> 195 Ptarget.setE(std::sqrt( Mtarget2 +
>> 196 Ptarget.x()*Ptarget.x()+
>> 197 Ptarget.y()*Ptarget.y()+
>> 198 PZcms2));
>> 199 }
>> 200
>> 201 G4double maxPtSquare = PZcms2;
>> 202
>> 203 //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
>> 204 //G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
>> 205 // G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
>> 206
>> 207 // G4cout << " Projectile Xplus / Xminus : " <<
>> 208 // Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
>> 209 // G4cout << " Target Xplus / Xminus : " <<
>> 210 // Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
>> 211
>> 212 G4LorentzVector Qmomentum;
>> 213 G4double Qminus, Qplus;
>> 214
>> 215 // /* Vova
>> 216 G4int whilecount=0;
>> 217 do {
>> 218 // Generate pt
>> 219
>> 220 if (whilecount++ >= 500 && (whilecount%100)==0)
>> 221 // G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping"
>> 222 // << ", loop count/ maxPtSquare : "
>> 223 // << whilecount << " / " << maxPtSquare << G4endl;
>> 224 if (whilecount > 1000 )
>> 225 {
>> 226 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
>> 227 return false; // Ignore this interaction
>> 228 }
>> 229
>> 230 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
>> 231
>> 232 //G4cout << "generated Pt " << Qmomentum << G4endl;
>> 233 //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
>> 234 //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
>> 235
>> 236 // Momentum transfer
>> 237 /* // Uzhi
>> 238 G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
>> 239 G4double Xmax=1.;
>> 240 G4double Xplus =ChooseX(Xmin,Xmax);
>> 241 G4double Xminus=ChooseX(Xmin,Xmax);
>> 242
>> 243 // G4cout << " X-plus " << Xplus << G4endl;
>> 244 // G4cout << " X-minus " << Xminus << G4endl;
>> 245
>> 246 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
>> 247 G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
>> 248 G4double Qminus= pt2 / Xplus /Pprojectile.plus();
>> 249 */ // Uzhi *
>> 250
>> 251 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
>> 252 ProjMassT2=Mprojectile2+Pt2;
>> 253 ProjMassT =std::sqrt(ProjMassT2);
>> 254
>> 255 TargMassT2=Mtarget2+Pt2;
>> 256 TargMassT =std::sqrt(TargMassT2);
>> 257
>> 258 PZcms2=(S*S+ProjMassT2*ProjMassT2+
>> 259 TargMassT2*TargMassT2-
>> 260 2.*S*ProjMassT2-2.*S*TargMassT2-
>> 261 2.*ProjMassT2*TargMassT2)/4./S;
>> 262 if(PZcms2 < 0 ) {PZcms2=0;};
>> 263 PZcms =std::sqrt(PZcms2);
>> 264
>> 265 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
>> 266 G4double PMinusMax=SqrtS-TargMassT;
>> 267
>> 268 PMinusNew=ChooseP(PMinusMin,PMinusMax);
>> 269 Qminus=PMinusNew-Pprojectile.minus();
>> 270
>> 271 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
>> 272 G4double TPlusMax=SqrtS-ProjMassT;
>> 273
>> 274 TPlusNew=ChooseP(TPlusMin, TPlusMax);
>> 275 Qplus=-(TPlusNew-Ptarget.plus());
>> 276
>> 277 Qmomentum.setPz( (Qplus-Qminus)/2 );
>> 278 Qmomentum.setE( (Qplus+Qminus)/2 );
>> 279
>> 280 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
>> 281 // G4cout << "pt2" << pt2 << G4endl;
>> 282 // G4cout << "Qmomentum " << Qmomentum << G4endl;
>> 283 // G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
>> 284 // " / " << (Ptarget-Qmomentum).mag() << G4endl;
>> 285 /* // Uzhi
>> 286 } while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 ||
>> 287 (Ptarget-Qmomentum).mag2() <= Mtarget2 );
>> 288 */ // Uzhi *
>> 289
>> 290
>> 291 } while (( (Pprojectile+Qmomentum).mag2() < Mprojectile2 || // Uzhi No without excitation
>> 292 (Ptarget -Qmomentum).mag2() < Mtarget2 ) || // Uzhi
>> 293 ( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // Uzhi No double Diffraction
>> 294 (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );// Uzhi
>> 295
>> 296 if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi Projectile diffraction
>> 297 {
>> 298 G4double TMinusNew=SqrtS-PMinusNew;
>> 299 Qminus=Ptarget.minus()-TMinusNew;
>> 300 TPlusNew=TargMassT2/TMinusNew;
>> 301 Qplus=Ptarget.plus()-TPlusNew;
>> 302
>> 303 Qmomentum.setPz( (Qplus-Qminus)/2 );
>> 304 Qmomentum.setE( (Qplus+Qminus)/2 );
>> 305 }
>> 306 else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi Target diffraction
>> 307 {
>> 308 G4double PPlusNew=SqrtS-TPlusNew;
>> 309 Qplus=PPlusNew-Pprojectile.plus();
>> 310 PMinusNew=ProjMassT2/PPlusNew;
>> 311 Qminus=PMinusNew-Pprojectile.minus();
>> 312
>> 313 Qmomentum.setPz( (Qplus-Qminus)/2 );
>> 314 Qmomentum.setE( (Qplus+Qminus)/2 );
>> 315 };
>> 316
>> 317 Pprojectile += Qmomentum;
>> 318 Ptarget -= Qmomentum;
>> 319
>> 320 // Vova
>> 321
>> 322 /*
>> 323 Pprojectile.setPz(0.);
>> 324 Pprojectile.setE(SqrtS-M0target);
>> 325
>> 326 Ptarget.setPz(0.);
>> 327 Ptarget.setE(M0target);
>> 328 */
>> 329
>> 330 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
>> 331 //G4cout << "Ptarget with Q : " << Ptarget << G4endl;
>> 332
>> 333 // G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
>> 334 // G4cout << "Target back: " << toLab * Ptarget << G4endl;
>> 335
>> 336 // Transform back and update SplitableHadron Participant.
>> 337 Pprojectile.transform(toLab);
>> 338 Ptarget.transform(toLab);
>> 339
>> 340 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;
>> 341 //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;
>> 342
>> 343 //G4cout << "Target mass " << Ptarget.mag() << G4endl;
>> 344
>> 345 target->Set4Momentum(Ptarget);
295 346
296 return true; << 347 //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;
>> 348
>> 349 projectile->Set4Momentum(Pprojectile);
>> 350
>> 351 return true;
297 } 352 }
298 353
299 354
300 G4ExcitedString * G4QGSDiffractiveExcitation:: 355 G4ExcitedString * G4QGSDiffractiveExcitation::
301 String(G4VSplitableHadron * hadron, G4bool isP 356 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
302 { 357 {
303 hadron->SplitUp(); << 358 hadron->SplitUp();
304 G4Parton *start= hadron->GetNextParton(); << 359 G4Parton *start= hadron->GetNextParton();
305 if ( start==NULL) << 360 if ( start==NULL)
306 { G4cout << " G4QGSDiffractiveExcitation::St << 361 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
307 return NULL; << 362 return NULL;
308 } << 363 }
309 G4Parton *end = hadron->GetNextParton(); << 364 G4Parton *end = hadron->GetNextParton();
310 if ( end==NULL) << 365 if ( end==NULL)
311 { G4cout << " G4QGSDiffractiveExcitation::St << 366 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
312 return NULL; << 367 return NULL;
313 } << 368 }
314 << 369
315 G4ExcitedString * string; << 370 G4ExcitedString * string;
316 if ( isProjectile ) << 371 if ( isProjectile )
317 { << 372 {
318 string= new G4ExcitedString(end,start, +1) << 373 string= new G4ExcitedString(end,start, +1);
319 } else { << 374 } else {
320 string= new G4ExcitedString(start,end, -1) << 375 string= new G4ExcitedString(start,end, -1);
321 } << 376 }
322 << 377
323 string->SetPosition(hadron->GetPosition()); << 378 string->SetPosition(hadron->GetPosition());
324 << 379
325 // momenta of string ends << 380 // momenta of string ends
326 << 381 G4double ptSquared= hadron->Get4Momentum().perp2();
327 G4double maxAvailMomentumSquared=sqr(hadron- << 382 G4double transverseMassSquared= hadron->Get4Momentum().plus()
328 << 383 * hadron->Get4Momentum().minus();
329 G4double widthOfPtSquare = 0.5*sqr(GeV); << 384
330 G4ThreeVector pt=GaussianPt(widthOfPtSquare, << 385
331 << 386 G4double maxAvailMomentumSquared=
332 G4LorentzVector Pstart(G4LorentzVector(pt,0. << 387 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
333 G4LorentzVector Pend; << 388
334 Pend.setPx(hadron->Get4Momentum().px() - pt. << 389 G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ??????????????????
335 Pend.setPy(hadron->Get4Momentum().py() - pt. << 390 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
336 << 391
337 G4double tm1=hadron->Get4Momentum().minus() << 392 G4LorentzVector Pstart(G4LorentzVector(pt,0.));
338 ( Pend.perp2()-Pstart.perp2() ) / had << 393 G4LorentzVector Pend;
339 << 394 Pend.setPx(hadron->Get4Momentum().px() - pt.x());
340 G4double tm2= std::sqrt( std::max(0., sqr(tm << 395 Pend.setPy(hadron->Get4Momentum().py() - pt.y());
341 4. * Pend.perp2() * hadron->Get4Momen << 396
342 / hadron->Get4Momentum().plus() )); << 397 G4double tm1=hadron->Get4Momentum().minus() +
343 << 398 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
344 G4int Sign= isProjectile ? -1 : 1; << 399
345 << 400 G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
346 G4double endMinus = 0.5 * (tm1 + Sign*tm2); << 401 4. * Pend.perp2() * hadron->Get4Momentum().minus()
347 G4double startMinus= hadron->Get4Momentum(). << 402 / hadron->Get4Momentum().plus() ));
348 << 403
349 G4double startPlus= Pstart.perp2() / startM << 404 G4int Sign= isProjectile ? -1 : 1;
350 G4double endPlus = hadron->Get4Momentum().p << 405
351 << 406 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
352 Pstart.setPz(0.5*(startPlus - startMinus)); << 407 G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
353 Pstart.setE(0.5*(startPlus + startMinus)); << 408
354 << 409 G4double startPlus= Pstart.perp2() / startMinus;
355 Pend.setPz(0.5*(endPlus - endMinus)); << 410 G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
356 Pend.setE(0.5*(endPlus + endMinus)); << 411
357 << 412 Pstart.setPz(0.5*(startPlus - startMinus));
358 start->Set4Momentum(Pstart); << 413 Pstart.setE(0.5*(startPlus + startMinus));
359 end->Set4Momentum(Pend); << 414
360 << 415 Pend.setPz(0.5*(endPlus - endMinus));
361 #ifdef debugQGSdiffExictation << 416 Pend.setE(0.5*(endPlus + endMinus));
362 G4cout << " generated string flavors << 417
363 G4cout << " generated string momenta: qu << 418 start->Set4Momentum(Pstart);
364 G4cout << " generated string momenta: Diqu << 419 end->Set4Momentum(Pend);
365 G4cout << " sum of ends << 420
366 G4cout << " Original << 421 #ifdef G4_FTFDEBUG
367 #endif << 422 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
>> 423 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
>> 424 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
>> 425 G4cout << " sum of ends " << Pstart+Pend << G4endl;
>> 426 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
>> 427 #endif
368 428
369 return string; << 429 return string;
370 } 430 }
371 431
372 432
373 // --------- private methods ----------------- 433 // --------- private methods ----------------------
374 434
375 G4double G4QGSDiffractiveExcitation::ChooseP(G << 435 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
376 { 436 {
377 // choose an x between Xmin and Xmax with P( << 437 // choose an x between Xmin and Xmax with P(x) ~ 1/x
378 // to be improved... <<
379 438
380 G4double range=Pmax-Pmin; << 439 // to be improved...
381 440
382 if ( Pmin <= 0. || range <=0. ) << 441 G4double range=Pmax-Pmin; // Uzhi
383 { << 442
384 G4cout << " Pmin, range : " << Pmin << " , << 443 if ( Pmin <= 0. || range <=0. )
385 throw G4HadronicException(__FILE__, __LINE << 444 {
386 } << 445 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
387 << 446 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
388 G4double P; << 447 }
389 P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmi <<
390 //debug-hpw cout << "DiffractiveX "<<x<<G4en <<
391 return P; <<
392 } <<
393 448
>> 449 G4double P;
>> 450 /* // Uzhi
>> 451 do {
>> 452 x=Xmin + G4UniformRand() * range;
>> 453 } while ( Xmin/x < G4UniformRand() );
>> 454 */ // Uzhi
394 455
395 G4ThreeVector G4QGSDiffractiveExcitation::Gaus << 456 P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); // Uzhi
>> 457
>> 458 //debug-hpw cout << "DiffractiveX "<<x<<G4endl;
>> 459 return P;
>> 460 }
>> 461
>> 462 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi
396 { // @@ this method is used in FTF 463 { // @@ this method is used in FTFModel as well. Should go somewhere common!
397 G4double Pt2; <<
398 464
399 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand << 465 G4double Pt2;
>> 466 /* // Uzhi
>> 467 do {
>> 468 pt2=widthSquare * std::log( G4UniformRand() );
>> 469 } while ( pt2 > maxPtSquare);
>> 470 */ // Uzhi
>> 471
>> 472 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
400 473
401 G4double Pt=std::sqrt(Pt2); << 474 G4double Pt=std::sqrt(Pt2);
402 475
403 G4double phi=G4UniformRand() * twopi; << 476 G4double phi=G4UniformRand() * twopi;
404 477
405 return G4ThreeVector (Pt*std::cos(phi), Pt*s << 478 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
406 } 479 }
407 480