WCSim
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
WCSimPrimaryGeneratorAction.cc
Go to the documentation of this file.
1 #include "WCSimPrimaryGeneratorAction.hh"
2 #include "WCSimDetectorConstruction.hh"
3 #include "WCSimPrimaryGeneratorMessenger.hh"
4 #include "G4RunManager.hh"
5 #include "G4Event.hh"
6 #include "G4ParticleGun.hh"
7 #include "G4GeneralParticleSource.hh"
8 #include "G4ParticleTable.hh"
9 #include "G4IonTable.hh"
10 #include "G4ParticleDefinition.hh"
11 #include "G4ThreeVector.hh"
12 #include "globals.hh"
13 #include "Randomize.hh"
14 #include <fstream>
15 #include <vector>
16 #include <string>
17 
18 #include "G4Navigator.hh"
19 #include "G4TransportationManager.hh"
20 
21 #include "G4PhysicalConstants.hh"
22 #include "G4SystemOfUnits.hh"
23 
24 using std::vector;
25 using std::string;
26 using std::fstream;
27 
28 vector<string> tokenize( string separators, string input );
29 
30 inline vector<string> readInLine(fstream& inFile, int lineSize, char* inBuf)
31 {
32  // Read in line break it up into tokens
33  // Any line starting with # is ignored
34  while(true)
35  {
36  if (inFile.getline(inBuf,lineSize))
37  {
38  if(inBuf[0]!='#')
39  return tokenize(" $\r", inBuf);
40  }
41  else
42  {
43  vector<string> nullLine;
44  return nullLine;
45  }
46  }
47 }
48 
49 inline double atof( const string& S ) {return std::atof( S.c_str() );}
50 inline int atoi( const string& S ) {return std::atoi( S.c_str() );}
51 
52 WCSimPrimaryGeneratorAction::WCSimPrimaryGeneratorAction(
53  WCSimDetectorConstruction* myDC)
54  :myDetector(myDC), vectorFileName("")
55 {
56  //T. Akiri: Initialize GPS to allow for the laser use
57  MyGPS = new G4GeneralParticleSource();
58 
59  // Initialize to zero
60  mode[0] = 0;
61  nvtxs = 0;
62  for( Int_t u=0; u<MAX_N_VERTICES; u++){
63  vtxsvol[u] = 0;
64  vtxs[u] = G4ThreeVector(0.,0.,0.);
65  }
66  nuEnergy = 0.;
67  _counterRock=0; // counter for generated in Rock
68  _counterCublic=0; // counter generated
69 
70  //---Set defaults. Do once at beginning of session.
71 
72  G4int n_particle = 1;
73  particleGun = new G4ParticleGun(n_particle);
74  particleGun->SetParticleEnergy(1.0*GeV);
75  particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.0));
76 
77  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
78 // G4IonTable* ionTable = G4IonTable::GetIonTable();
79  G4String particleName;
80  particleGun->
81  SetParticleDefinition(particleTable->FindParticle(particleName="mu+"));
82 
83  particleGun->
84  SetParticlePosition(G4ThreeVector(0.*m,0.*m,0.*m));
85 
86  messenger = new WCSimPrimaryGeneratorMessenger(this);
87  useMulineEvt = true;
88  useGunEvt = false;
89  useLaserEvt = false;
90  useGPSEvt = false;
91  useRadioactiveEvt = false;
92  useRadonEvt = false;
93 
94  // Radioactive and Radon generator variables:
95  radioactive_sources.clear();
96  myRn222Generator = 0;
97  fRnScenario = 0;
98  fRnSymmetry = 1;
99  // Time units for vertices
100  fTimeUnit=CLHEP::nanosecond;
101 }
102 
103 WCSimPrimaryGeneratorAction::~WCSimPrimaryGeneratorAction()
104 {
105  if (IsGeneratingVertexInRock()){
106  G4cout << "Fraction of Rock volume is : " << G4endl;
107  G4cout << " Random number generated in Rock / in Cublic = "
108  << _counterRock << "/" << _counterCublic
109  << " = " << _counterRock/(G4double)_counterCublic << G4endl;
110  }
111  inputFile.close();
112  delete particleGun;
113  delete MyGPS; //T. Akiri: Delete the GPS variable
114  delete messenger;
115 }
116 
117 void WCSimPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
118 {
119  // We will need a particle table
120  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
121  G4IonTable* ionTable = G4IonTable::GetIonTable();
122  // Temporary kludge to turn on/off vector text format
123  G4bool useNuanceTextFormat = true;
124  // Do for every event
125  if (useMulineEvt)
126  {
127  if ( !inputFile.is_open() )
128  {
129  G4cout << "Set a vector file using the command /mygen/vecfile name"
130  << G4endl;
131  exit(-1);
132  }
133  // The original documentation describing the nuance text format can be found here:
134  // http://neutrino.phy.duke.edu/nuance-format/
135  //
136  // Information specific to WCSim can be found in the file Nuance_MC_Format.txt in
137  // the doc directory.
138  // The format must be strictly adhered to for it to be processed correctly.
139  // The lines and their meanings from begin through info are fixed, and then
140  // a variable number of tracks may follow.
141  //
142  if (useNuanceTextFormat)
143  {
144  const int lineSize=100;
145  char inBuf[lineSize];
146  vector<string> token(1);
147 
148  token = readInLine(inputFile, lineSize, inBuf);
149 
150  if (token.size() == 0 || token[0] == "stop" )
151  {
152  G4cout << "End of nuance vector file - run terminated..."<< G4endl;
153  G4RunManager::GetRunManager()-> AbortRun();
154  }
155  else if (token[0] != "begin" )
156  {
157  G4cout << "unexpected line begins with " << token[0] << " we were expecting \" begin \" "<<G4endl;
158  }
159  else // normal parsing begins here
160  {
161  // Read the nuance line
162  // should be nuance <value>
163  // but could be just
164  // nuance
165  // if value is given set mode to equal it.
166 
167  token = readInLine(inputFile, lineSize, inBuf);
168  int iVertex=0;
169  while(token[0]=="nuance" && iVertex < MAX_N_VERTICES)
170  {
171  if(token.size()>1)
172  mode[iVertex] = atoi(token[1]);
173  // Read the Vertex line
174  token = readInLine(inputFile, lineSize, inBuf);
175  vtxs[iVertex] = G4ThreeVector(atof(token[1])*cm,
176  atof(token[2])*cm,
177  atof(token[3])*cm);
178  G4double VertexTime=atof(token[4])*fTimeUnit;
179  vertexTimes[iVertex]=VertexTime;
180  // true : Generate vertex in Rock , false : Generate vertex in WC tank
181  SetGenerateVertexInRock(false);
182  // Next we read the incoming neutrino and target
183  // First, the neutrino line
184  token=readInLine(inputFile, lineSize, inBuf);
185  beampdgs[iVertex] = atoi(token[1]);
186  beamenergies[iVertex] = atof(token[2])*MeV;
187  beamdirs[iVertex] = G4ThreeVector(atof(token[3]),
188  atof(token[4]),
189  atof(token[5]));
190  SetBeamEnergy(beamenergies[iVertex]);
191  SetBeamDir(beamdirs[iVertex]);
192  // Now read the target line
193  token=readInLine(inputFile, lineSize, inBuf);
194  targetpdgs[iVertex] = atoi(token[1]);
195  targetenergies[iVertex] = atof(token[2])*MeV;
196  targetdirs[iVertex] = G4ThreeVector(atof(token[3]),
197  atof(token[4]),
198  atof(token[5]));
199  // Read the info line, basically a dummy
200  token=readInLine(inputFile, lineSize, inBuf);
201  vecRecNumber = atoi(token[2]);
202  // Now read the outgoing particles
203  // These we will simulate.
204  while ( token=readInLine(inputFile, lineSize, inBuf),
205  token[0] == "track" )
206  {
207  // We are only interested in the particles
208  // that leave the nucleus, tagged by "0"
209  if ( token[6] == "0")
210  {
211  G4int pdgid = atoi(token[1]);
212  G4double tempEnergy = atof(token[2])*MeV;
213  G4ThreeVector dir = G4ThreeVector(atof(token[3]),
214  atof(token[4]),
215  atof(token[5]));
216  //must handle the case of an ion seperatly from other particles
217  //check PDG code if we have an ion.
218  //PDG code format for ions ±10LZZZAAAI
219  char strPDG[11];
220  char strA[10]={0};
221  char strZ[10]={0};
222  long int A=0,Z=0;
223  if(abs(pdgid) >= 1000000000)
224  {
225  //ion
226  sprintf(strPDG,"%i",abs(pdgid));
227  strncpy(strZ, &strPDG[3], 3);
228  strncpy(strA, &strPDG[6], 3);
229  strA[3]='\0';
230  strZ[3]='\0';
231  A=atoi(strA);
232  Z=atoi(strZ);
233  G4ParticleDefinition* ion;
234  ion = ionTable->GetIon(Z, A, 0.);
235  particleGun->SetParticleDefinition(ion);
236  particleGun->SetParticleCharge(0);
237  }
238  else {
239  //not ion
240  particleGun->
241  SetParticleDefinition(particleTable->
242  FindParticle(pdgid));
243  }
244  G4double mass =
245  particleGun->GetParticleDefinition()->GetPDGMass();
246 
247  G4double ekin = tempEnergy - mass;
248 
249  particleGun->SetParticleEnergy(ekin);
250  particleGun->SetParticlePosition(vtxs[iVertex]);
251  particleGun->SetParticleMomentumDirection(dir);
252  particleGun->SetParticleTime(VertexTime);
253  particleGun->GeneratePrimaryVertex(anEvent);
254 
255  }
256  }
257  iVertex++;
258  if(iVertex > MAX_N_VERTICES)
259  G4cout<<" CAN NOT DEAL WITH MORE THAN "<<MAX_N_VERTICES<<" VERTICES - TRUNCATING EVENT HERE "<<G4endl;
260  }
261  nvtxs=iVertex;
262  SetNvtxs(nvtxs);
263 
264  }
265  }
266  else
267  { // old muline format
268  inputFile >> nuEnergy >> energy >> xPos >> yPos >> zPos
269  >> xDir >> yDir >> zDir;
270 
271  G4double random_z = ((myDetector->GetWaterTubePosition())
272  - .5*(myDetector->GetWaterTubeLength())
273  + 1.*m + 15.0*m*G4UniformRand())/m;
274  zPos = random_z;
275  G4ThreeVector vtx = G4ThreeVector(xPos, yPos, random_z);
276  G4ThreeVector dir = G4ThreeVector(xDir,yDir,zDir);
277 
278  particleGun->SetParticleEnergy(energy*MeV);
279  particleGun->SetParticlePosition(vtx);
280  particleGun->SetParticleMomentumDirection(dir);
281  particleGun->GeneratePrimaryVertex(anEvent);
282  }
283  }
284 
285  else if (useGunEvt)
286  { // manual gun operation
287  particleGun->GeneratePrimaryVertex(anEvent);
288 
289  //To prevent occasional seg fault from an un assigned targetpdg
290  targetpdgs[0] = 2212; //ie. proton
291 
292  G4ThreeVector P =anEvent->GetPrimaryVertex()->GetPrimary()->GetMomentum();
293  G4ThreeVector vtx=anEvent->GetPrimaryVertex()->GetPosition();
294  G4double mass =anEvent->GetPrimaryVertex()->GetPrimary()->GetMass();
295  G4int pdg =anEvent->GetPrimaryVertex()->GetPrimary()->GetPDGcode();
296 
297  char strPDG[11];
298  char strA[10]={0};
299  char strZ[10]={0};
300 
301 
302  long int A=0,Z=0;
303  // A=strotl(strPDG,&str);
304  if(abs(pdg) >= 1000000000)
305  {
306  //ion
307  sprintf(strPDG,"%i",abs(pdg));
308  strncpy(strZ, &strPDG[3], 3);
309  strncpy(strA, &strPDG[6], 3);
310  strA[3]='\0';
311  strZ[3]='\0';
312  A=atoi(strA);
313  Z=atoi(strZ);
314 
315  G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon(Z, A, 0);
316  ion->SetPDGStable(false);
317  ion->SetPDGLifeTime(0.);
318 
319  G4ParticleDefinition* ion2 = G4IonTable::GetIonTable()->GetIon(Z, A, 0);
320  std::cout<<"ion2 "<<ion2->GetPDGLifeTime()<<"\n";
321  }
322 
323  G4ThreeVector dir = P.unit();
324  G4double E = std::sqrt((P.dot(P))+(mass*mass));
325 
326  SetVtx(vtx);
327  SetBeamEnergy(E);
328  SetBeamDir(dir);
329  SetBeamPDG(pdg);
330  }
331  else if (useLaserEvt)
332  {
333  targetpdgs[0] = 2212; //ie. proton
334  //T. Akiri: Create the GPS LASER event
335  MyGPS->GeneratePrimaryVertex(anEvent);
336 
337  G4ThreeVector P =anEvent->GetPrimaryVertex()->GetPrimary()->GetMomentum();
338  G4ThreeVector vtx =anEvent->GetPrimaryVertex()->GetPosition();
339  G4int pdg =anEvent->GetPrimaryVertex()->GetPrimary()->GetPDGcode();
340 
341  // G4ThreeVector dir = P.unit();
342  G4double E = std::sqrt((P.dot(P)));
343 
344  //SetVtx(vtx);
345  SetBeamEnergy(E);
346  //SetBeamDir(dir);
347  SetBeamPDG(pdg);
348  }
349  else if (useGPSEvt)
350  {
351  MyGPS->GeneratePrimaryVertex(anEvent);
352 
353  G4ThreeVector P =anEvent->GetPrimaryVertex()->GetPrimary()->GetMomentum();
354  G4ThreeVector vtx =anEvent->GetPrimaryVertex()->GetPosition();
355  G4double mass =anEvent->GetPrimaryVertex()->GetPrimary()->GetMass();
356  G4int pdg =anEvent->GetPrimaryVertex()->GetPrimary()->GetPDGcode();
357 
358  G4ThreeVector dir = P.unit();
359  G4double E = std::sqrt((P.dot(P))+(mass*mass));
360 
361  SetVtx(vtx);
362  SetBeamEnergy(E);
363  SetBeamDir(dir);
364  SetBeamPDG(pdg);
365  }
366  else if (useRadioactiveEvt)
367  {
368 
369  // initialize GPS properties
370  MyGPS->ClearAll();
371 
372  MyGPS->SetMultipleVertex(true);
373 
374  std::vector<WCSimPmtInfo*> *pmts=NULL;
375 
376  std::vector<struct radioactive_source>::iterator it;
377 
378  for ( it = radioactive_sources.begin(); it != radioactive_sources.end(); it++ ){
379  G4String IsotopeName = it->IsotopeName;
380  G4String IsotopeLocation = it->IsotopeLocation;
381  G4double IsotopeActivity = it->IsotopeActivity;
382 
383  double average= IsotopeActivity * GetRadioactiveTimeWindow();
384  if (IsotopeLocation.compareTo("PMT") == 0){
385  pmts = myDetector->Get_Pmts();
386  average *= pmts->size();
387  }
388 
389  // random poisson number of vertices based on average
390  int n_vertices = CLHEP::RandPoisson::shoot(average);
391 
392  // n_vertices = 1; // qqq
393 
394  for(int u=0; u<n_vertices; u++){
395 
396  MyGPS->AddaSource(1.);
397 
398  MyGPS->SetCurrentSourceto(MyGPS->GetNumberofSource() - 1);
399 
400  if (IsotopeName.compareTo("Tl208") == 0){
401  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 81, 208, 0));
402  }else if (IsotopeName.compareTo("Bi214") == 0){
403  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 83, 214, 0));
404  }else if (IsotopeName.compareTo("K40") == 0){
405  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 19, 40, 0));
406  }else if (IsotopeName.compareTo("Rn220") == 0){
407  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 86, 220, 0));
408  }else if (IsotopeName.compareTo("Po216") == 0){
409  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 84, 216, 0));
410  }else if (IsotopeName.compareTo("Pb212") == 0){
411  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 82, 212, 0));
412  }else if (IsotopeName.compareTo("Bi212") == 0){
413  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 83, 212, 0));
414  }else if (IsotopeName.compareTo("Po212") == 0){
415  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 84, 212, 0));
416  }else if (IsotopeName.compareTo("Rn222") == 0){
417  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 86, 222, 0));
418  }else if (IsotopeName.compareTo("Po218") == 0){
419  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 84, 218, 0));
420  }else if (IsotopeName.compareTo("At218") == 0){
421  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 85, 218, 0));
422  }else if (IsotopeName.compareTo("Pb214") == 0){
423  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 82, 214, 0));
424  }else if (IsotopeName.compareTo("Po214") == 0){
425  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 84, 214, 0));
426  }else if (IsotopeName.compareTo("Tl210") == 0){
427  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 81, 210, 0));
428  }else if (IsotopeName.compareTo("Pb210") == 0){
429  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 82, 210, 0));
430  }else if (IsotopeName.compareTo("Bi210") == 0){
431  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 83, 210, 0));
432  }else if (IsotopeName.compareTo("Po210") == 0){
433  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 84, 210, 0));
434  }else if (IsotopeName.compareTo("Hg206") == 0){
435  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 80, 206, 0));
436  }else if (IsotopeName.compareTo("Tl206") == 0){
437  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 81, 206, 0));
438  }else if (IsotopeName.compareTo("Rn219") == 0){
439  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 86, 219, 0));
440  }else if (IsotopeName.compareTo("Po215") == 0){
441  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 84, 215, 0));
442  }else if (IsotopeName.compareTo("At215") == 0){
443  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 85, 215, 0));
444  }else if (IsotopeName.compareTo("Pb211") == 0){
445  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 82, 211, 0));
446  }else if (IsotopeName.compareTo("Bi211") == 0){
447  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 83, 211, 0));
448  }else if (IsotopeName.compareTo("Po211") == 0){
449  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 84, 211, 0));
450  }else if (IsotopeName.compareTo("Tl207") == 0){
451  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 81, 207, 0));
452  }else if (IsotopeName.compareTo("Th232") == 0){
453  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 90, 232, 0));
454  }else if (IsotopeName.compareTo("Ra228") == 0){
455  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 88, 228, 0));
456  }else if (IsotopeName.compareTo("Ac228") == 0){
457  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 89, 228, 0));
458  }else if (IsotopeName.compareTo("Th228") == 0){
459  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 90, 228, 0));
460  }else if (IsotopeName.compareTo("Ra224") == 0){
461  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 88, 224, 0));
462  }else if (IsotopeName.compareTo("U238") == 0){
463  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 92, 238, 0));
464  }else if (IsotopeName.compareTo("Th234") == 0){
465  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 90, 234, 0));
466  }else if (IsotopeName.compareTo("Pa234") == 0){
467  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 91, 234, 0));
468  }else if (IsotopeName.compareTo("U234") == 0){
469  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 92, 234, 0));
470  }else if (IsotopeName.compareTo("Th230") == 0){
471  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 90, 230, 0));
472  }else if (IsotopeName.compareTo("Ra226") == 0){
473  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 88, 226, 0));
474  }else if (IsotopeName.compareTo("U235") == 0){
475  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 92, 235, 0));
476  }else if (IsotopeName.compareTo("Th231") == 0){
477  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 90, 231, 0));
478  }else if (IsotopeName.compareTo("Pa231") == 0){
479  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 91, 231, 0));
480  }else if (IsotopeName.compareTo("Ac227") == 0){
481  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 89, 227, 0));
482  }else if (IsotopeName.compareTo("Th227") == 0){
483  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 90, 227, 0));
484  }else if (IsotopeName.compareTo("Fr223") == 0){
485  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 87, 223, 0));
486  }else if (IsotopeName.compareTo("Ra223") == 0){
487  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 88, 223, 0));
488  }else if (IsotopeName.compareTo("At219") == 0){
489  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 85, 219, 0));
490  }else if (IsotopeName.compareTo("Bi215") == 0){
491  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 83, 215, 0));
492  }
493 
494  if (IsotopeLocation.compareTo("water") == 0){
495  MyGPS->GetCurrentSource()->GetEneDist()->SetEnergyDisType("Mono");
496  MyGPS->GetCurrentSource()->GetEneDist()->SetMonoEnergy(0.);
497  MyGPS->GetCurrentSource()->GetPosDist()->SetPosDisType("Point");
498  MyGPS->GetCurrentSource()->GetPosDist()->SetCentreCoords(G4ThreeVector(0, 0, 0));
499  MyGPS->GetCurrentSource()->GetPosDist()->SetPosDisType("Volume");
500  MyGPS->GetCurrentSource()->GetPosDist()->SetPosDisShape("Cylinder");
501  G4String WCIDCollectionName = myDetector->GetIDCollectionName();
502  WCSimPMTObject *PMT = myDetector->GetPMTPointer(WCIDCollectionName);
503  MyGPS->GetCurrentSource()->GetPosDist()->SetRadius(myDetector->GetGeo_Dm(3)*CLHEP::cm - 2.*PMT->GetRadius());
504  MyGPS->GetCurrentSource()->GetPosDist()->SetHalfZ(myDetector->GetGeo_Dm(2)*CLHEP::cm/2. - 2.*PMT->GetRadius());
505  MyGPS->GetCurrentSource()->GetPosDist()->SetPosRot1(G4ThreeVector(1, 0, 0));
506  MyGPS->GetCurrentSource()->GetPosDist()->SetPosRot2(G4ThreeVector(0, 1, 0));
507 
508  }
509  else if (IsotopeLocation.compareTo("PMT") == 0){
510  int npmts = pmts->size();
511  int random_pmt_id = CLHEP::RandFlat::shootInt(1,npmts);
512  WCSimPmtInfo* pmtinfo = (WCSimPmtInfo*)pmts->at( random_pmt_id - 1 );
513  G4ThreeVector random_pmt_center(pmtinfo->Get_transx()*CLHEP::cm, pmtinfo->Get_transy()*CLHEP::cm, pmtinfo->Get_transz()*CLHEP::cm);
514  double random_cos_theta = CLHEP::RandFlat::shoot(0., 1.);
515  double random_sin_theta = sqrt(1. - pow(random_cos_theta,2));
516  random_sin_theta *= (CLHEP::RandFlat::shootBit() == 0 ? -1 : 1);
517  double random_phi = CLHEP::RandFlat::shoot(0., 2.*CLHEP::pi*CLHEP::rad);
518  G4String WCIDCollectionName = myDetector->GetIDCollectionName();
519  WCSimPMTObject *PMT = myDetector->GetPMTPointer(WCIDCollectionName);
520  double PMT_radius = PMT->GetRadius();
521  double glassThickness = PMT->GetPMTGlassThickness();
522  double expose = PMT->GetExposeHeight();
523  double sphereRadius = (expose*expose+ PMT_radius*PMT_radius)/(2*expose);
524  double Rmin = sphereRadius-glassThickness;
525  double Rmax = sphereRadius;
526  double random_R = CLHEP::RandFlat::shoot(Rmin, Rmax);
527  G4ThreeVector orientation(pmtinfo->Get_orienx(), pmtinfo->Get_orieny(), pmtinfo->Get_orienz());
528  G4ThreeVector axis_1 = orientation.orthogonal();
529  G4ThreeVector axis_2 = orientation.cross(axis_1);
530  G4ThreeVector position = random_pmt_center + random_R*(orientation*random_cos_theta + axis_1*random_sin_theta*cos(random_phi) + axis_2*random_sin_theta*sin(random_phi));
531 
532  //G4cout << " random id " << random_pmt_id << " of " << npmts << " costheta " << random_cos_theta << " sintheta " << random_sin_theta << " phi " << random_phi << " WCIDCollectionName " << WCIDCollectionName << " PMT_radius " << PMT_radius << " expose " << expose << " sphereRadius " << sphereRadius << " Rmin " << Rmin << " Rmax " << Rmax << " random_R " << random_R << " orientation (" << orientation.x() << ", " << orientation.y() << ", " << orientation.z() << ") center (" << random_pmt_center.x() << ", " << random_pmt_center.y() << ", " << random_pmt_center.z() << ") position (" << position.x() << ", " << position.y() << ", " << position.z() << ") " << G4endl;
533 
534  MyGPS->GetCurrentSource()->GetEneDist()->SetEnergyDisType("Mono");
535  MyGPS->GetCurrentSource()->GetEneDist()->SetMonoEnergy(0.);
536  MyGPS->GetCurrentSource()->GetPosDist()->SetPosDisType("Point");
537  MyGPS->GetCurrentSource()->GetPosDist()->SetCentreCoords(position);
538  }
539 
540  }
541 
542 // G4cout << " is " << IsotopeName << " of " << radioactive_sources.size() << " loc " << IsotopeLocation << " a " << IsotopeActivity << " nv " << n_vertices << G4endl;
543 
544  }
545 
546  G4int number_of_sources = MyGPS->GetNumberofSource();
547 
548  // this will generate several primary vertices
549  MyGPS->GeneratePrimaryVertex(anEvent);
550 
551  SetNvtxs(number_of_sources);
552  for( G4int u=0; u<number_of_sources; u++){
553  targetpdgs[u] = 2212; //ie. proton
554 
555  G4ThreeVector P =anEvent->GetPrimaryVertex(u)->GetPrimary()->GetMomentum();
556  G4ThreeVector vtx =anEvent->GetPrimaryVertex(u)->GetPosition();
557  G4int pdg =anEvent->GetPrimaryVertex(u)->GetPrimary()->GetPDGcode();
558 
559  // G4ThreeVector dir = P.unit();
560  G4double E = std::sqrt((P.dot(P)));
561 
562 // G4cout << " vertex " << u << " of " << number_of_sources << " (" << vtx.x() << ", " << vtx.y() << ", " << vtx.z() << ")" << G4endl;
563 
564  SetVtxs(u,vtx);
565  SetBeamEnergy(E,u);
566  // SetBeamDir(dir);
567  SetBeamPDG(pdg,u);
568  }
569 
570  }
571  else if (useRadonEvt)
572  { //G. Pronost: Add Radon (adaptation of Radioactive event)
573 
574  // Currently only one generator is possible
575  // In order to have several, we need to find a solution for the fitting graphes (which are static currently)
576  // Idea: array of fitting graphes? (each new generators having a specific ID)
577  if ( !myRn222Generator ) {
578  myRn222Generator = new WCSimGenerator_Radioactivity(myDetector);
579  myRn222Generator->Configuration(fRnScenario);
580  }
581 
582  //G4cout << " Generate radon events " << G4endl;
583  // initialize GPS properties
584  MyGPS->ClearAll();
585 
586  MyGPS->SetMultipleVertex(true);
587 
588 
589  std::vector<struct radioactive_source>::iterator it;
590 
591  G4String IsotopeName = "Rn222";
592  G4double IsotopeActivity = myRn222Generator->GetMeanActivity() * 1e-3; // mBq to Bq
593  G4double iEventAvg = IsotopeActivity * GetRadioactiveTimeWindow();
594 
595  //G4cout << " Average " << iEventAvg << G4endl;
596  // random poisson number of vertices based on average
597  int n_vertices = CLHEP::RandPoisson::shoot(iEventAvg);
598 
599  if ( n_vertices < 1 ) {
600  n_vertices = 1;
601  }
602 
603  for(int u=0; u<n_vertices; u++){
604 
605  MyGPS->AddaSource(1.);
606  MyGPS->SetCurrentSourceto(MyGPS->GetNumberofSource() - 1);
607 
608  // Bi214 (source of electron in Rn222 decay chain, assumed to be in equilibrium)
609  MyGPS->SetParticleDefinition(G4IonTable::GetIonTable()->GetIon( 83, 214, 0));
610 
611  // Get position (first position take few seconds to be produced, there after there is no trouble)
612  //G4cout << "GetRandomVertex" << G4endl;
613  G4ThreeVector position = myRn222Generator->GetRandomVertex(fRnSymmetry);
614  //G4cout << "Done: " << position << G4endl;
615  // energy
616  MyGPS->GetCurrentSource()->GetEneDist()->SetEnergyDisType("Mono");
617  MyGPS->GetCurrentSource()->GetEneDist()->SetMonoEnergy(0.);
618 
619  // position
620  MyGPS->GetCurrentSource()->GetPosDist()->SetPosDisType("Point");
621  MyGPS->GetCurrentSource()->GetPosDist()->SetCentreCoords(position);
622 
623  //G4cout << u << " is " << IsotopeName << " loc " << position << G4endl;
624 
625  }
626  G4int number_of_sources = MyGPS->GetNumberofSource();
627 
628  // this will generate several primary vertices
629  MyGPS->GeneratePrimaryVertex(anEvent);
630 
631  SetNvtxs(number_of_sources);
632  for( G4int u=0; u<number_of_sources; u++){
633  targetpdgs[u] = 2212; //ie. proton
634 
635  G4ThreeVector P =anEvent->GetPrimaryVertex(u)->GetPrimary()->GetMomentum();
636  G4ThreeVector vtx =anEvent->GetPrimaryVertex(u)->GetPosition();
637  G4int pdg =anEvent->GetPrimaryVertex(u)->GetPrimary()->GetPDGcode();
638 
639  // G4ThreeVector dir = P.unit();
640  G4double E = std::sqrt((P.dot(P)));
641 
642  //G4cout << " vertex " << u << " of " << number_of_sources << " (" << vtx.x() << ", " << vtx.y() << ", " << vtx.z() << ") with pdg: " << pdg << G4endl;
643 
644  SetVtxs(u,vtx);
645  SetBeamEnergy(E,u);
646  // SetBeamDir(dir);
647  SetBeamPDG(pdg,u);
648  }
649 
650  }
651 }
652 
653 void WCSimPrimaryGeneratorAction::SaveOptionsToOutput(WCSimRootOptions * wcopt)
654 {
655  if(useMulineEvt)
656  wcopt->SetVectorFileName(vectorFileName);
657  else
658  wcopt->SetVectorFileName("");
659  wcopt->SetGeneratorType(GetGeneratorTypeString());
660 }
661 
662 G4String WCSimPrimaryGeneratorAction::GetGeneratorTypeString()
663 {
664  if(useMulineEvt)
665  return "muline";
666  else if(useGunEvt)
667  return "gun";
668  else if(useGPSEvt)
669  return "gps";
670  else if(useLaserEvt)
671  return "laser";
672  return "";
673 }
674 
675 // Returns a vector with the tokens
676 vector<string> tokenize( string separators, string input )
677 {
678  std::size_t startToken = 0, endToken; // Pointers to the token pos
679  vector<string> tokens; // Vector to keep the tokens
680 
681  if( separators.size() > 0 && input.size() > 0 )
682  {
683 
684  while( startToken < input.size() )
685  {
686  // Find the start of token
687  startToken = input.find_first_not_of( separators, startToken );
688 
689  // If found...
690  if( startToken != input.npos )
691  {
692  // Find end of token
693  endToken = input.find_first_of( separators, startToken );
694  if( endToken == input.npos )
695  // If there was no end of token, assign it to the end of string
696  endToken = input.size();
697 
698  // Extract token
699  tokens.push_back( input.substr( startToken, endToken - startToken ) );
700 
701  // Update startToken
702  startToken = endToken;
703  }
704  }
705  }
706 
707  return tokens;
708 }
709 
tokenize
vector< string > tokenize(string separators, string input)
Definition: WCSimPrimaryGeneratorAction.cc:676
WCSimPrimaryGeneratorAction::myRn222Generator
WCSimGenerator_Radioactivity * myRn222Generator
Definition: WCSimPrimaryGeneratorAction.hh:100
WCSimPmtInfo::Get_orienx
Double_t Get_orienx()
Definition: WCSimPmtInfo.hh:36
WCSimPmtInfo::Get_transy
Double_t Get_transy()
Definition: WCSimPmtInfo.hh:34
WCSimPrimaryGeneratorAction::messenger
WCSimPrimaryGeneratorMessenger * messenger
Definition: WCSimPrimaryGeneratorAction.hh:81
WCSimPmtInfo::Get_orieny
Double_t Get_orieny()
Definition: WCSimPmtInfo.hh:37
WCSimGenerator_Radioactivity::GetMeanActivity
G4double GetMeanActivity()
Definition: WCSimGenerator_Radioactivity.hh:27
WCSimPrimaryGeneratorAction::SetBeamDir
void SetBeamDir(G4ThreeVector i, G4int n=0)
Definition: WCSimPrimaryGeneratorAction.hh:42
WCSimGenerator_Radioactivity::GetRandomVertex
G4ThreeVector GetRandomVertex(G4int tSymNumber)
Definition: WCSimGenerator_Radioactivity.cc:634
WCSimDetectorConstruction::GetPMTPointer
WCSimPMTObject * GetPMTPointer(G4String CollectionName)
Definition: WCSimDetectorConstruction.hh:112
WCSimPrimaryGeneratorMessenger.hh
WCSimPrimaryGeneratorAction::myDetector
WCSimDetectorConstruction * myDetector
Definition: WCSimPrimaryGeneratorAction.hh:78
WCSimPrimaryGeneratorAction::radioactive_sources
std::vector< struct radioactive_source > radioactive_sources
Definition: WCSimPrimaryGeneratorAction.hh:96
WCSimRootOptions::SetGeneratorType
void SetGeneratorType(string iGeneratorType)
Definition: WCSimRootOptions.hh:137
WCSimPrimaryGeneratorAction::SetBeamEnergy
void SetBeamEnergy(G4double i, G4int n=0)
Definition: WCSimPrimaryGeneratorAction.hh:41
WCSimPmtInfo::Get_transz
Double_t Get_transz()
Definition: WCSimPmtInfo.hh:35
WCSimPrimaryGeneratorAction::vtxs
G4ThreeVector vtxs[MAX_N_VERTICES]
Definition: WCSimPrimaryGeneratorAction.hh:111
WCSimPmtInfo::Get_transx
Double_t Get_transx()
Definition: WCSimPmtInfo.hh:33
WCSimPrimaryGeneratorAction::SetVtxs
void SetVtxs(G4int i, G4ThreeVector v)
Definition: WCSimPrimaryGeneratorAction.hh:45
WCSimGenerator_Radioactivity::Configuration
void Configuration(G4int iScenario, G4double dLifeTime=0)
Definition: WCSimGenerator_Radioactivity.cc:264
WCSimPrimaryGeneratorAction::beamdirs
G4ThreeVector beamdirs[MAX_N_PRIMARIES]
Definition: WCSimPrimaryGeneratorAction.hh:115
WCSimGenerator_Radioactivity
Definition: WCSimGenerator_Radioactivity.hh:17
WCSimPrimaryGeneratorAction::targetenergies
G4double targetenergies[MAX_N_PRIMARIES]
Definition: WCSimPrimaryGeneratorAction.hh:116
readInLine
vector< string > readInLine(fstream &inFile, int lineSize, char *inBuf)
Definition: WCSimPrimaryGeneratorAction.cc:30
WCSimRootOptions::SetVectorFileName
void SetVectorFileName(string iVectorFileName)
Definition: WCSimRootOptions.hh:136
WCSimDetectorConstruction.hh
WCSimPrimaryGeneratorAction::SaveOptionsToOutput
void SaveOptionsToOutput(WCSimRootOptions *wcopt)
Definition: WCSimPrimaryGeneratorAction.cc:653
WCSimPrimaryGeneratorAction::WCSimPrimaryGeneratorAction
WCSimPrimaryGeneratorAction(WCSimDetectorConstruction *)
Definition: WCSimPrimaryGeneratorAction.cc:52
WCSimPMTObject::GetExposeHeight
virtual G4double GetExposeHeight()=0
atof
double atof(const string &S)
Definition: WCSimPrimaryGeneratorAction.cc:49
WCSimDetectorConstruction::Get_Pmts
std::vector< WCSimPmtInfo * > * Get_Pmts()
Definition: WCSimDetectorConstruction.hh:153
WCSimPrimaryGeneratorAction::beamenergies
G4double beamenergies[MAX_N_PRIMARIES]
Definition: WCSimPrimaryGeneratorAction.hh:116
MAX_N_VERTICES
static const int MAX_N_VERTICES
Definition: jhfNtuple.h:3
atoi
int atoi(const string &S)
Definition: WCSimPrimaryGeneratorAction.cc:50
WCSimPMTObject::GetRadius
virtual G4double GetRadius()=0
WCSimPMTObject::GetPMTGlassThickness
virtual G4double GetPMTGlassThickness()=0
WCSimPmtInfo::Get_orienz
Double_t Get_orienz()
Definition: WCSimPmtInfo.hh:38
WCSimPrimaryGeneratorAction.hh
WCSimPrimaryGeneratorAction::targetdirs
G4ThreeVector targetdirs[MAX_N_PRIMARIES]
Definition: WCSimPrimaryGeneratorAction.hh:115
Generated on Tue Nov 3 2020 13:31:01 for WCSim by   doxygen 1.8.5