G4TauLeptonicDecayChannel.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 //
00027 // $Id$
00028 //
00029 // 
00030 // ------------------------------------------------------------
00031 //      GEANT 4 class header file
00032 //
00033 //      History: first implementation, based on object model of
00034 //      30 May  1997 H.Kurashige
00035 //
00036 //      Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige 
00037 // ------------------------------------------------------------
00038 
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DecayProducts.hh"
00043 #include "G4VDecayChannel.hh"
00044 #include "G4TauLeptonicDecayChannel.hh"
00045 #include "Randomize.hh"
00046 #include "G4LorentzVector.hh"
00047 #include "G4LorentzRotation.hh"
00048 
00049 
00050 G4TauLeptonicDecayChannel::G4TauLeptonicDecayChannel()
00051   :G4VDecayChannel()
00052 {
00053 }
00054 
00055 
00056 G4TauLeptonicDecayChannel::G4TauLeptonicDecayChannel(
00057                                                      const G4String& theParentName, 
00058                                                      G4double        theBR,
00059                                                      const G4String& theLeptonName)
00060                    :G4VDecayChannel("Tau Leptonic Decay",1)
00061 {
00062   // set names for daughter particles
00063   if (theParentName == "tau+") {
00064     SetBR(theBR);
00065     SetParent("tau+");
00066     SetNumberOfDaughters(3);
00067     if ((theLeptonName=="e-"||theLeptonName=="e+")){
00068       SetDaughter(0, "e+");
00069       SetDaughter(1, "nu_e");
00070       SetDaughter(2, "anti_nu_tau");
00071     } else {
00072       SetDaughter(0, "mu+");
00073       SetDaughter(1, "nu_mu");
00074       SetDaughter(2, "anti_nu_tau");
00075     } 
00076   } else if (theParentName == "tau-") {
00077     SetBR(theBR);
00078     SetParent("tau-");
00079     SetNumberOfDaughters(3);
00080     if ((theLeptonName=="e-"||theLeptonName=="e+")){
00081       SetDaughter(0, "e-");
00082       SetDaughter(1, "anti_nu_e");
00083       SetDaughter(2, "nu_tau");
00084     } else {
00085       SetDaughter(0, "mu-");
00086       SetDaughter(1, "anti_nu_mu");
00087       SetDaughter(2, "nu_tau");
00088     } 
00089   } else {
00090 #ifdef G4VERBOSE
00091     if (GetVerboseLevel()>0) {
00092       G4cout << "G4TauLeptonicDecayChannel:: constructor :";
00093       G4cout << " parent particle is not tau but ";
00094       G4cout << theParentName << G4endl;
00095     }
00096 #endif
00097   }
00098 }
00099 
00100 G4TauLeptonicDecayChannel::~G4TauLeptonicDecayChannel()
00101 {
00102 }
00103 
00104 G4TauLeptonicDecayChannel::G4TauLeptonicDecayChannel(const G4TauLeptonicDecayChannel &right):
00105   G4VDecayChannel(right)
00106 {
00107 }
00108 
00109 G4TauLeptonicDecayChannel & G4TauLeptonicDecayChannel::operator=(const G4TauLeptonicDecayChannel & right)
00110 {
00111   if (this != &right) { 
00112     kinematics_name = right.kinematics_name;
00113     verboseLevel = right.verboseLevel;
00114     rbranch = right.rbranch;
00115 
00116     // copy parent name
00117     parent_name = new G4String(*right.parent_name);
00118 
00119     // clear daughters_name array
00120     ClearDaughtersName();
00121 
00122     // recreate array
00123     numberOfDaughters = right.numberOfDaughters;
00124     if ( numberOfDaughters >0 ) {
00125       if (daughters_name !=0) ClearDaughtersName();
00126       daughters_name = new G4String*[numberOfDaughters];
00127       //copy daughters name
00128       for (G4int index=0; index < numberOfDaughters; index++) {
00129           daughters_name[index] = new G4String(*right.daughters_name[index]);
00130       }
00131     }
00132   }
00133   return *this;
00134 }
00135 
00136 G4DecayProducts *G4TauLeptonicDecayChannel::DecayIt(G4double) 
00137 {
00138   // this version neglects muon polarization 
00139   //              assumes the pure V-A coupling
00140   //              gives incorrect energy spectrum for neutrinos
00141 #ifdef G4VERBOSE
00142   if (GetVerboseLevel()>1) G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
00143 #endif
00144 
00145   if (parent == 0) FillParent();  
00146   if (daughters == 0) FillDaughters();
00147  
00148   // parent mass
00149   G4double parentmass = parent->GetPDGMass();
00150 
00151   //daughters'mass
00152   G4double daughtermass[3]; 
00153   for (G4int index=0; index<3; index++){
00154     daughtermass[index] = daughters[index]->GetPDGMass();
00155   }
00156 
00157    //create parent G4DynamicParticle at rest
00158   G4ThreeVector dummy;
00159   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00160   //create G4Decayproducts
00161   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00162   delete parentparticle;
00163 
00164   // calculate daughter momentum
00165   G4double daughtermomentum[3];
00166 
00167   // calcurate lepton momentum
00168   G4double pmax = (parentmass*parentmass-daughtermass[0]*daughtermass[0])/2./parentmass;
00169   G4double p, e;
00170   G4double r;
00171   do {
00172     // determine momentum/energy
00173     r = G4UniformRand();
00174     p = pmax*G4UniformRand();
00175     e = std::sqrt(p*p + daughtermass[0]*daughtermass[0]);
00176   } while (r > spectrum(p,e,parentmass,daughtermass[0]) );
00177 
00178   //create daughter G4DynamicParticle 
00179   // daughter 0 (lepton)
00180   daughtermomentum[0] = p;
00181   G4double costheta, sintheta, phi, sinphi, cosphi; 
00182   costheta = 2.*G4UniformRand()-1.0;
00183   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00184   phi  = twopi*G4UniformRand()*rad;
00185   sinphi = std::sin(phi);
00186   cosphi = std::cos(phi);
00187   G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
00188   G4DynamicParticle * daughterparticle 
00189          = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
00190   products->PushProducts(daughterparticle);
00191 
00192   // daughter 1 ,2 (nutrinos)
00193   // create neutrinos in the C.M frame of two neutrinos
00194   G4double energy2 = parentmass-e; 
00195   G4double vmass   = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
00196   G4double beta = -1.0*daughtermomentum[0]/energy2;
00197   G4double costhetan = 2.*G4UniformRand()-1.0;
00198   G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00199   G4double phin  = twopi*G4UniformRand()*rad;
00200   G4double sinphin = std::sin(phin);
00201   G4double cosphin = std::cos(phin);
00202 
00203   G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
00204   G4DynamicParticle * daughterparticle1 
00205          = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
00206   G4DynamicParticle * daughterparticle2
00207          = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
00208 
00209   // boost to the muon rest frame
00210   G4LorentzVector p4;
00211   p4 = daughterparticle1->Get4Momentum();
00212   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
00213   daughterparticle1->Set4Momentum(p4);
00214   p4 = daughterparticle2->Get4Momentum();
00215   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
00216   daughterparticle2->Set4Momentum(p4);
00217   products->PushProducts(daughterparticle1);
00218   products->PushProducts(daughterparticle2);
00219   daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
00220   daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
00221 
00222 
00223  // output message
00224 #ifdef G4VERBOSE
00225   if (GetVerboseLevel()>1) {
00226     G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
00227     G4cout << "  create decay products in rest frame " <<G4endl;
00228     products->DumpInfo();
00229   }
00230 #endif
00231   return products;
00232 }
00233 
00234 
00235 
00236 
00237 G4double G4TauLeptonicDecayChannel::spectrum(G4double p,
00238                                              G4double e,
00239                                              G4double mtau,
00240                                              G4double ml)
00241 {
00242   G4double f1;
00243   f1 = 3.0*e*(mtau*mtau+ml*ml)-4.0*mtau*e*e-2.0*mtau*ml*ml;
00244   return p*(f1)/(mtau*mtau*mtau*mtau)/(0.6);
00245 }
00246 
00247 
00248 

Generated on Mon May 27 17:49:56 2013 for Geant4 by  doxygen 1.4.7