G4eeTo3PiModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4eeTo3PiModel.cc 66996 2013-01-29 14:50:52Z gcosmo $
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class header file
00031 //
00032 //
00033 // File name:     G4eeTo3PiModel
00034 //
00035 // Author:        Vladimir Ivanchenko
00036 //
00037 // Creation date: 25.10.2003
00038 //
00039 // Modifications:
00040 //
00041 //
00042 // -------------------------------------------------------------------
00043 //
00044 
00045 
00046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00048 
00049 #include "G4eeTo3PiModel.hh"
00050 #include "Randomize.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4PionPlus.hh"
00053 #include "G4PionMinus.hh"
00054 #include "G4PionZero.hh"
00055 #include "G4DynamicParticle.hh"
00056 #include "G4PhysicsVector.hh"
00057 #include "G4PhysicsLinearVector.hh"
00058 #include "G4eeCrossSections.hh"
00059 #include "G4RandomDirection.hh"
00060 #include <complex>
00061 
00062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00063 
00064 using namespace std;
00065 
00066 G4eeTo3PiModel::G4eeTo3PiModel(G4eeCrossSections* cr):
00067   cross(cr)
00068 {
00069   massPi  = G4PionPlus::PionPlus()->GetPDGMass();
00070   massPi0 = G4PionZero::PionZero()->GetPDGMass();
00071   massOm  = 782.62*MeV;
00072   massPhi = 1019.46*MeV;
00073   gcash   = 0.0;
00074   gmax    = 1.0;
00075 }
00076 
00077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00078 
00079 G4eeTo3PiModel::~G4eeTo3PiModel()
00080 {}
00081 
00082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00083 
00084 G4double G4eeTo3PiModel::ThresholdEnergy() const
00085 {
00086   return std::max(LowEnergy(),2.0*massPi + massPi0);
00087 }
00088 
00089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00090 
00091 G4double G4eeTo3PiModel::PeakEnergy() const
00092 {
00093   G4double e = massOm;
00094   if(HighEnergy() > massPhi) e = massPhi; 
00095   return e;
00096 }
00097 
00098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00099 
00100 G4double G4eeTo3PiModel::ComputeCrossSection(G4double e) const
00101 {
00102   G4double ee = std::min(HighEnergy(),e);
00103   return cross->CrossSection3pi(ee);
00104 }
00105 
00106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00107 
00108 G4PhysicsVector* G4eeTo3PiModel::PhysicsVector(G4double emin, 
00109                                                G4double emax) const
00110 {
00111   G4double tmin = std::max(emin, ThresholdEnergy());
00112   G4double tmax = std::max(tmin, emax);
00113   G4int nbins = (G4int)((tmax - tmin)/(1.*MeV));
00114   G4PhysicsVector* v = new G4PhysicsLinearVector(emin,emax,nbins);
00115   v->SetSpline(true);
00116   return v;
00117 }
00118 
00119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00120 
00121 void G4eeTo3PiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
00122             G4double e, const G4ThreeVector& direction) 
00123 {
00124   if(e < ThresholdEnergy()) return;
00125 
00126   G4double x0 = massPi0/e;
00127   G4double x1 = massPi/e;
00128 
00129   G4LorentzVector w0, w1, w2;
00130   G4ThreeVector dir0, dir1;
00131   G4double e0, p0, e2, p, gg, m01, m02, m12;
00132 
00133   // max pi0 energy
00134   G4double edel  = 0.5*e*(1.0 + x0*x0 - 4.0*x1*x1) - massPi0;
00135 
00136   do {
00137     // pi0 sample
00138     e0 = edel*G4UniformRand() + massPi0;
00139     p0 = sqrt(e0 - massPi0*massPi0);
00140     dir0 = G4RandomDirection();
00141     w0 = G4LorentzVector(p0*dir0.x(),p0*dir0.y(),p0*dir0.z(),e0);
00142 
00143     // pi+pi- pair
00144     w1 = G4LorentzVector(-p0*dir0.x(),-p0*dir0.y(),-p0*dir0.z(),e-e0);
00145     G4ThreeVector bst = w1.boostVector();
00146     e2 = 0.25*w1.m2();
00147 
00148     // pi+ 
00149     p = sqrt(e2 - massPi*massPi);
00150     dir1 = G4RandomDirection();
00151     w2 = G4LorentzVector(p*dir1.x(),p*dir1.y(),p*dir1.z(),sqrt(e2));
00152     w2.boost(bst);
00153     G4double px2 = w2.x();
00154     G4double py2 = w2.y();
00155     G4double pz2 = w2.z();
00156 
00157     // pi- 
00158     w1 -= w2;
00159     G4double px1 = w1.x();
00160     G4double py1 = w1.y();
00161     G4double pz1 = w1.z();
00162 
00163     m01 = w0*w1;
00164     m02 = w0*w2;
00165     m12 = w1*w2;
00166 
00167     G4double px = py1*pz2 - py2*pz1;
00168     G4double py = pz1*px2 - pz2*px1;
00169     G4double pz = px1*py2 - px2*py1;
00170 
00171     gg = (px*px + py*py + pz*pz)*
00172       norm( 1.0/cross->DpRho(m01) +  1.0/cross->DpRho(m02)
00173             + 1.0/cross->DpRho(m12) );
00174     if(gg > gmax) {
00175       G4cout << "G4eeTo3PiModel::SampleSecondaries WARNING matrix element g= "
00176              << gg << " > " << gmax << " (majoranta)" << G4endl;
00177     }
00178     if(gg > gcash) gcash = gg;
00179     
00180   } while( gmax*G4UniformRand() > gg );
00181 
00182   w0.rotateUz(direction);
00183   w1.rotateUz(direction);
00184   w2.rotateUz(direction);
00185 
00186   // create G4DynamicParticle objects 
00187   G4DynamicParticle* dp0 = 
00188      new G4DynamicParticle(G4PionZero::PionZero(), w0);
00189   G4DynamicParticle* dp1 = 
00190      new G4DynamicParticle(G4PionPlus::PionPlus(), w1);
00191   G4DynamicParticle* dp2 = 
00192      new G4DynamicParticle(G4PionMinus::PionMinus(), w2);
00193   newp->push_back(dp0);
00194   newp->push_back(dp1);
00195   newp->push_back(dp2);
00196 }
00197 
00198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00199 

Generated on Mon May 27 17:48:05 2013 for Geant4 by  doxygen 1.4.7