source: trunk/source/event/src/G4SPSRandomGenerator.cc @ 1249

Last change on this file since 1249 was 816, checked in by garnier, 16 years ago

import all except CVS
File size: 21.4 KB
Line 
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// MODULE:        G4SPSRandomGenerator.cc
29//
30// Version:      1.0
31// Date:         5/02/04
32// Author:       Fan Lei
33// Organisation: QinetiQ ltd.
34// Customer:     ESA/ESTEC
35//
36///////////////////////////////////////////////////////////////////////////////
37//
38// CHANGE HISTORY
39// --------------
40//
41//
42// Version 1.0, 05/02/2004, Fan Lei, Created.
43//    Based on the G4GeneralParticleSource class in Geant4 v6.0
44//
45///////////////////////////////////////////////////////////////////////////////
46//
47#include "G4PrimaryParticle.hh"
48#include "G4Event.hh"
49#include "Randomize.hh"
50#include <cmath>
51#include "G4TransportationManager.hh"
52#include "G4VPhysicalVolume.hh"
53#include "G4PhysicalVolumeStore.hh"
54#include "G4ParticleTable.hh"
55#include "G4ParticleDefinition.hh"
56#include "G4IonTable.hh"
57#include "G4Ions.hh"
58#include "G4TrackingManager.hh"
59#include "G4Track.hh"
60
61#include "G4SPSRandomGenerator.hh"
62
63//G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
64
65G4SPSRandomGenerator::G4SPSRandomGenerator()
66{
67  // Initialise all variables
68
69  // Bias variables
70  XBias = false;
71  IPDFXBias = false;
72  YBias = false;
73  IPDFYBias = false;
74  ZBias = false;
75  IPDFZBias = false;
76  ThetaBias = false;
77  IPDFThetaBias = false;
78  PhiBias = false;
79  IPDFPhiBias = false;
80  EnergyBias = false;
81  IPDFEnergyBias = false;
82  PosThetaBias = false;
83  IPDFPosThetaBias = false;
84  PosPhiBias = false;
85  IPDFPosPhiBias = false;
86  bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4] = bweights[5] = bweights[6] = bweights[7]= bweights[8]= 1. ;
87  verbosityLevel = 0 ;
88}
89
90G4SPSRandomGenerator::~G4SPSRandomGenerator()
91{}
92
93//G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
94//{
95//  if (instance == 0) instance = new G4SPSRandomGenerator();
96//  return instance;
97//}
98
99// Biasing methods
100
101void G4SPSRandomGenerator::SetXBias(G4ThreeVector input)
102{ 
103  G4double ehi, val;
104  ehi = input.x();
105  val = input.y();
106  XBiasH.InsertValues(ehi, val);
107  XBias = true;
108}
109
110void G4SPSRandomGenerator::SetYBias(G4ThreeVector input)
111{
112  G4double ehi, val;
113  ehi = input.x();
114  val = input.y();
115  YBiasH.InsertValues(ehi, val);
116  YBias = true;
117}
118
119void G4SPSRandomGenerator::SetZBias(G4ThreeVector input)
120{
121  G4double ehi, val;
122  ehi = input.x();
123  val = input.y();
124  ZBiasH.InsertValues(ehi, val);
125  ZBias = true;
126}
127
128void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input)
129{
130  G4double ehi, val;
131  ehi = input.x();
132  val = input.y();
133  ThetaBiasH.InsertValues(ehi, val);
134  ThetaBias = true;
135}
136
137void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input)
138{
139  G4double ehi, val;
140  ehi = input.x();
141  val = input.y();
142  PhiBiasH.InsertValues(ehi, val);
143  PhiBias = true;
144}
145
146void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input)
147{
148  G4double ehi, val;
149  ehi = input.x();
150  val = input.y();
151  EnergyBiasH.InsertValues(ehi, val);
152  EnergyBias = true;
153}
154
155void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input)
156{
157  G4double ehi, val;
158  ehi = input.x();
159  val = input.y();
160  PosThetaBiasH.InsertValues(ehi, val);
161  PosThetaBias = true;
162}
163
164void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input)
165{
166  G4double ehi, val;
167  ehi = input.x();
168  val = input.y();
169  PosPhiBiasH.InsertValues(ehi, val);
170  PosPhiBias = true;
171}
172
173void G4SPSRandomGenerator::ReSetHist(G4String atype)
174{
175  if ( atype == "biasx") {
176    XBias = false ;
177    IPDFXBias = false;
178    XBiasH = IPDFXBiasH = ZeroPhysVector ;}
179  else if ( atype == "biasy") {
180    YBias = false ;
181    IPDFYBias = false;
182    YBiasH = IPDFYBiasH = ZeroPhysVector ;} 
183  else if ( atype == "biasz") {
184    ZBias = false ;
185    IPDFZBias = false;
186    ZBiasH = IPDFZBiasH = ZeroPhysVector ;}
187  else if ( atype == "biast") {
188    ThetaBias = false ;
189    IPDFThetaBias = false;
190    ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector ;}
191  else if ( atype == "biasp") {
192    PhiBias = false ;
193    IPDFPhiBias = false;
194    PhiBiasH = IPDFPhiBiasH = ZeroPhysVector ;}
195  else if ( atype == "biase") {
196    EnergyBias = false ;
197    IPDFEnergyBias = false;
198    EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector ;}
199  else if ( atype == "biaspt") {
200    PosThetaBias = false ;
201    IPDFPosThetaBias = false;
202    PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector ;}
203  else if ( atype == "biaspp") {
204    PosPhiBias = false ;
205    IPDFPosPhiBias = false;
206    PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector ;}
207  else {
208    G4cout << "Error, histtype not accepted " << G4endl;
209  }
210}
211
212G4double G4SPSRandomGenerator::GenRandX()
213{
214  if(verbosityLevel >= 1)
215    G4cout << "In GenRandX" << G4endl;
216  if(XBias == false)
217    {
218      // X is not biased
219      G4double rndm = G4UniformRand();
220      return(rndm);
221    }
222  else
223    {
224      // X is biased
225      if(IPDFXBias == false)
226        {
227          // IPDF has not been created, so create it
228          G4double bins[1024],vals[1024], sum;
229          G4int ii;
230          G4int maxbin = G4int(XBiasH.GetVectorLength());
231          bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
232          vals[0] = XBiasH(size_t(0));
233          sum = vals[0];
234          for(ii=1;ii<maxbin;ii++)
235            {
236              bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
237              vals[ii] = XBiasH(size_t(ii)) + vals[ii-1];
238              sum = sum + XBiasH(size_t(ii));
239            }
240
241          for(ii=0;ii<maxbin;ii++)
242            {
243              vals[ii] = vals[ii]/sum;
244              IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
245            }
246          // Make IPDFXBias = true
247          IPDFXBias = true;
248        }
249      // IPDF has been create so carry on
250      G4double rndm = G4UniformRand();
251
252      // Calculate the weighting: Find the bin that the determined
253      // rndm is in and the weigthing will be the difference in the
254      // natural probability (from the x-axis) divided by the
255      // difference in the biased probability (the area).
256      size_t numberOfBin = IPDFXBiasH.GetVectorLength();
257      G4int biasn1 = 0;
258      G4int biasn2 = numberOfBin/2;
259      G4int biasn3 = numberOfBin - 1;
260      while (biasn1 != biasn3 - 1) {
261      if (rndm > IPDFXBiasH(biasn2))
262         biasn1 = biasn2;
263      else
264         biasn3 = biasn2;
265      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
266      }
267      // retrieve the areas and then the x-axis values
268      bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
269      G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
270      G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
271      G4double NatProb = xaxisu - xaxisl;
272      //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl;
273      //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl;
274      bweights[0] = NatProb/bweights[0];
275      if(verbosityLevel >= 1)
276        G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
277      return(IPDFXBiasH.GetEnergy(rndm));
278    }
279}
280
281G4double G4SPSRandomGenerator::GenRandY()
282{
283  if(verbosityLevel >= 1)
284    G4cout << "In GenRandY" << G4endl;
285  if(YBias == false)
286    {
287      // Y is not biased
288      G4double rndm = G4UniformRand();
289      return(rndm);
290    }
291  else
292    {
293      // Y is biased
294      if(IPDFYBias == false)
295        {
296          // IPDF has not been created, so create it
297          G4double bins[1024],vals[1024], sum;
298          G4int ii;
299          G4int maxbin = G4int(YBiasH.GetVectorLength());
300          bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
301          vals[0] = YBiasH(size_t(0));
302          sum = vals[0];
303          for(ii=1;ii<maxbin;ii++)
304            {
305              bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
306              vals[ii] = YBiasH(size_t(ii)) + vals[ii-1];
307              sum = sum + YBiasH(size_t(ii));
308            }
309
310          for(ii=0;ii<maxbin;ii++)
311            {
312              vals[ii] = vals[ii]/sum;
313              IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
314            }
315          // Make IPDFYBias = true
316          IPDFYBias = true;
317        }
318      // IPDF has been create so carry on
319      G4double rndm = G4UniformRand();
320      size_t numberOfBin = IPDFYBiasH.GetVectorLength();
321      G4int biasn1 = 0;
322      G4int biasn2 = numberOfBin/2;
323      G4int biasn3 = numberOfBin - 1;
324      while (biasn1 != biasn3 - 1) {
325      if (rndm > IPDFYBiasH(biasn2))
326         biasn1 = biasn2;
327      else
328         biasn3 = biasn2;
329      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
330      }
331      bweights[1] =  IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
332      G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
333      G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
334      G4double NatProb = xaxisu - xaxisl;
335      bweights[1] = NatProb/bweights[1];
336      if(verbosityLevel >= 1)
337        G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
338      return(IPDFYBiasH.GetEnergy(rndm));
339    }
340}
341
342G4double G4SPSRandomGenerator::GenRandZ()
343{
344  if(verbosityLevel >= 1)
345    G4cout << "In GenRandZ" << G4endl;
346  if(ZBias == false)
347    {
348      // Z is not biased
349      G4double rndm = G4UniformRand();
350      return(rndm);
351    }
352  else
353    {
354      // Z is biased
355      if(IPDFZBias == false)
356        {
357          // IPDF has not been created, so create it
358          G4double bins[1024],vals[1024], sum;
359          G4int ii;
360          G4int maxbin = G4int(ZBiasH.GetVectorLength());
361          bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
362          vals[0] = ZBiasH(size_t(0));
363          sum = vals[0];
364          for(ii=1;ii<maxbin;ii++)
365            {
366              bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
367              vals[ii] = ZBiasH(size_t(ii)) + vals[ii-1];
368              sum = sum + ZBiasH(size_t(ii));
369            }
370
371          for(ii=0;ii<maxbin;ii++)
372            {
373              vals[ii] = vals[ii]/sum;
374              IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
375            }
376          // Make IPDFZBias = true
377          IPDFZBias = true;
378        }
379      // IPDF has been create so carry on
380      G4double rndm = G4UniformRand();
381      //      size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm);
382      size_t numberOfBin = IPDFZBiasH.GetVectorLength();
383      G4int biasn1 = 0;
384      G4int biasn2 = numberOfBin/2;
385      G4int biasn3 = numberOfBin - 1;
386      while (biasn1 != biasn3 - 1) {
387      if (rndm > IPDFZBiasH(biasn2))
388         biasn1 = biasn2;
389      else
390         biasn3 = biasn2;
391      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
392      }
393      bweights[2] =  IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
394      G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
395      G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
396      G4double NatProb = xaxisu - xaxisl;
397      bweights[2] = NatProb/bweights[2];
398      if(verbosityLevel >= 1)
399        G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
400      return(IPDFZBiasH.GetEnergy(rndm));
401    }
402}
403
404G4double G4SPSRandomGenerator::GenRandTheta()
405{
406  if(verbosityLevel >= 1)
407    {
408      G4cout << "In GenRandTheta" << G4endl;
409      G4cout << "Verbosity " << verbosityLevel << G4endl;
410    }
411  if(ThetaBias == false)
412    {
413      // Theta is not biased
414      G4double rndm = G4UniformRand();
415      return(rndm);
416    }
417  else
418    {
419      // Theta is biased
420      if(IPDFThetaBias == false)
421        {
422          // IPDF has not been created, so create it
423          G4double bins[1024],vals[1024], sum;
424          G4int ii;
425          G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
426          bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
427          vals[0] = ThetaBiasH(size_t(0));
428          sum = vals[0];
429          for(ii=1;ii<maxbin;ii++)
430            {
431              bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
432              vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii-1];
433              sum = sum + ThetaBiasH(size_t(ii));
434            }
435
436          for(ii=0;ii<maxbin;ii++)
437            {
438              vals[ii] = vals[ii]/sum;
439              IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
440            }
441          // Make IPDFThetaBias = true
442          IPDFThetaBias = true;
443        }
444      // IPDF has been create so carry on
445      G4double rndm = G4UniformRand();
446      //      size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
447      size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
448      G4int biasn1 = 0;
449      G4int biasn2 = numberOfBin/2;
450      G4int biasn3 = numberOfBin - 1;
451      while (biasn1 != biasn3 - 1) {
452      if (rndm > IPDFThetaBiasH(biasn2))
453         biasn1 = biasn2;
454      else
455         biasn3 = biasn2;
456      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
457      }
458      bweights[3] =  IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
459      G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
460      G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
461      G4double NatProb = xaxisu - xaxisl;
462      bweights[3] = NatProb/bweights[3];
463      if(verbosityLevel >= 1)
464        G4cout << "Theta bin weight " << bweights[3] << " " << rndm << G4endl;
465      return(IPDFThetaBiasH.GetEnergy(rndm));
466    }
467}
468
469G4double G4SPSRandomGenerator::GenRandPhi()
470{
471  if(verbosityLevel >= 1)
472    G4cout << "In GenRandPhi" << G4endl;
473  if(PhiBias == false)
474    {
475      // Phi is not biased
476      G4double rndm = G4UniformRand();
477      return(rndm);
478    }
479  else
480    {
481      // Phi is biased
482      if(IPDFPhiBias == false)
483        {
484          // IPDF has not been created, so create it
485          G4double bins[1024],vals[1024], sum;
486          G4int ii;
487          G4int maxbin = G4int(PhiBiasH.GetVectorLength());
488          bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
489          vals[0] = PhiBiasH(size_t(0));
490          sum = vals[0];
491          for(ii=1;ii<maxbin;ii++)
492            {
493              bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
494              vals[ii] = PhiBiasH(size_t(ii)) + vals[ii-1];
495              sum = sum + PhiBiasH(size_t(ii));
496            }
497
498          for(ii=0;ii<maxbin;ii++)
499            {
500              vals[ii] = vals[ii]/sum;
501              IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
502            }
503          // Make IPDFPhiBias = true
504          IPDFPhiBias = true;
505        }
506      // IPDF has been create so carry on
507      G4double rndm = G4UniformRand();
508      //      size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm);
509      size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
510      G4int biasn1 = 0;
511      G4int biasn2 = numberOfBin/2;
512      G4int biasn3 = numberOfBin - 1;
513      while (biasn1 != biasn3 - 1) {
514      if (rndm > IPDFPhiBiasH(biasn2))
515         biasn1 = biasn2;
516      else
517         biasn3 = biasn2;
518      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
519      }
520      bweights[4] =  IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
521      G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
522      G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
523      G4double NatProb = xaxisu - xaxisl;
524      bweights[4] = NatProb/bweights[4];
525      if(verbosityLevel >= 1)
526        G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
527      return(IPDFPhiBiasH.GetEnergy(rndm));
528    }
529}
530
531G4double G4SPSRandomGenerator::GenRandEnergy()
532{
533  if(verbosityLevel >= 1)
534    G4cout << "In GenRandEnergy" << G4endl;
535  if(EnergyBias == false)
536    {
537      // Energy is not biased
538      G4double rndm = G4UniformRand();
539      return(rndm);
540    }
541  else {
542    // ENERGY is biased
543    if(IPDFEnergyBias == false) {
544      // IPDF has not been created, so create it
545      G4double bins[1024],vals[1024], sum;
546      G4int ii;
547      G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
548      bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
549      vals[0] = EnergyBiasH(size_t(0));
550      sum = vals[0];
551      for(ii=1;ii<maxbin;ii++) {
552        bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
553        vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii-1];
554        sum = sum + EnergyBiasH(size_t(ii));
555      }
556     
557      for(ii=0;ii<maxbin;ii++) {
558        vals[ii] = vals[ii]/sum;
559        IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
560      }
561      // Make IPDFEnergyBias = true
562      IPDFEnergyBias = true;
563    }
564    // IPDF has been create so carry on
565    G4double rndm = G4UniformRand();
566    //  size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
567    size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
568    G4int biasn1 = 0;
569    G4int biasn2 = numberOfBin/2;
570    G4int biasn3 = numberOfBin - 1;
571    while (biasn1 != biasn3 - 1) {
572      if (rndm > IPDFEnergyBiasH(biasn2))
573        biasn1 = biasn2;
574      else
575        biasn3 = biasn2;
576      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
577    }
578    bweights[5] =  IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
579    G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
580    G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
581    G4double NatProb = xaxisu - xaxisl;
582    bweights[5] = NatProb/bweights[5];
583    if(verbosityLevel >= 1)
584      G4cout << "Energy bin weight " << bweights[5] << " " << rndm << G4endl;
585    return(IPDFEnergyBiasH.GetEnergy(rndm));
586  }
587}
588
589
590G4double G4SPSRandomGenerator::GenRandPosTheta()
591{
592  if(verbosityLevel >= 1)
593    {
594      G4cout << "In GenRandPosTheta" << G4endl;
595      G4cout << "Verbosity " << verbosityLevel << G4endl;
596    }
597  if(PosThetaBias == false)
598    {
599      // Theta is not biased
600      G4double rndm = G4UniformRand();
601      return(rndm);
602    }
603  else
604    {
605      // Theta is biased
606      if(IPDFPosThetaBias == false)
607        {
608          // IPDF has not been created, so create it
609          G4double bins[1024],vals[1024], sum;
610          G4int ii;
611          G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
612          bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
613          vals[0] = PosThetaBiasH(size_t(0));
614          sum = vals[0];
615          for(ii=1;ii<maxbin;ii++)
616            {
617              bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
618              vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii-1];
619              sum = sum + PosThetaBiasH(size_t(ii));
620            }
621
622          for(ii=0;ii<maxbin;ii++)
623            {
624              vals[ii] = vals[ii]/sum;
625              IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
626            }
627          // Make IPDFThetaBias = true
628          IPDFPosThetaBias = true;
629        }
630      // IPDF has been create so carry on
631      G4double rndm = G4UniformRand();
632      //      size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
633      size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
634      G4int biasn1 = 0;
635      G4int biasn2 = numberOfBin/2;
636      G4int biasn3 = numberOfBin - 1;
637      while (biasn1 != biasn3 - 1) {
638      if (rndm > IPDFPosThetaBiasH(biasn2))
639         biasn1 = biasn2;
640      else
641         biasn3 = biasn2;
642      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
643      }
644      bweights[6] =  IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
645      G4double xaxisl = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
646      G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
647      G4double NatProb = xaxisu - xaxisl;
648      bweights[6] = NatProb/bweights[6];
649      if(verbosityLevel >= 1)
650        G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm << G4endl;
651      return(IPDFPosThetaBiasH.GetEnergy(rndm));
652    }
653}
654
655G4double G4SPSRandomGenerator::GenRandPosPhi()
656{
657  if(verbosityLevel >= 1)
658    G4cout << "In GenRandPosPhi" << G4endl;
659  if(PosPhiBias == false)
660    {
661      // PosPhi is not biased
662      G4double rndm = G4UniformRand();
663      return(rndm);
664    }
665  else
666    {
667      // PosPhi is biased
668      if(IPDFPosPhiBias == false)
669        {
670          // IPDF has not been created, so create it
671          G4double bins[1024],vals[1024], sum;
672          G4int ii;
673          G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
674          bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
675          vals[0] = PosPhiBiasH(size_t(0));
676          sum = vals[0];
677          for(ii=1;ii<maxbin;ii++)
678            {
679              bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
680              vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii-1];
681              sum = sum + PosPhiBiasH(size_t(ii));
682            }
683
684          for(ii=0;ii<maxbin;ii++)
685            {
686              vals[ii] = vals[ii]/sum;
687              IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
688            }
689          // Make IPDFPosPhiBias = true
690          IPDFPosPhiBias = true;
691        }
692      // IPDF has been create so carry on
693      G4double rndm = G4UniformRand();
694      //      size_t weight_bin_no = IPDFPosPhiBiasH.FindValueBinLocation(rndm);
695      size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
696      G4int biasn1 = 0;
697      G4int biasn2 = numberOfBin/2;
698      G4int biasn3 = numberOfBin - 1;
699      while (biasn1 != biasn3 - 1) {
700      if (rndm > IPDFPosPhiBiasH(biasn2))
701         biasn1 = biasn2;
702      else
703         biasn3 = biasn2;
704      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
705      }
706      bweights[7] =  IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
707      G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
708      G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
709      G4double NatProb = xaxisu - xaxisl;
710      bweights[7] = NatProb/bweights[7];
711      if(verbosityLevel >= 1)
712        G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm << G4endl;
713      return(IPDFPosPhiBiasH.GetEnergy(rndm));
714    }
715}
716
717
718
719
720
Note: See TracBrowser for help on using the repository browser.