source: trunk/source/processes/hadronic/models/util/src/G4GeneralPhaseSpaceDecay.cc @ 1196

Last change on this file since 1196 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 20.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27//
28// $Id: G4GeneralSpaceDecay.cc,v 1.0 1998/05/21
29// ----------------------------------------------------------------
30//      GEANT 4 class header file
31//
32//      History: first implementation, A. Feliciello, 21st May 1998
33//
34//      Note: this class is a generalization of the
35//            G4PhaseSpaceDecayChannel one
36// ----------------------------------------------------------------
37
38#include "G4ParticleDefinition.hh"
39#include "G4DecayProducts.hh"
40#include "G4VDecayChannel.hh"
41#include "G4GeneralPhaseSpaceDecay.hh"
42#include "Randomize.hh"
43#include "G4LorentzVector.hh"
44#include "G4LorentzRotation.hh"
45#include "G4ios.hh"
46
47
48G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(G4int Verbose) : 
49                          G4VDecayChannel("Phase Space", Verbose),
50                          theDaughterMasses(0)
51{
52  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
53}
54
55G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
56                                                   G4double        theBR,
57                                                   G4int           theNumberOfDaughters,
58                                                   const G4String& theDaughterName1,
59                                                   const G4String& theDaughterName2,
60                                                   const G4String& theDaughterName3) :
61                                   G4VDecayChannel("Phase Space",
62                                                   theParentName,theBR,
63                                                   theNumberOfDaughters,
64                                                   theDaughterName1,
65                                                   theDaughterName2,
66                                                   theDaughterName3),
67                                   theDaughterMasses(0)
68{
69  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
70 
71  //   Set the parent particle (resonance) mass to the (default) PDG vale
72  if (parent != NULL)
73     {
74      parentmass = parent->GetPDGMass();
75     }
76}
77
78G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
79                                                   G4double        theParentMass,
80                                                   G4double        theBR,
81                                                   G4int           theNumberOfDaughters,
82                                                   const G4String& theDaughterName1,
83                                                   const G4String& theDaughterName2,
84                                                   const G4String& theDaughterName3) :
85                                   G4VDecayChannel("Phase Space",
86                                                   theParentName,theBR,
87                                                   theNumberOfDaughters,
88                                                   theDaughterName1,
89                                                   theDaughterName2,
90                                                   theDaughterName3),
91                                   parentmass(theParentMass),
92                                   theDaughterMasses(0)
93{
94  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
95}
96
97G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
98                                                   G4double        theParentMass,
99                                                   G4double        theBR,
100                                                   G4int           theNumberOfDaughters,
101                                                   const G4String& theDaughterName1,
102                                                   const G4String& theDaughterName2,
103                                                   const G4String& theDaughterName3,
104                                                   const G4double *masses) :
105                                   G4VDecayChannel("Phase Space",
106                                                   theParentName,theBR,
107                                                   theNumberOfDaughters,
108                                                   theDaughterName1,
109                                                   theDaughterName2,
110                                                   theDaughterName3),
111                                   parentmass(theParentMass),
112                                   theDaughterMasses(masses)
113{
114  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
115}
116
117G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay()
118{
119}
120
121G4DecayProducts *G4GeneralPhaseSpaceDecay::DecayIt(G4double) 
122{
123  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
124  G4DecayProducts * products = NULL;
125 
126  if (parent == NULL) FillParent(); 
127  if (daughters == NULL) FillDaughters();
128 
129  switch (numberOfDaughters){
130  case 0:
131    if (GetVerboseLevel()>0) {
132      G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
133      G4cout << " daughters not defined " <<G4endl;
134    }
135    break;
136  case 1:
137    products =  OneBodyDecayIt();
138    break;
139  case 2:
140    products =  TwoBodyDecayIt();
141    break;
142  case 3:
143    products =  ThreeBodyDecayIt();
144    break;
145  default:
146    products =  ManyBodyDecayIt();
147    break;
148  }
149  if ((products == NULL) && (GetVerboseLevel()>0)) {
150    G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
151    G4cout << *parent_name << " can not decay " << G4endl;
152    DumpInfo();
153  }
154  return products;
155}
156
157G4DecayProducts *G4GeneralPhaseSpaceDecay::OneBodyDecayIt()
158{
159  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt()"<<G4endl;
160
161//  G4double daughtermass = daughters[0]->GetPDGMass();
162
163  //create parent G4DynamicParticle at rest
164  G4ParticleMomentum dummy;
165  G4DynamicParticle * parentparticle = new G4DynamicParticle(parent, dummy, 0.0);
166
167  //create G4Decayproducts
168  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
169  delete parentparticle;
170
171  //create daughter G4DynamicParticle at rest
172  G4DynamicParticle * daughterparticle = new G4DynamicParticle(daughters[0], dummy, 0.0);
173  products->PushProducts(daughterparticle);
174
175  if (GetVerboseLevel()>1) 
176    {
177     G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt ";
178     G4cout << "  create decay products in rest frame " <<G4endl;
179     products->DumpInfo();
180    }
181  return products;
182}
183
184G4DecayProducts *G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()
185{
186  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()"<<G4endl;
187 
188  //daughters'mass
189  G4double daughtermass[2]; 
190  G4double daughtermomentum;
191  if ( theDaughterMasses )
192  { 
193     daughtermass[0]= *(theDaughterMasses);
194     daughtermass[1] = *(theDaughterMasses+1);
195  } else {   
196     daughtermass[0] = daughters[0]->GetPDGMass();
197     daughtermass[1] = daughters[1]->GetPDGMass();
198  }
199 
200//  G4double sumofdaughtermass =  daughtermass[0] + daughtermass[1];
201
202  //create parent G4DynamicParticle at rest
203  G4ParticleMomentum dummy;
204  G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
205
206  //create G4Decayproducts  @@GF why dummy parentparticle?
207  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
208  delete parentparticle;
209
210  //calculate daughter momentum
211  daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
212  G4double costheta = 2.*G4UniformRand()-1.0;
213  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
214  G4double phi  = twopi*G4UniformRand()*rad;
215  G4ParticleMomentum direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
216
217  //create daughter G4DynamicParticle
218  G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum); 
219  G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0],Etotal, direction*daughtermomentum);
220  products->PushProducts(daughterparticle);
221  Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum);
222  daughterparticle = new G4DynamicParticle( daughters[1],Etotal, direction*(-1.0*daughtermomentum));
223  products->PushProducts(daughterparticle);
224
225  if (GetVerboseLevel()>1) 
226    {
227     G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ";
228     G4cout << "  create decay products in rest frame " <<G4endl;
229     products->DumpInfo();
230    }
231  return products;
232}
233
234G4DecayProducts *G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()
235// algorism of this code is originally written in GDECA3 of GEANT3
236{
237  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()"<<G4endl;
238
239  //daughters'mass
240  G4double daughtermass[3]; 
241  G4double sumofdaughtermass = 0.0;
242  for (G4int index=0; index<3; index++)
243    {
244     if ( theDaughterMasses )
245     { 
246         daughtermass[index]= *(theDaughterMasses+index);
247     } else {   
248         daughtermass[index] = daughters[index]->GetPDGMass();
249     }   
250     sumofdaughtermass += daughtermass[index];
251    }
252 
253  //create parent G4DynamicParticle at rest
254  G4ParticleMomentum dummy;
255  G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
256
257  //create G4Decayproducts
258  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
259  delete parentparticle;
260
261  //calculate daughter momentum
262  //  Generate two
263  G4double rd1, rd2, rd;
264  G4double daughtermomentum[3];
265  G4double momentummax=0.0, momentumsum = 0.0;
266  G4double energy;
267
268  do 
269    {
270     rd1 = G4UniformRand();
271     rd2 = G4UniformRand();
272     if (rd2 > rd1) 
273       {
274        rd  = rd1;
275        rd1 = rd2;
276        rd2 = rd;
277       } 
278     momentummax = 0.0;
279     momentumsum = 0.0;
280     // daughter 0
281
282     energy = rd2*(parentmass - sumofdaughtermass);
283     daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
284     if ( daughtermomentum[0] >momentummax )momentummax =  daughtermomentum[0];
285     momentumsum  +=  daughtermomentum[0];
286
287     // daughter 1
288     energy = (1.-rd1)*(parentmass - sumofdaughtermass);
289     daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
290     if ( daughtermomentum[1] >momentummax )momentummax =  daughtermomentum[1];
291     momentumsum  +=  daughtermomentum[1];
292
293     // daughter 2
294     energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
295     daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
296     if ( daughtermomentum[2] >momentummax )momentummax =  daughtermomentum[2];
297     momentumsum  +=  daughtermomentum[2];
298    } while (momentummax >  momentumsum - momentummax );
299
300  // output message
301  if (GetVerboseLevel()>1) {
302    G4cout << "     daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
303    G4cout << "     daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
304    G4cout << "     daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
305    G4cout << "   momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
306  }
307
308  //create daughter G4DynamicParticle
309  G4double costheta, sintheta, phi, sinphi, cosphi; 
310  G4double costhetan, sinthetan, phin, sinphin, cosphin; 
311  costheta = 2.*G4UniformRand()-1.0;
312  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
313  phi  = twopi*G4UniformRand()*rad;
314  sinphi = std::sin(phi);
315  cosphi = std::cos(phi);
316  G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
317  G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]);
318  G4DynamicParticle * daughterparticle
319         = new G4DynamicParticle( daughters[0], Etotal, direction0*daughtermomentum[0]);
320  products->PushProducts(daughterparticle);
321
322  costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
323  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
324  phin  = twopi*G4UniformRand()*rad;
325  sinphin = std::sin(phin);
326  cosphin = std::cos(phin);
327  G4ParticleMomentum direction2;
328  direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 
329  direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 
330  direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
331  Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.mag2());
332  daughterparticle = new G4DynamicParticle( daughters[2],Etotal, direction2*(daughtermomentum[2]/direction2.mag()));
333  products->PushProducts(daughterparticle);
334  G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
335  Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.mag2() );
336  daughterparticle = 
337       new G4DynamicParticle(daughters[1], Etotal, mom);
338  products->PushProducts(daughterparticle);
339
340  if (GetVerboseLevel()>1) {
341     G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ";
342     G4cout << "  create decay products in rest frame " <<G4endl;
343     products->DumpInfo();
344  }
345  return products;
346}
347
348G4DecayProducts *G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()
349// algorism of this code is originally written in FORTRAN by M.Asai
350//*****************************************************************
351//  NBODY
352//   N-body phase space Monte-Carlo generator
353//  Makoto Asai
354//   Hiroshima Institute of Technology
355//    (asai@kekvax.kek.jp)
356//  Revised release : 19/Apr/1995
357//
358{
359  //return value
360  G4DecayProducts *products;
361
362  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()"<<G4endl;
363
364  //daughters'mass
365  G4double *daughtermass = new G4double[numberOfDaughters]; 
366  G4double sumofdaughtermass = 0.0;
367  for (G4int index=0; index<numberOfDaughters; index++){
368    daughtermass[index] = daughters[index]->GetPDGMass();
369    sumofdaughtermass += daughtermass[index];
370  }
371 
372  //Calculate daughter momentum
373  G4double *daughtermomentum = new G4double[numberOfDaughters];
374  G4ParticleMomentum direction; 
375  G4DynamicParticle **daughterparticle;
376  G4double *sm = new G4double[numberOfDaughters];
377  G4double tmas;
378  G4double weight = 1.0;
379  G4int    numberOfTry = 0; 
380  G4int index1, index2;
381
382  do {
383    //Generate rundom number in descending order
384    G4double temp;
385    G4double *rd = new G4double[numberOfDaughters];
386    rd[0] = 1.0;
387    for(index1 =1; index1 < numberOfDaughters -1; index1++)
388      rd[index1] = G4UniformRand(); 
389    rd[ numberOfDaughters -1] = 0.0;
390    for(index1 =1; index1 < numberOfDaughters -1; index1++) {
391      for(index2 = index1+1; index2 < numberOfDaughters; index2++) {
392        if (rd[index1] < rd[index2]){
393          temp         = rd[index1];
394          rd[index1]   = rd[index2];
395          rd[index2]   = temp;
396        }
397      }
398    }
399   
400    //calcurate virtual mass
401    tmas = parentmass -  sumofdaughtermass;
402    temp =  sumofdaughtermass; 
403    for(index1 =0; index1 < numberOfDaughters; index1++) {
404      sm[index1] = rd[index1]*tmas + temp;
405      temp -= daughtermass[index1];
406      if (GetVerboseLevel()>1) {
407        G4cout << index1 << "  rundom number:" << rd[index1]; 
408        G4cout << "   virtual mass:" << sm[index1]/GeV << "[GeV/c/c]" <<G4endl; 
409      }
410    }
411    delete [] rd;
412
413    //Calculate daughter momentum
414    weight = 1.0;
415    index1 =numberOfDaughters-1;
416    daughtermomentum[index1]= Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]);
417    if (GetVerboseLevel()>1) {
418      G4cout << "     daughter " << index1 << ":" << *daughters_name[index1];
419      G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
420    }
421    for(index1 =numberOfDaughters-2; index1>=0; index1--) {
422      // calculate
423      daughtermomentum[index1]= Pmx( sm[index1],daughtermass[index1], sm[index1 +1]);
424      if(daughtermomentum[index1] < 0.0) {
425        // !!! illegal momentum !!!
426        if (GetVerboseLevel()>0) {
427          G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
428          G4cout << "     can not calculate daughter momentum " <<G4endl;
429          G4cout << "     parent:" << *parent_name;
430          G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
431          G4cout << "     daughter " << index1 << ":" << *daughters_name[index1];
432          G4cout << " mass:" << daughtermass[index1]/GeV << "[GeV/c/c]" ;
433          G4cout << " mass:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
434        }
435        delete [] sm;
436        delete [] daughtermass;
437        delete [] daughtermomentum;
438        return NULL;   // Error detection
439
440      } else {
441        // calculate weight of this events
442        weight *=  daughtermomentum[index1]/sm[index1];
443        if (GetVerboseLevel()>1) {
444          G4cout << "     daughter " << index1 << ":" << *daughters_name[index1];
445          G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
446        }
447      }
448    }
449    if (GetVerboseLevel()>1) {
450      G4cout << "    weight: " << weight <<G4endl;
451    }
452   
453    // exit if number of Try exceeds 100
454    if (numberOfTry++ >100) {
455      if (GetVerboseLevel()>0) {
456        G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: ";
457        G4cout << " can not determine Decay Kinematics " << G4endl;
458      }
459      delete [] sm;
460      delete [] daughtermass;
461      delete [] daughtermomentum;
462      return NULL;  // Error detection
463    }
464  } while ( weight > G4UniformRand());
465  if (GetVerboseLevel()>1) {
466      G4cout << "Start calulation of daughters momentum vector "<<G4endl;
467  }
468 
469  G4double costheta, sintheta, phi;
470  G4double beta;
471  daughterparticle = new G4DynamicParticle*[numberOfDaughters];
472
473  index1 = numberOfDaughters -2;
474  costheta = 2.*G4UniformRand()-1.0;
475  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
476  phi  = twopi*G4UniformRand()*rad;
477  direction.setZ(costheta);
478  direction.setY(sintheta*std::sin(phi));
479  direction.setX(sintheta*std::cos(phi));
480  daughterparticle[index1] = new G4DynamicParticle( daughters[index1], direction*daughtermomentum[index1] );
481  daughterparticle[index1+1] = new G4DynamicParticle( daughters[index1+1], direction*(-1.0*daughtermomentum[index1]) );
482
483  for (index1 = numberOfDaughters -3;  index1 >= 0; index1--) {
484    //calculate momentum direction
485    costheta = 2.*G4UniformRand()-1.0;
486    sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
487    phi  = twopi*G4UniformRand()*rad;
488    direction.setZ(costheta);
489    direction.setY(sintheta*std::sin(phi));
490    direction.setX(sintheta*std::cos(phi));
491
492    // boost already created particles
493    beta = daughtermomentum[index1];
494    beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] );
495    for (G4int index2 = index1+1; index2<numberOfDaughters; index2++) {
496      G4LorentzVector p4;
497      // make G4LorentzVector for secondaries
498      p4 = daughterparticle[index2]->Get4Momentum();
499
500      // boost secondaries to  new frame
501      p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
502
503      // change energy/momentum
504      daughterparticle[index2]->Set4Momentum(p4);
505    }
506    //create daughter G4DynamicParticle
507    daughterparticle[index1]= new G4DynamicParticle( daughters[index1], direction*(-1.0*daughtermomentum[index1]));
508  }
509
510  //create G4Decayproducts
511  G4DynamicParticle *parentparticle; 
512  direction.setX(1.0);  direction.setY(0.0); direction.setZ(0.0);
513  parentparticle = new G4DynamicParticle( parent, direction, 0.0);
514  products = new G4DecayProducts(*parentparticle);
515  delete parentparticle;
516  for (index1 = 0; index1<numberOfDaughters; index1++) {
517    products->PushProducts(daughterparticle[index1]);
518  }
519  if (GetVerboseLevel()>1) { 
520    G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
521    G4cout << "  create decay products in rest frame " << G4endl;
522    products->DumpInfo();
523  }
524
525  delete [] daughterparticle;
526  delete [] daughtermomentum;
527  delete [] daughtermass;
528  delete [] sm;
529 
530  return products;
531}
532
533
534
535
536
Note: See TracBrowser for help on using the repository browser.