ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLPiNToEtaChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLPiNToEtaChannel.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 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 #include "G4INCLPiNToEtaChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 
45 namespace G4INCL {
46 
47  PiNToEtaChannel::PiNToEtaChannel(Particle *p1, Particle *p2)
48  : particle1(p1), particle2(p2)
49  {
50 
51  }
52 
53  PiNToEtaChannel::~PiNToEtaChannel(){
54 
55  }
56 
57  void PiNToEtaChannel::fillFinalState(FinalState *fs) {
58  Particle * nucleon;
59  Particle * pion;
60  if(particle1->isNucleon()) {
61  nucleon = particle1;
62  pion = particle2;
63  } else {
64  nucleon = particle2;
65  pion = particle1;
66  }
67 
68 
69  G4int iso=ParticleTable::getIsospin(nucleon->getType())+ParticleTable::getIsospin(pion->getType());
70 // assert(iso == 1 || iso == -1);
71  if (iso == 1) {
72  nucleon->setType(Proton);
73  }
74  else if (iso == -1) {
75  nucleon->setType(Neutron);
76  }
77  pion->setType(Eta);
78  G4double sh=nucleon->getEnergy()+pion->getEnergy();
79  G4double mn=nucleon->getMass();
80  G4double me=pion->getMass();
81  G4double en=(sh*sh+mn*mn-me*me)/(2*sh);
82  nucleon->setEnergy(en);
83  G4double ee=std::sqrt(en*en-mn*mn+me*me);
84  pion->setEnergy(ee);
85  G4double pn=std::sqrt(en*en-mn*mn);
86 
87 // real distribution (from PRC 78, 025204 (2008))
88 
89 
90  G4double ECM=G4INCL::KinematicsUtils::totalEnergyInCM(particle1,particle2);
91 
92  const G4double pi=std::acos(-1.0);
93  G4double x1;
94  G4double u1;
95  G4double fteta;
96  G4double teta;
97  G4double fi;
98 
99  if (ECM < 1650.) {
100 // below 1650 MeV - angular distribution (x=cos(theta): ax^2+bx+c
101 
102  G4double f1= -0.0000288627*ECM*ECM+0.09155289*ECM-72.25436; // f(1) that is the maximum (fit on experimental data)
103  G4double b1=(f1-(f1/(1.5-0.5*std::pow((ECM-1580.)/95.,2))))/2.; // ideas: 1) f(-1)=0.5f(1); 2) "power term" flattens the distribution away from ECM=1580 MeV
104  G4double a1=2.5*b1; // minimum at cos(theta) = -0.2
105  G4double c1=f1-3.5*b1;
106 
107  G4double interg1=2.*a1/3. +2.*c1; // (integral to normalize)
108 
109  G4int passe1=0;
110  while (passe1==0) {
111  // Sample x from -1 to 1
112  x1=Random::shoot();
113  if (Random::shoot() > 0.5) x1=-x1;
114 
115  // Sample u from 0 to 1
116  u1=Random::shoot();
117  fteta=(a1*x1*x1+b1*x1+c1)/interg1;
118  // The condition
119  if (u1*f1/interg1 < fteta) {
120  teta=std::acos(x1);
121  passe1=1;
122  }
123  }
124  }
125  else {
126 // above 1650 MeV - angular distribution (x=cos(theta): (ax^2+bx+c)*(0.5+(arctan(10*(x+dev)))/pi) + vert
127 
128  G4double a2=-0.29;
129  G4double b2=0.348; // ax^2+bx+c: around cos(theta)=0.6 with maximum at 0.644963 (value = 0.1872666)
130  G4double c2=0.0546;
131  G4double dev=-0.2; // tail close to zero from "dev" down to -1
132  G4double vert=0.04; // to avoid negative differential cross sections
133 
134  G4double interg2=0.1716182902205207; // with the above given parameters! (integral to normalize)
135  const G4double f2=1.09118088; // maximum (integral taken into account)
136 
137  G4int passe2=0;
138  while (passe2==0) {
139  // Sample x from -1 to 1
140  x1=Random::shoot();
141  if (Random::shoot() > 0.5) x1=-x1;
142 
143  // Sample u from 0 to 1
144  u1=Random::shoot();
145  fteta=((a2*x1*x1+b2*x1+c2)*(0.5+(std::atan(10*(x1+dev)))/pi) + vert)/interg2;
146  // The condition
147  if (u1*f2 < fteta) {
148  teta=std::acos(x1);
149  passe2=1;
150  }
151  }
152  }
153 
154  fi=(2.0*pi)*Random::shoot();
155 
156  ThreeVector mom_nucleon(
157  pn*std::sin(teta)*std::cos(fi),
158  pn*std::sin(teta)*std::sin(fi),
159  pn*std::cos(teta)
160  );
161 // end real distribution
162 
163  nucleon->setMomentum(-mom_nucleon);
164  pion->setMomentum(mom_nucleon);
165 
166  fs->addModifiedParticle(nucleon);
167  fs->addModifiedParticle(pion);
168  }
169 
170 }