ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PionDecayMakeSpin.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PionDecayMakeSpin.cc
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 = nullptr;
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  G4DynamicParticle* aSecondary = (*products)[index];
97  const G4ParticleDefinition* aSecondaryDef = aSecondary->GetDefinition();
98 
99  if (aSecondaryDef == aMuonPlus ||
100  aSecondaryDef == aMuonMinus ) {
101  // Muon+ or Muon-
102  aMuon = aSecondary;
103  emu = aSecondary->GetTotalEnergy();
104  p_muon = aSecondary->GetMomentum();
105  } else if (aSecondaryDef == aNeutrinoMu ||
106  aSecondaryDef == aAntiNeutrinoMu ) {
107  // Muon-Neutrino / Muon-Anti-Neutrino
108  eneutrino = aSecondary->GetTotalEnergy();
109  p_neutrino = aSecondary->GetMomentum();
110  }
111  }
112  }
113 
114  // This routine deals only with decays with a
115  // muon and mu-(anti)neutrinos in the final state
116  if (aMuon == nullptr) return;
117  if (eneutrino==0||emu==0) return;
118 
119  G4ThreeVector spin(0,0,0);
120 
121  const G4DynamicParticle* theParentParticle = products->GetParentParticle();
122 
123  G4double amass = theParentParticle->GetMass();
124  G4double emmu = aMuonPlus->GetPDGMass();
125 
126  if (numberOfSecondaries == 2 ) {
127  G4double scale = - (eneutrino - ( p_muon * p_neutrino )/(emu+emmu));
128 
129  p_muon = scale * p_muon;
130  p_neutrino = emmu * p_neutrino;
131  spin = p_muon + p_neutrino;
132 
133  scale = 2./(amass*amass-emmu*emmu);
134  spin = scale * spin;
135 
136  if (aParticle->GetCharge() < 0.0) spin = -spin;
137 
138  } else {
139  spin = G4RandomDirection();
140 
141  }
142 
143  spin = spin.unit();
144 
145  aMuon->SetPolarization(spin.x(),spin.y(),spin.z());
146 
147  return;
148 }
149 
150 void G4PionDecayMakeSpin::ProcessDescription(std::ostream& outFile) const
151 {
152  outFile << GetProcessName()
153  << ": Decay of mesons that can decay into a muon \n"
154  << " i.e. pi+, pi-, K+, K- and K0_long \n"
155  << " kinematics of daughters are dertermined by DecayChannels \n"
156  << " polarization of daughter particles are take into account. \n";
157 }
158