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

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 21.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4DiffractiveExcitation.cc,v 1.1 2007/05/25 06:56:53 vuzhinsk Exp $
28// ------------------------------------------------------------
29// GEANT 4 class implemetation file
30//
31// ---------------- G4DiffractiveExcitation --------------
32// by Gunter Folger, October 1998.
33// diffractive Excitation used by strings models
34// Take a projectile and a target
35// excite the projectile and target
36// Essential changed by V. Uzhinsky in November - December 2006
37// in order to put it in a correspondence with original FRITIOF
38// model. Variant of FRITIOF with nucleon de-excitation is implemented.
39// Other changes by V.Uzhinsky in May 2007 were introduced to fit
40// meson-nucleon interactions
41// ---------------------------------------------------------------------
42
43
44#include "globals.hh"
45#include "Randomize.hh"
46
47#include "G4DiffractiveExcitation.hh"
48#include "G4LorentzRotation.hh"
49#include "G4ThreeVector.hh"
50#include "G4ParticleDefinition.hh"
51#include "G4VSplitableHadron.hh"
52#include "G4ExcitedString.hh"
53//#include "G4ios.hh"
54
55G4DiffractiveExcitation::G4DiffractiveExcitation() // Uzhi
56{
57}
58
59G4bool G4DiffractiveExcitation::
60 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
61{
62
63 G4LorentzVector Pprojectile=projectile->Get4Momentum();
64
65// -------------------- Projectile parameters -----------------------------------
66 G4bool PutOnMassShell=0;
67
68// With de-excitation
69// G4double M0projectile=projectile->GetDefinition()->GetPDGMass();
70// Without de-excitation
71 G4double M0projectile = Pprojectile.mag();
72
73 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
74 {
75 PutOnMassShell=1;
76 M0projectile=projectile->GetDefinition()->GetPDGMass();
77 }
78
79 G4double Mprojectile2 = M0projectile * M0projectile;
80
81 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();
82 G4int absPDGcode=std::abs(PDGcode);
83 G4double ProjectileDiffCut;
84 G4double ProjectileBackgroundWeight; // Uzhi May 2007
85
86 G4double TargetBackgroundWeight; // Uzhi May 2007
87
88 G4double AveragePt2;
89
90 if( absPDGcode > 1000 ) //------Projectile is baryon --------
91 {
92 ProjectileDiffCut = 1.1; // GeV
93 ProjectileBackgroundWeight=0.; // Uzhi May 2007
94 AveragePt2 = 0.3; // GeV^2
95 }
96 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion -----------
97 {
98 ProjectileDiffCut = 0.6; // GeV Uzhi May 2007 1.0->0.6
99 ProjectileBackgroundWeight=0.9; // Uzhi May 2007
100 AveragePt2 = 0.3; // GeV^2
101 }
102 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
103 {
104 ProjectileDiffCut = 0.7; // GeV 1.1
105 ProjectileBackgroundWeight=0.5; // Uzhi May 2007 ???
106 AveragePt2 = 0.3; // GeV^2
107 }
108 else //------Projectile is undefined,
109 //------Nucleon assumed
110 {
111 ProjectileDiffCut = 1.1; // GeV
112 ProjectileBackgroundWeight=0.; // Uzhi May 2007
113 AveragePt2 = 0.3; // GeV^2
114 };
115
116 ProjectileDiffCut = ProjectileDiffCut * GeV;
117 AveragePt2 = AveragePt2 * GeV*GeV;
118
119// -------------------- Target parameters ----------------------------------------------
120 G4LorentzVector Ptarget=target->Get4Momentum();
121
122 G4double M0target = Ptarget.mag();
123
124 if(M0target < target->GetDefinition()->GetPDGMass())
125 {
126 PutOnMassShell=1;
127 M0target=target->GetDefinition()->GetPDGMass();
128 }
129
130 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2();
131 // for AA-inter.
132
133 G4double NuclearNucleonDiffCut = 1.1*GeV;
134 TargetBackgroundWeight=0.; // Uzhi May 2007
135
136 G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
137 G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
138
139// Transform momenta to cms and then rotate parallel to z axis;
140
141 G4LorentzVector Psum;
142 Psum=Pprojectile+Ptarget;
143
144 G4LorentzRotation toCms(-1*Psum.boostVector());
145
146 G4LorentzVector Ptmp=toCms*Pprojectile;
147
148 if ( Ptmp.pz() <= 0. )
149 {
150 // "String" moving backwards in CMS, abort collision !!
151 //G4cout << " abort Collision!! " << G4endl;
152 return false;
153 }
154
155 toCms.rotateZ(-1*Ptmp.phi());
156 toCms.rotateY(-1*Ptmp.theta());
157
158 G4LorentzRotation toLab(toCms.inverse());
159
160 Pprojectile.transform(toCms);
161 Ptarget.transform(toCms);
162
163 G4double Pt2;
164 G4double ProjMassT2, ProjMassT;
165 G4double TargMassT2, TargMassT;
166 G4double PZcms2, PZcms;
167 G4double PMinusNew, TPlusNew;
168
169 G4double S=Psum.mag2();
170 G4double SqrtS=std::sqrt(S);
171
172 if(absPDGcode > 1000 && SqrtS < 2200*MeV)
173 {return false;} // The model cannot work for
174 // p+p-interactions
175 // at Plab < 1.3 GeV/c.
176
177 if(( absPDGcode == 211 || PDGcode == 111) && SqrtS < 1600*MeV)
178 {return false;} // The model cannot work for
179 // Pi+p-interactions
180 // at Plab < 1. GeV/c.
181
182 if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV)
183 {return false;} // The model cannot work for
184 // K+p-interactions
185 // at Plab < ??? GeV/c. ???
186
187 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
188 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
189 if(PZcms2 < 0)
190 {return false;} // It can be in an interaction with off-shell nuclear nucleon
191
192 PZcms = std::sqrt(PZcms2);
193
194 if(PutOnMassShell)
195 {
196 if(Pprojectile.z() > 0.)
197 {
198 Pprojectile.setPz( PZcms);
199 Ptarget.setPz( -PZcms);
200 }
201 else
202 {
203 Pprojectile.setPz(-PZcms);
204 Ptarget.setPz( PZcms);
205 };
206
207 Pprojectile.setE(std::sqrt(Mprojectile2+
208 Pprojectile.x()*Pprojectile.x()+
209 Pprojectile.y()*Pprojectile.y()+
210 PZcms2));
211 Ptarget.setE(std::sqrt( Mtarget2 +
212 Ptarget.x()*Ptarget.x()+
213 Ptarget.y()*Ptarget.y()+
214 PZcms2));
215 }
216
217 G4double maxPtSquare = PZcms2;
218
219//G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
220//G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
221// G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
222
223// G4cout << " Projectile Xplus / Xminus : " <<
224// Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
225// G4cout << " Target Xplus / Xminus : " <<
226// Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
227
228 G4LorentzVector Qmomentum;
229 G4double Qminus, Qplus;
230
231 G4int whilecount=0;
232 do {
233// Generate pt
234
235 if (whilecount++ >= 500 && (whilecount%100)==0)
236// G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
237// << ", loop count/ maxPtSquare : "
238// << whilecount << " / " << maxPtSquare << G4endl;
239 if (whilecount > 1000 )
240 {
241 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
242 return false; // Ignore this interaction
243 }
244
245 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
246
247//G4cout << "generated Pt " << Qmomentum << G4endl;
248//G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
249//G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
250
251// Momentum transfer
252/* // Uzhi
253 G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
254 G4double Xmax=1.;
255 G4double Xplus =ChooseX(Xmin,Xmax);
256 G4double Xminus=ChooseX(Xmin,Xmax);
257
258// G4cout << " X-plus " << Xplus << G4endl;
259// G4cout << " X-minus " << Xminus << G4endl;
260
261 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
262 G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
263 G4double Qminus= pt2 / Xplus /Pprojectile.plus();
264*/ // Uzhi *
265
266 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
267 ProjMassT2=Mprojectile2+Pt2;
268 ProjMassT =std::sqrt(ProjMassT2);
269
270 TargMassT2=Mtarget2+Pt2;
271 TargMassT =std::sqrt(TargMassT2);
272
273 PZcms2=(S*S+ProjMassT2*ProjMassT2+
274 TargMassT2*TargMassT2-
275 2.*S*ProjMassT2-2.*S*TargMassT2-
276 2.*ProjMassT2*TargMassT2)/4./S;
277 if(PZcms2 < 0 ) {PZcms2=0;};
278 PZcms =std::sqrt(PZcms2);
279
280 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
281 G4double PMinusMax=SqrtS-TargMassT;
282
283 if(G4UniformRand() < ProjectileBackgroundWeight) // Uzhi May 2007
284 {
285 PMinusNew=PMinusMax-(PMinusMax-PMinusMin)*G4UniformRand(); // Uzhi May 2007
286 } // Uzhi May 2007
287 else // Uzhi May 2007
288 { // Uzhi May 2007
289 PMinusNew=ChooseP(PMinusMin, PMinusMax); // Uzhi May 2007
290 }; // Uzhi May 2007
291// PMinusNew=ChooseP(PMinusMin, PMinusMax); // Uzhi May 2007
292 Qminus=PMinusNew-Pprojectile.minus();
293
294 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
295 G4double TPlusMax=SqrtS-PMinusNew; // Uzhi May 2007
296// G4double TPlusMax=SqrtS-ProjMassT; // Uzhi May 2007
297
298 if(G4UniformRand() < TargetBackgroundWeight) // Uzhi May 2007
299 {
300 TPlusNew=TPlusMax-(TPlusMax-TPlusMin)*G4UniformRand(); // Uzhi May 2007
301 } // Uzhi May 2007
302 else // Uzhi May 2007
303 { // Uzhi May 2007
304 TPlusNew=ChooseP(TPlusMin, TPlusMax); // Uzhi May 2007
305 }; // Uzhi May 2007
306
307// TPlusNew=ChooseP(TPlusMin, TPlusMax); // Uzhi May 2007
308 Qplus=-(TPlusNew-Ptarget.plus());
309
310 Qmomentum.setPz( (Qplus-Qminus)/2 );
311 Qmomentum.setE( (Qplus+Qminus)/2 );
312
313//G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
314// G4cout << "pt2" << pt2 << G4endl;
315// G4cout << "Qmomentum " << Qmomentum << G4endl;
316// G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
317// " / " << (Ptarget-Qmomentum).mag() << G4endl;
318
319 } while (
320( (Pprojectile+Qmomentum).mag2() < Mprojectile2 || // Uzhi No without excitation
321 (Ptarget -Qmomentum).mag2() < Mtarget2 ) || // Uzhi
322( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // Uzhi No double Diffraction
323 (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );// Uzhi
324
325 if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi
326//Projectile diffraction
327 {
328 G4double TMinusNew=SqrtS-PMinusNew;
329 Qminus=Ptarget.minus()-TMinusNew;
330 TPlusNew=TargMassT2/TMinusNew;
331 Qplus=Ptarget.plus()-TPlusNew;
332
333 Qmomentum.setPz( (Qplus-Qminus)/2 );
334 Qmomentum.setE( (Qplus+Qminus)/2 );
335 }
336 else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi
337//Target diffraction
338 {
339 G4double PPlusNew=SqrtS-TPlusNew;
340 Qplus=PPlusNew-Pprojectile.plus();
341 PMinusNew=ProjMassT2/PPlusNew;
342 Qminus=PMinusNew-Pprojectile.minus();
343
344 Qmomentum.setPz( (Qplus-Qminus)/2 );
345 Qmomentum.setE( (Qplus+Qminus)/2 );
346 };
347
348 Pprojectile += Qmomentum;
349 Ptarget -= Qmomentum;
350
351/*
352Pprojectile.setPz(0.);
353Pprojectile.setE(SqrtS-M0target);
354
355Ptarget.setPz(0.);
356Ptarget.setE(M0target);
357*/
358
359//G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
360//G4cout << "Ptarget with Q : " << Ptarget << G4endl;
361
362// G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
363// G4cout << "Target back: " << toLab * Ptarget << G4endl;
364
365// Transform back and update SplitableHadron Participant.
366 Pprojectile.transform(toLab);
367 Ptarget.transform(toLab);
368
369//G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;
370//G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;
371
372//G4cout << "Target mass " << Ptarget.mag() << G4endl;
373
374 target->Set4Momentum(Ptarget);
375
376//G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;
377
378 projectile->Set4Momentum(Pprojectile);
379
380 return true;
381}
382
383
384G4ExcitedString * G4DiffractiveExcitation::
385 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
386{
387 hadron->SplitUp();
388 G4Parton *start= hadron->GetNextParton();
389 if ( start==NULL)
390 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
391 return NULL;
392 }
393 G4Parton *end = hadron->GetNextParton();
394 if ( end==NULL)
395 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
396 return NULL;
397 }
398
399 G4ExcitedString * string;
400 if ( isProjectile )
401 {
402 string= new G4ExcitedString(end,start, +1);
403 } else {
404 string= new G4ExcitedString(start,end, -1);
405 }
406
407 string->SetPosition(hadron->GetPosition());
408
409// momenta of string ends
410 G4double ptSquared= hadron->Get4Momentum().perp2();
411 G4double transverseMassSquared= hadron->Get4Momentum().plus()
412 * hadron->Get4Momentum().minus();
413
414
415 G4double maxAvailMomentumSquared=
416 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
417
418 G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ??????????????????
419 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
420
421 G4LorentzVector Pstart(G4LorentzVector(pt,0.));
422 G4LorentzVector Pend;
423 Pend.setPx(hadron->Get4Momentum().px() - pt.x());
424 Pend.setPy(hadron->Get4Momentum().py() - pt.y());
425
426 G4double tm1=hadron->Get4Momentum().minus() +
427 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
428
429 G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
430 4. * Pend.perp2() * hadron->Get4Momentum().minus()
431 / hadron->Get4Momentum().plus() ));
432
433 G4int Sign= isProjectile ? -1 : 1;
434
435 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
436 G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
437
438 G4double startPlus= Pstart.perp2() / startMinus;
439 G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
440
441 Pstart.setPz(0.5*(startPlus - startMinus));
442 Pstart.setE(0.5*(startPlus + startMinus));
443
444 Pend.setPz(0.5*(endPlus - endMinus));
445 Pend.setE(0.5*(endPlus + endMinus));
446
447 start->Set4Momentum(Pstart);
448 end->Set4Momentum(Pend);
449
450#ifdef G4_FTFDEBUG
451 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
452 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
453 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
454 G4cout << " sum of ends " << Pstart+Pend << G4endl;
455 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
456#endif
457
458 return string;
459}
460
461
462// --------- private methods ----------------------
463
464G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
465{
466// choose an x between Xmin and Xmax with P(x) ~ 1/x
467
468// to be improved...
469
470 G4double range=Pmax-Pmin; // Uzhi
471
472 if ( Pmin <= 0. || range <=0. )
473 {
474 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
475 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
476 }
477
478 G4double P;
479/* // Uzhi
480 do {
481 x=Xmin + G4UniformRand() * range;
482 } while ( Xmin/x < G4UniformRand() );
483*/ // Uzhi
484
485 P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); // Uzhi
486
487//debug-hpw cout << "DiffractiveX "<<x<<G4endl;
488 return P;
489}
490
491G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi
492{ // @@ this method is used in FTFModel as well. Should go somewhere common!
493
494 G4double Pt2;
495/* // Uzhi
496 do {
497 pt2=widthSquare * std::log( G4UniformRand() );
498 } while ( pt2 > maxPtSquare);
499*/ // Uzhi
500
501 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
502
503 G4double Pt=std::sqrt(Pt2);
504
505 G4double phi=G4UniformRand() * twopi;
506
507 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
508}
509
510G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
511{
512 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
513}
514
515
516G4DiffractiveExcitation::~G4DiffractiveExcitation()
517{
518}
519
520
521const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
522{
523 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
524 return *this;
525}
526
527
528int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
529{
530 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
531 return false;
532}
533
534int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
535{
536 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
537 return true;
538}
Note: See TracBrowser for help on using the repository browser.