source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4BigBanger.cc @ 962

Last change on this file since 962 was 962, checked in by garnier, 15 years ago

update processes
File size: 10.8 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#include "G4BigBanger.hh"
27#include "G4InuclNuclei.hh"
28#include "G4ParticleLargerEkin.hh"
29#include "G4LorentzConvertor.hh"
30#include <algorithm>
31
32typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
33
34G4BigBanger::G4BigBanger()
35  : verboseLevel(1) {
36  if (verboseLevel > 3) {
37    G4cout << " >>> G4BigBanger::G4BigBanger" << G4endl;
38  }
39}
40
41G4CollisionOutput G4BigBanger::collide(G4InuclParticle* /*bullet*/,
42                                       G4InuclParticle* target) {
43
44  if (verboseLevel > 3) {
45    G4cout << " >>> G4BigBanger::collide" << G4endl;
46  }
47
48  // primitive explosion model A -> nucleons to prevent too exotic evaporation
49
50  const G4double small_ekin = 1.0e-6;
51
52  G4CollisionOutput output;
53  G4CascadeMomentum totscm;
54  G4CascadeMomentum totlab;
55
56  if(G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
57 
58    G4double A = nuclei_target->getA();
59    G4double Z = nuclei_target->getZ();
60    const G4CascadeMomentum& PEX = nuclei_target->getMomentum();
61    G4double EEXS = nuclei_target->getExitationEnergy();
62    G4InuclElementaryParticle dummy(small_ekin, 1);
63    G4LorentzConvertor toTheNucleiSystemRestFrame;
64
65    toTheNucleiSystemRestFrame.setBullet(dummy.getMomentum(), dummy.getMass());
66    toTheNucleiSystemRestFrame.setTarget(PEX, nuclei_target->getMass());
67    toTheNucleiSystemRestFrame.toTheTargetRestFrame();
68
69    G4double etot = 0.001 * (EEXS - bindingEnergy(A, Z));
70
71    if (verboseLevel > 2) {
72      G4cout << " BigBanger: target " << G4endl;
73      nuclei_target->printParticle(); 
74      G4cout << " BigBanger: a " << A << " z " << Z << " eexs " << EEXS << " etot " <<
75        etot << " nm " << nuclei_target->getMass() << G4endl;
76    }
77 
78    std::vector<G4InuclElementaryParticle> particles =     
79      generateBangInSCM(etot, A, Z, dummy.getParticleMass(1), dummy.getParticleMass(2));
80
81    if (verboseLevel > 2) {
82      G4cout << " particles " << particles.size() << G4endl;
83      for(G4int i = 0; i < G4int(particles.size()); i++) 
84        particles[i].printParticle();
85    }
86    if(!particles.empty()) { // convert back to Lab
87      //      if (verboseLevel > 2) {
88      // not used    G4CascadeMomentum totscm;
89      // not used    G4CascadeMomentum totlab;
90      //      }
91      particleIterator ipart;
92
93      for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
94        if (verboseLevel > 2) {
95          const G4CascadeMomentum& mom_scm = ipart->getMomentum();
96
97          for(G4int i = 0; i < 4; i++) totscm[i] += mom_scm[i];
98        }
99        G4CascadeMomentum mom = 
100          toTheNucleiSystemRestFrame.backToTheLab(ipart->getMomentum());
101        ipart->setMomentum(mom); 
102
103        if (verboseLevel > 2) {
104          mom = ipart->getMomentum();
105          for(G4int i = 0; i < 4; i++) totlab[i] += mom[i];
106        }
107      };
108      std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
109      if (verboseLevel > 2) {
110        G4cout << " In SCM: total outgoing momentum " << G4endl
111               << " E " << totscm[0] << " px " << totscm[1]
112               << " py " << totscm[2] << " pz " << totscm[3] << G4endl; 
113        G4cout << " In Lab: mom cons " << G4endl
114               << " E " << PEX[0] + 0.001 * EEXS - totlab[0] 
115               << " px " << PEX[1] - totlab[1]
116               << " py " << PEX[2] - totlab[2] 
117               << " pz " << PEX[3] - totlab[3] << G4endl; 
118      }
119    }; 
120    output.addOutgoingParticles(particles);
121  }
122  else {
123    G4cout << " BigBanger -> try to bang not nuclei " << G4endl;
124  }; 
125
126  return output;
127}                   
128
129std::vector<G4InuclElementaryParticle>             
130G4BigBanger::generateBangInSCM(G4double etot, 
131                               G4double a, 
132                               G4double z, 
133                               G4double mp,
134                               G4double mn) const {
135
136  if (verboseLevel > 3) {
137    G4cout << " >>> G4BigBanger::generateBangInSCM" << G4endl;
138  }
139
140  const G4double ang_cut = 0.9999;
141  const G4int itry_max = 1000;
142 
143  G4int ia = int(a + 0.1);
144  G4int iz = int(z + 0.1);
145
146  if (verboseLevel > 2) {
147    G4cout << " ia " << ia << " iz " << iz << G4endl;
148  }
149  std::vector<G4InuclElementaryParticle> particles;
150 
151  if(ia == 1) {
152    // abnormal situation
153    G4double m = iz > 0 ? mp : mn;
154    G4double pmod = std::sqrt((etot + 2.0 * m) * etot);
155    G4CascadeMomentum mom;
156    std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
157    G4double FI = randomPHI();
158    G4double Pt = pmod * COS_SIN.second;
159
160    mom[1] = Pt * std::cos(FI);
161    mom[2] = Pt * std::sin(FI);
162    mom[3] = Pt * COS_SIN.first;   
163
164    G4int knd = iz > 0 ? 1 : 2;
165
166    //    particles.push_back(G4InuclElementaryParticle(mom, knd));
167    particles.push_back(G4InuclElementaryParticle(mom, knd, 8)); // modelId included
168
169    return particles;
170  }; 
171     
172  std::vector<G4double> pmod = generateMomentumModules(etot, a, z, mp, mn);
173  G4bool bad = true;
174  G4int itry = 0;
175
176  while(bad && itry < itry_max) {
177    itry++;
178    std::vector<G4CascadeMomentum> scm_momentums;
179    G4CascadeMomentum tot_mom;
180
181    if(ia == 2) {
182      G4CascadeMomentum mom;
183      std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
184      double FI = randomPHI();
185      double Pt = pmod[0] * COS_SIN.second;
186
187      mom[1] = Pt * std::cos(FI);
188      mom[2] = Pt * std::sin(FI);
189      mom[3] = Pt * COS_SIN.first;   
190
191      for(G4int j = 1; j < 4; j++) tot_mom[j] += mom[j];                 
192
193      scm_momentums.push_back(mom);
194
195      G4CascadeMomentum mom1;
196
197      for(G4int i = 1; i < 4; i++) mom1[i] = - mom[i];
198
199      scm_momentums.push_back(mom1); 
200      bad = false;
201    }
202    else {
203      for(G4int i = 0; i < ia - 2; i++) {
204        G4CascadeMomentum mom;
205        std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
206        G4double FI = randomPHI();
207        G4double Pt = pmod[i] * COS_SIN.second;
208
209        mom[1] = Pt * std::cos(FI);
210        mom[2] = Pt * std::sin(FI);
211        mom[3] = Pt * COS_SIN.first;   
212
213        for(G4int j = 1; j < 4; j++) tot_mom[j] += mom[j];               
214
215        scm_momentums.push_back(mom);
216      };
217
218      //                handle last two
219      G4double tot_mod = std::sqrt(tot_mom[1] * tot_mom[1] + 
220                              tot_mom[2] * tot_mom[2] + 
221                              tot_mom[3] * tot_mom[3]); 
222      G4double ct = -0.5 * (tot_mod * tot_mod + pmod[ia - 2] * pmod[ia - 2] -
223                            pmod[ia - 1] * pmod[ia - 1]) / tot_mod / pmod[ia - 2];
224
225      if (verboseLevel > 2) {
226        G4cout << " ct last " << ct << G4endl;
227      }
228 
229      if(std::fabs(ct) < ang_cut) {
230        G4CascadeMomentum mom2 = generateWithFixedTheta(ct, pmod[ia - 2]);
231        //       rotate to the normal system
232        G4CascadeMomentum apr = tot_mom;
233        G4int i;
234        for(i = 1; i < 4; i++) apr[i] /= tot_mod;
235        G4double a_tr = std::sqrt(apr[1] * apr[1] + apr[2] * apr[2]);
236        G4CascadeMomentum mom;
237        mom[1] = mom2[3] * apr[1] + ( mom2[1] * apr[2] + mom2[2] * apr[3] * apr[1]) / a_tr; // ::: replace with clhep tools?
238        mom[2] = mom2[3] * apr[2] + (-mom2[1] * apr[1] + mom2[2] * apr[3] * apr[2]) / a_tr;     
239        mom[3] = mom2[3] * apr[3] - mom2[2] * a_tr;     
240        scm_momentums.push_back(mom);
241        //               and the last one
242        G4CascadeMomentum mom1;
243        for(i = 1; i < 4; i++) mom1[i] = - mom[i] - tot_mom[i];
244        scm_momentums.push_back(mom1); 
245        bad = false;
246      };
247    };   
248    if(!bad) {
249      for(G4int i = 0; i < ia; i++) {
250        G4int knd = i < iz ? 1 : 2;
251
252        particles.push_back(G4InuclElementaryParticle(scm_momentums[i], knd));
253      };
254    };
255  }; 
256  if (verboseLevel > 2) {
257    if(itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl;
258  }
259
260  return particles;
261 
262}
263           
264std::vector<G4double> G4BigBanger::generateMomentumModules(G4double etot, 
265                                                             G4double a, 
266                                                             G4double z, 
267                                                             G4double mp, 
268                                                             G4double mn) const {
269
270
271  if (verboseLevel > 3) {
272    G4cout << " >>> G4BigBanger::generateMomentumModules" << G4endl;
273  }
274
275  G4int ia = int(a + 0.1);
276  G4int iz = int(z + 0.1);
277  std::vector<G4double> pmod;
278  G4double xtot = 0.0;
279  G4double promax = maxProbability(a);
280 
281  G4int i;
282  for(i = 0; i < ia; i++) { 
283    G4double x = generateX(ia, a, promax);
284
285    if (verboseLevel > 2) {
286      G4cout << " i " << i << " x " << x << G4endl;
287    }
288    pmod.push_back(x);
289    xtot += x;
290  };
291  for(i = 0; i < ia; i++) {
292    G4double m = i < iz ? mp : mn;
293
294    pmod[i] = pmod[i] * etot / xtot;
295    pmod[i] = std::sqrt(pmod[i] * (pmod[i] + 2.0 * m));
296
297    if (verboseLevel > 2) {
298      G4cout << " i " << i << " pmod " << pmod[i] << G4endl;
299    }
300  };
301
302  return pmod; 
303}
304
305G4double G4BigBanger::xProbability(G4double x, 
306                                   G4int ia) const {
307
308
309  if (verboseLevel > 3) {
310    G4cout << " >>> G4BigBanger::xProbability" << G4endl;
311  }
312
313  G4int ihalf = ia / 2;
314  G4double ekpr = 0.0;
315
316  if(x < 1.0 || x > 0.0) {
317    ekpr = x * x;
318
319    if(2 * ihalf == ia) { // even A
320      ekpr *= std::sqrt(1.0 - x) * std::pow((1.0 - x), G4int(G4double(3 * ia - 6) / 2.0)); 
321    }
322    else {
323      ekpr *= std::pow((1.0 - x), G4int(G4double(3 * ia - 5) / 2.0));
324    };
325  }; 
326 
327  return ekpr;
328}
329
330G4double G4BigBanger::maxProbability(G4double a) const {
331
332  if (verboseLevel > 3) {
333    G4cout << " >>> G4BigBanger::maxProbability" << G4endl;
334  }
335
336  return xProbability(1.0 / (a - 1.0) / 1.5, G4int(a + 0.1));
337}
338
339G4double G4BigBanger::generateX(G4int ia, 
340                                G4double a, 
341                                G4double promax) const {
342
343  if (verboseLevel > 3) {
344    G4cout << " >>> G4BigBanger::generateX" << G4endl;
345  }
346
347  const G4int itry_max = 1000;
348  G4int itry = 0;
349  G4double x;
350 
351  while(itry < itry_max) {
352    itry++;
353    x = inuclRndm();
354
355    if(xProbability(x, ia) >= promax * inuclRndm()) return x;
356  };
357  if (verboseLevel > 2) {
358    G4cout << " BigBanger -> can not generate x " << G4endl;
359  }
360
361  return maxProbability(a);
362}
Note: See TracBrowser for help on using the repository browser.