ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
B02ImportanceDetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file B02ImportanceDetectorConstruction.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 //
28 //
29 //
30 //
31 
32 #include "globals.hh"
33 #include <sstream>
34 
35 #include "B02ImportanceDetectorConstruction.hh"
36 
37 #include "G4Material.hh"
38 #include "G4Tubs.hh"
39 #include "G4LogicalVolume.hh"
40 #include "G4ThreeVector.hh"
41 #include "G4PVPlacement.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 
45 // For Primitive Scorers
46 #include "G4SDManager.hh"
47 #include "G4MultiFunctionalDetector.hh"
48 #include "G4SDParticleFilter.hh"
49 #include "G4PSNofCollision.hh"
50 #include "G4PSPopulation.hh"
51 #include "G4PSTrackCounter.hh"
52 #include "G4PSTrackLength.hh"
53 
54 // for importance biasing
55 #include "G4IStore.hh"
56 
57 // for weight window technique
58 #include "G4WeightWindowStore.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 B02ImportanceDetectorConstruction::
63 B02ImportanceDetectorConstruction(G4String worldName)
64 :G4VUserParallelWorld(worldName),fLogicalVolumeVector()
65 {
66  // Construct();
67 }
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
71 B02ImportanceDetectorConstruction::~B02ImportanceDetectorConstruction()
72 {
73  fLogicalVolumeVector.clear();
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77 
78 void B02ImportanceDetectorConstruction::Construct()
79 {
80  G4cout << " constructing parallel world " << G4endl;
81 
82  G4Material* dummyMat = 0;
83 
84  //GetWorld methods create a clone of the mass world to the parallel world (!)
85  // via the transportation manager
86  fGhostWorld = GetWorld();
87  G4cout << " B02ImportanceDetectorConstruction:: ghostWorldName = "
88  << fGhostWorld->GetName() << G4endl;
89  G4LogicalVolume* worldLogical = fGhostWorld->GetLogicalVolume();
90  fLogicalVolumeVector.push_back(worldLogical);
91 
92  // fPVolumeStore.AddPVolume(G4GeometryCell(*pWorldVolume, 0));
93  fPVolumeStore.AddPVolume(G4GeometryCell(*fGhostWorld, 0));
94 
95  // creating 18 slobs of 10 cm thicknes
96 
97  G4double innerRadiusShield = 0*cm;
98  G4double outerRadiusShield = 100*cm;
99  G4double heightShield = 5*cm;
100  G4double startAngleShield = 0*deg;
101  G4double spanningAngleShield = 360*deg;
102 
103  G4Tubs *aShield = new G4Tubs("aShield",
104  innerRadiusShield,
105  outerRadiusShield,
106  heightShield,
107  startAngleShield,
108  spanningAngleShield);
109 
110  // logical parallel cells
111 
112  G4LogicalVolume *aShield_log_imp =
113  new G4LogicalVolume(aShield, dummyMat, "aShield_log_imp");
114  fLogicalVolumeVector.push_back(aShield_log_imp);
115 
116  // physical parallel cells
117  G4String name = "none";
118  G4int i = 1;
119  G4double startz = -85*cm;
120  // for (i=1; i<=18; ++i) {
121  for (i=1; i<=18; i++) {
122 
123  name = GetCellName(i);
124 
125  G4double pos_x = 0*cm;
126  G4double pos_y = 0*cm;
127  G4double pos_z = startz + (i-1) * (2*heightShield);
128  G4VPhysicalVolume *pvol =
129  new G4PVPlacement(0,
130  G4ThreeVector(pos_x, pos_y, pos_z),
131  aShield_log_imp,
132  name,
133  worldLogical,
134  false,
135  i);
136  // 0);
137  G4GeometryCell cell(*pvol, i);
138  // G4GeometryCell cell(*pvol, 0);
139  fPVolumeStore.AddPVolume(cell);
140  }
141 
142  // filling the rest of the world volumr behind the concrete with
143  // another slob which should get the same importance value as the
144  // last slob
145  innerRadiusShield = 0*cm;
146  // outerRadiusShield = 110*cm; exceeds world volume!!!!
147  outerRadiusShield = 100*cm;
148  // heightShield = 10*cm;
149  heightShield = 5*cm;
150  startAngleShield = 0*deg;
151  spanningAngleShield = 360*deg;
152 
153  G4Tubs *aRest = new G4Tubs("Rest",
154  innerRadiusShield,
155  outerRadiusShield,
156  heightShield,
157  startAngleShield,
158  spanningAngleShield);
159 
160  G4LogicalVolume *aRest_log =
161  new G4LogicalVolume(aRest, dummyMat, "aRest_log");
162 
163  fLogicalVolumeVector.push_back(aRest_log);
164 
165  name = GetCellName(19);
166 
167  G4double pos_x = 0*cm;
168  G4double pos_y = 0*cm;
169  // G4double pos_z = 100*cm;
170  G4double pos_z = 95*cm;
171  G4VPhysicalVolume *pvol =
172  new G4PVPlacement(0,
173  G4ThreeVector(pos_x, pos_y, pos_z),
174  aRest_log,
175  name,
176  worldLogical,
177  false,
178  19);
179  // 0);
180  G4GeometryCell cell(*pvol, 19);
181  // G4GeometryCell cell(*pvol, 0);
182  fPVolumeStore.AddPVolume(cell);
183 
184  SetSensitive();
185 
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
189 
190 const G4VPhysicalVolume &B02ImportanceDetectorConstruction::
191 GetPhysicalVolumeByName(const G4String& name) const {
192  return *fPVolumeStore.GetPVolume(name);
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196 
197 G4String B02ImportanceDetectorConstruction::ListPhysNamesAsG4String(){
198  G4String names(fPVolumeStore.GetPNames());
199  return names;
200 }
201 
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
203 
204 G4String B02ImportanceDetectorConstruction::GetCellName(G4int i) {
205  std::ostringstream os;
206  os << "cell_";
207  if (i<10) {
208  os << "0";
209  }
210  os << i;
211  G4String name = os.str();
212  return name;
213 }
214 
215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
216 
217 G4GeometryCell B02ImportanceDetectorConstruction::GetGeometryCell(G4int i){
218  G4String name(GetCellName(i));
219  const G4VPhysicalVolume *p=0;
220  p = fPVolumeStore.GetPVolume(name);
221  if (p) {
222  return G4GeometryCell(*p,0);
223  }
224  else {
225  G4cout << "B02ImportanceDetectorConstruction::GetGeometryCell: " << G4endl
226  << " couldn't get G4GeometryCell" << G4endl;
227  return G4GeometryCell(*fGhostWorld,-2);
228  }
229 }
230 
231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
232 
233 G4VPhysicalVolume &B02ImportanceDetectorConstruction::
234 GetWorldVolumeAddress() const{
235  return *fGhostWorld;
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239 
240 G4VPhysicalVolume *B02ImportanceDetectorConstruction::GetWorldVolume() {
241  return fGhostWorld;
242 }
243 
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
245 
246 void B02ImportanceDetectorConstruction::SetSensitive(){
247 
248  // -------------------------------------------------
249  // The collection names of defined Primitives are
250  // 0 ConcreteSD/Collisions
251  // 1 ConcreteSD/CollWeight
252  // 2 ConcreteSD/Population
253  // 3 ConcreteSD/TrackEnter
254  // 4 ConcreteSD/SL
255  // 5 ConcreteSD/SLW
256  // 6 ConcreteSD/SLWE
257  // 7 ConcreteSD/SLW_V
258  // 8 ConcreteSD/SLWE_V
259  // -------------------------------------------------
260 
261  // moved to ConstructSD() for MT compliance
262 
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 void B02ImportanceDetectorConstruction::ConstructSD()
267 {
268 
269  G4SDManager* SDman = G4SDManager::GetSDMpointer();
270  //
271  // Sensitive Detector Name
272  G4String concreteSDname = "ConcreteSD";
273 
274  //------------------------
275  // MultiFunctionalDetector
276  //------------------------
277  //
278  // Define MultiFunctionalDetector with name.
279  G4MultiFunctionalDetector* MFDet =
280  new G4MultiFunctionalDetector(concreteSDname);
281  SDman->AddNewDetector( MFDet ); // Register SD to SDManager
282 
283  G4String fltName,particleName;
284  G4SDParticleFilter* neutronFilter =
285  new G4SDParticleFilter(fltName="neutronFilter", particleName="neutron");
286 
287  MFDet->SetFilter(neutronFilter);
288 
289  for (std::vector<G4LogicalVolume *>::iterator it =
290  fLogicalVolumeVector.begin();
291  it != fLogicalVolumeVector.end(); it++){
292  // (*it)->SetSensitiveDetector(MFDet);
293  SetSensitiveDetector((*it)->GetName(), MFDet);
294  }
295 
296  G4String psName;
297  G4PSNofCollision* scorer0 = new G4PSNofCollision(psName="Collisions");
298  MFDet->RegisterPrimitive(scorer0);
299 
300  G4PSNofCollision* scorer1 = new G4PSNofCollision(psName="CollWeight");
301  scorer1->Weighted(true);
302  MFDet->RegisterPrimitive(scorer1);
303 
304  G4PSPopulation* scorer2 = new G4PSPopulation(psName="Population");
305  MFDet->RegisterPrimitive(scorer2);
306 
307  G4PSTrackCounter* scorer3 =
308  new G4PSTrackCounter(psName="TrackEnter",fCurrent_In);
309  MFDet->RegisterPrimitive(scorer3);
310 
311  G4PSTrackLength* scorer4 = new G4PSTrackLength(psName="SL");
312  MFDet->RegisterPrimitive(scorer4);
313 
314  G4PSTrackLength* scorer5 = new G4PSTrackLength(psName="SLW");
315  scorer5->Weighted(true);
316  MFDet->RegisterPrimitive(scorer5);
317 
318  G4PSTrackLength* scorer6 = new G4PSTrackLength(psName="SLWE");
319  scorer6->Weighted(true);
320  scorer6->MultiplyKineticEnergy(true);
321  MFDet->RegisterPrimitive(scorer6);
322 
323  G4PSTrackLength* scorer7 = new G4PSTrackLength(psName="SLW_V");
324  scorer7->Weighted(true);
325  scorer7->DivideByVelocity(true);
326  MFDet->RegisterPrimitive(scorer7);
327 
328  G4PSTrackLength* scorer8 = new G4PSTrackLength(psName="SLWE_V");
329  scorer8->Weighted(true);
330  scorer8->MultiplyKineticEnergy(true);
331  scorer8->DivideByVelocity(true);
332  MFDet->RegisterPrimitive(scorer8);
333 }
334 
335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
336 G4VIStore* B02ImportanceDetectorConstruction::CreateImportanceStore()
337 {
338 
339  G4cout << " B02ImportanceDetectorConstruction:: Creating Importance Store "
340  << G4endl;
341  if (!fPVolumeStore.Size())
342  {
343  G4Exception("B02ImportanceDetectorConstruction::CreateImportanceStore"
344  ,"exampleB02_0001",RunMustBeAborted
345  ,"no physical volumes created yet!");
346  }
347 
348  // creating and filling the importance store
349 
350  // G4IStore *istore = new G4IStore(*fWorldVolume);
351 
352  G4IStore *istore = G4IStore::GetInstance(GetName());
353 
354  G4GeometryCell gWorldVolumeCell(GetWorldVolumeAddress(), 0);
355 
356  G4double imp =1;
357 
358  istore->AddImportanceGeometryCell(1, gWorldVolumeCell);
359 
360  // set importance values and create scorers
361  G4int cell(1);
362  for (cell=1; cell<=18; cell++) {
363  G4GeometryCell gCell = GetGeometryCell(cell);
364  G4cout << " adding cell: " << cell
365  << " replica: " << gCell.GetReplicaNumber()
366  << " name: " << gCell.GetPhysicalVolume().GetName() << G4endl;
367  imp = std::pow(2.0,cell-1);
368 
369  G4cout << "Going to assign importance: " << imp << ", to volume: "
370  << gCell.GetPhysicalVolume().GetName() << G4endl;
371  //x aIstore.AddImportanceGeometryCell(imp, gCell);
372  istore->AddImportanceGeometryCell(imp, gCell.GetPhysicalVolume(), cell);
373  }
374 
375  // creating the geometry cell and add both to the store
376  // G4GeometryCell gCell = GetGeometryCell(18);
377 
378  // create importance geometry cell pair for the "rest"cell
379  // with the same importance as the last concrete cell
380  G4GeometryCell gCell = GetGeometryCell(19);
381  // G4double imp = std::pow(2.0,18);
382  imp = std::pow(2.0,17);
383  istore->AddImportanceGeometryCell(imp, gCell.GetPhysicalVolume(), 19);
384 
385  return istore;
386 
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390 
391 G4VWeightWindowStore *B02ImportanceDetectorConstruction::
392 CreateWeightWindowStore()
393 {
394  G4cout << " B02ImportanceDetectorConstruction:: Creating Importance Store "
395  << G4endl;
396  if (!fPVolumeStore.Size())
397  {
398  G4Exception("B02ImportanceDetectorConstruction::CreateWeightWindowStore"
399  ,"exampleB02_0002",RunMustBeAborted
400  ,"no physical volumes created yet!");
401  }
402 
403  // creating and filling the importance store
404 
405  // G4IStore *istore = new G4IStore(*fWorldVolume);
406 
407  G4WeightWindowStore *wwstore = G4WeightWindowStore::GetInstance(GetName());
408 
409  // create one energy region covering the energies of the problem
410  //
411  std::set<G4double, std::less<G4double> > enBounds;
412  enBounds.insert(1 * GeV);
413  wwstore->SetGeneralUpperEnergyBounds(enBounds);
414 
415  G4int n = 0;
416  G4double lowerWeight =1;
417  std::vector<G4double> lowerWeights;
418 
419  lowerWeights.push_back(1);
420  G4GeometryCell gWorldCell(GetWorldVolumeAddress(),0);
421  wwstore->AddLowerWeights(gWorldCell, lowerWeights);
422 
423  G4int cell(1);
424  for (cell=1; cell<=18; cell++) {
425  G4GeometryCell gCell = GetGeometryCell(cell);
426  G4cout << " adding cell: " << cell
427  << " replica: " << gCell.GetReplicaNumber()
428  << " name: " << gCell.GetPhysicalVolume().GetName() << G4endl;
429 
430  lowerWeight = 1./std::pow(2., n++);
431  G4cout << "Going to assign lower weight: " << lowerWeight
432  << ", to volume: "
433  << gCell.GetPhysicalVolume().GetName() << G4endl;
434  lowerWeights.clear();
435  lowerWeights.push_back(lowerWeight);
436  wwstore->AddLowerWeights(gCell, lowerWeights);
437  }
438 
439  // the remaining part pf the geometry (rest) gets the same
440  // lower weight bound as the last conrete cell
441  //
442 
443  // create importance geometry cell pair for the "rest"cell
444  // with the same importance as the last concrete cell
445  G4GeometryCell gCell = GetGeometryCell(19);
446  wwstore->AddLowerWeights(gCell, lowerWeights);
447 
448  return wwstore;
449 
450 }