Geant4_10
G4PionDecayMakeSpin.cc
Go to the documentation of this file.
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 
29 #include "G4PionDecayMakeSpin.hh"
30 
31 #include "G4Decay.hh"
32 #include "G4DecayProducts.hh"
33 
34 #include "G4RandomDirection.hh"
35 
36 // constructor
37 
39  : G4Decay(processName)
40 {
41  // set Process Sub Type
42  SetProcessSubType(static_cast<int>(DECAY_PionMakeSpin));
43 
44 }
45 
47 
49  G4DecayProducts* products)
50 {
51  // This routine deals only with particles that can decay into a muon
52  // pi+, pi-, K+, K- and K0_long
53 
54  // get particle
55 
56  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
57  const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
58 
59  G4ParticleDefinition* aMuonPlus =
61  G4ParticleDefinition* aMuonMinus =
63  G4ParticleDefinition* aPionPlus =
65  G4ParticleDefinition* aPionMinus =
67  G4ParticleDefinition* aKaonPlus =
69  G4ParticleDefinition* aKaonMinus =
71  G4ParticleDefinition* aKaon0Long =
73  G4ParticleDefinition* aNeutrinoMu =
75  G4ParticleDefinition* aAntiNeutrinoMu =
77 
78  if( aParticleDef == aPionPlus ||
79  aParticleDef == aPionMinus ||
80  aParticleDef == aKaonPlus ||
81  aParticleDef == aKaonMinus ||
82  aParticleDef == aKaon0Long ) {
83  } else {
84  return;
85  }
86 
87  G4DynamicParticle* aMuon = NULL;
88 
89  G4double emu(0), eneutrino(0);
90  G4ThreeVector p_muon, p_neutrino;
91 
92  G4int numberOfSecondaries = products->entries();
93 
94  if (numberOfSecondaries > 0) {
95  for (G4int index=0; index < numberOfSecondaries; index++)
96  {
97  G4DynamicParticle* aSecondary = (*products)[index];
98  const G4ParticleDefinition* aSecondaryDef = aSecondary->GetDefinition();
99 
100  if (aSecondaryDef == aMuonPlus ||
101  aSecondaryDef == aMuonMinus ) {
102  // Muon+ or Muon-
103  aMuon = aSecondary;
104  emu = aSecondary->GetTotalEnergy();
105  p_muon = aSecondary->GetMomentum();
106  } else if (aSecondaryDef == aNeutrinoMu ||
107  aSecondaryDef == aAntiNeutrinoMu ) {
108  // Muon-Neutrino / Muon-Anti-Neutrino
109  eneutrino = aSecondary->GetTotalEnergy();
110  p_neutrino = aSecondary->GetMomentum();
111  }
112  }
113  }
114 
115  // This routine deals only with decays with a
116  // muon and mu-(anti)neutrinos in the final state
117 
118  if (!aMuon) return;
119  if (eneutrino==0||emu==0) return;
120 
121  G4ThreeVector spin(0,0,0);
122 
123  const G4DynamicParticle* theParentParticle = products->GetParentParticle();
124 
125  G4double amass = theParentParticle->GetMass();
126  G4double emmu = aMuonPlus->GetPDGMass();
127 
128  if (numberOfSecondaries == 2 ) {
129 
130  G4double scale = - (eneutrino - ( p_muon * p_neutrino )/(emu+emmu));
131 
132  p_muon = scale * p_muon;
133  p_neutrino = emmu * p_neutrino;
134  spin = p_muon + p_neutrino;
135 
136  scale = 2./(amass*amass-emmu*emmu);
137  spin = scale * spin;
138 
139  if (aParticle->GetCharge() < 0.0) spin = -spin;
140 
141  } else {
142 
143  spin = G4RandomDirection();
144 
145  }
146 
147  spin = spin.unit();
148 
149  aMuon->SetPolarization(spin.x(),spin.y(),spin.z());
150 
151  return;
152 }
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Int_t index
Definition: macro.C:9
G4double GetTotalEnergy() const
double x() const
virtual void DaughterPolarization(const G4Track &aTrack, G4DecayProducts *products)
const G4DynamicParticle * GetDynamicParticle() const
G4PionDecayMakeSpin(const G4String &processName="Decay")
G4ThreeVector G4RandomDirection()
const G4DynamicParticle * GetParentParticle() const
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
double z() const
G4double GetMass() const
G4double GetCharge() const
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
void SetPolarization(G4double polX, G4double polY, G4double polZ)
Double_t scale
Definition: plot.C:11
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
Hep3Vector unit() const
double y() const
G4int entries() const
double G4double
Definition: G4Types.hh:76
G4ThreeVector GetMomentum() const