ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CollisionComposite.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4CollisionComposite.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 #include "globals.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4CollisionComposite.hh"
31 #include "G4VCollision.hh"
32 #include "G4CollisionVector.hh"
33 #include "G4KineticTrack.hh"
34 #include "G4KineticTrackVector.hh"
35 #include "G4VCrossSectionSource.hh"
36 #include "G4HadTmpUtil.hh"
37 #include "G4AutoLock.hh"
38 
40 
41 const G4double G4CollisionComposite::theT[nPoints] =
42 {.01, .03, .05, .1, .15, .2, .3, .4, .5, .6, .7, .8, .9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 8.0, 10., 15, 20, 50, 100};
43 
45 {
47 }
48 
49 
51 {
53  std::for_each(components.begin(), components.end(), G4Delete());
54 }
55 
56 
58  const G4KineticTrack& trk2) const
59 {
60  G4double crossSect = 0.;
62  if (xSource != 0)
63  // There is a total cross section for this Collision
64  {
65  crossSect = xSource->CrossSection(trk1,trk2);
66  }
67  else
68  {
70  // waiting for mutable to enable buffering.
71  const_cast<G4CollisionComposite *>(this)->BufferCrossSection(trk1.GetDefinition(), trk2.GetDefinition());
72 // G4cerr << "Buffer filled, reying with sqrts = "<< (trk1.Get4Momentum()+trk2.Get4Momentum()).mag() <<G4endl;
73  crossSect = BufferedCrossSection(trk1,trk2);
74  }
75  return crossSect;
76 }
77 
78 
80  const G4KineticTrack& trk2) const
81 {
82  std::vector<G4double> cxCache;
83  G4double partialCxSum = 0.0;
84 
85  size_t i;
86  for (i=0; i<components.size(); i++)
87  {
88  G4double partialCx;
89 // cout << "comp" << i << " " << components[i]()->GetName();
90  if (components[i]->IsInCharge(trk1,trk2))
91  {
92  partialCx = components[i]->CrossSection(trk1,trk2);
93  }
94  else
95  {
96  partialCx = 0.0;
97  }
98 // cout << " cx=" << partialCx << endl;
99  partialCxSum += partialCx;
100  cxCache.push_back(partialCx);
101  }
102 
103  G4double random = G4UniformRand()*partialCxSum;
104  G4double running = 0;
105  for (i=0; i<cxCache.size(); i++)
106  {
107  running += cxCache[i];
108  if (running > random)
109  {
110  return components[i]->FinalState(trk1, trk2);
111  }
112  }
113 // G4cerr <<"in charge = "<<IsInCharge(trk1, trk2)<<G4endl;
114 // G4cerr <<"Cross-section = "<<CrossSection(trk1, trk2)/millibarn<<" "<<running<<" "<<cxCache.size()<<G4endl;
115 // G4cerr <<"Names = "<<trk1.GetDefinition()->GetParticleName()<<", "<<trk2.GetDefinition()->GetParticleName()<<G4endl;
116 // throw G4HadronicException(__FILE__, __LINE__, "G4CollisionComposite: no final state found!");
117  return NULL;
118 }
119 
120 
122  const G4KineticTrack& trk2) const
123 {
124  G4bool isInCharge = false;
125 
126  // The composite is in charge if any of its components is in charge
127 
128  const G4CollisionVector* comps = GetComponents();
129  if (comps)
130  {
131  G4CollisionVector::const_iterator iter;
132  for (iter = comps->begin(); iter != comps->end(); ++iter)
133  {
134  if ( ((*iter))->IsInCharge(trk1,trk2) ) isInCharge = true;
135  }
136  }
137 
138  return isInCharge;
139 }
140 
143 {
144  // check if already buffered
145  size_t i;
146  for(i=0; i<theBuffer.size(); i++)
147  {
148  if(theBuffer[i].InCharge(aP, bP)) return;
149  }
150 // G4cerr << "Buffering for "<<aP->GetParticleName()<<" "<<bP->GetParticleName()<<G4endl;
151 
152  // buffer the new one.
153  G4CrossSectionBuffer aNewBuff(aP, bP);
154  size_t maxE=nPoints;
155  for(size_t tt=0; tt<maxE; tt++)
156  {
157  G4double aT = theT[tt]*GeV;
158  G4double crossSect = 0;
159  // The total cross-section is summed over all the component channels
160 
161  //A.R. 28-Sep-2012 Fix reproducibility problem
162  // Assign the kinetic energy to the lightest of the
163  // two particles, instead to the first one always.
164  G4double atime = 0;
165  G4double btime = 0;
166  G4ThreeVector aPosition(0,0,0);
167  G4ThreeVector bPosition(0,0,0);
168  G4double aM = aP->GetPDGMass();
169  G4double bM = bP->GetPDGMass();
170  G4double aE = aM;
171  G4double bE = bM;
172  G4ThreeVector aMom(0,0,0);
173  G4ThreeVector bMom(0,0,0);
174  if ( aM <= bM ) {
175  aE += aT;
176  aMom = G4ThreeVector(0,0,std::sqrt(aE*aE-aM*aM));
177  } else {
178  bE += aT;
179  bMom = G4ThreeVector(0,0,std::sqrt(bE*bE-bM*bM));
180  }
181  G4LorentzVector a4Momentum(aE, aMom);
182  G4LorentzVector b4Momentum(bE, bMom);
183  G4KineticTrack a(aP, atime, aPosition, a4Momentum);
184  G4KineticTrack b(bP, btime, bPosition, b4Momentum);
185 
186  for (i=0; i<components.size(); i++)
187  {
188  if(components[i]->IsInCharge(a,b))
189  {
190  crossSect += components[i]->CrossSection(a,b);
191  }
192  }
193  G4double sqrts = (a4Momentum+b4Momentum).mag();
194  aNewBuff.push_back(sqrts, crossSect);
195  }
196  theBuffer.push_back(aNewBuff);
197 // theBuffer.back().Print();
198 }
199 
200 
202 BufferedCrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
203 {
204  for(size_t i=0; i<theBuffer.size(); i++)
205  {
206  if(theBuffer[i].InCharge(trk1.GetDefinition(), trk2.GetDefinition()))
207  {
208  return theBuffer[i].CrossSection(trk1, trk2);
209  }
210  }
211  throw G4HadronicException(__FILE__, __LINE__, "G4CollisionComposite::BufferedCrossSection - Blitz !!");
212  return 0;
213 }
214