Geant4  10.00.p03
XrayFluoAnalysisManager.cc
Go to the documentation of this file.
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 // $Id: XrayFluoAnalysisManager.cc
28 // GEANT4 tag $Name: xray_fluo-V03-02-00
29 //
30 // Author: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
31 //
32 // History:
33 // -----------
34 // 20 Jul 2007 A.Mantero, code cleaning and update. Replaced histos with tuple
35 // 11 Jul 2003 A.Mantero, code cleaning / Plotter-XML addiction
36 // Sep 2002 A.Mantero, AIDA3.0 Migration
37 // 06 Dec 2001 A.Pfeiffer updated for singleton
38 // 30 Nov 2001 Guy Barrand : migrate to AIDA-2.2.
39 // 28 Nov 2001 Elena Guardincerri Created
40 //
41 // -------------------------------------------------------------------
42 #ifdef G4ANALYSIS_USE
43 
44 #include "G4VProcess.hh"
45 #include <fstream>
46 #include "G4ios.hh"
47 #include "XrayFluoAnalysisManager.hh"
48 #include "G4Step.hh"
49 #include "XrayFluoDetectorConstruction.hh"
50 #include "G4VPhysicalVolume.hh"
51 #include "G4SystemOfUnits.hh"
52 
53 XrayFluoAnalysisManager* XrayFluoAnalysisManager::instance = 0;
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 
57 XrayFluoAnalysisManager::XrayFluoAnalysisManager()
58  :outputFileName("xrayfluo"), visPlotter(false), phaseSpaceFlag(false), physicFlag (false), persistencyType("root"),
59  deletePersistencyFile(true), gunParticleEnergies(0), gunParticleTypes(0), analysisFactory(0), tree(0), treeDet(0), histogramFactory(0), plotter(0)
60 {
61  //creating the messenger
62  analisysMessenger = new XrayFluoAnalysisMessenger(this);
63 
64  // Hooking an AIDA compliant analysis system.
65  // creating Analysis factroy, necessary to create/manage analysis
66  analysisFactory = AIDA_createAnalysisFactory();
67 
68  CreatePersistency(outputFileName,persistencyType);
69 
70  if(analysisFactory) {
71 
72 
73  // Creating the plotter factory
74 
75  if (visPlotter){
76  plotterFactory = analysisFactory->createPlotterFactory();
77  }
78  }
79 
80  G4cout << "XrayFluoAnalysisManager created" << G4endl;
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
85 XrayFluoAnalysisManager::~XrayFluoAnalysisManager()
86 {
87  delete histogramFactory;
88  histogramFactory=0;
89 
90  delete analysisFactory;
91  analysisFactory = 0;
92 
93  delete plotterFactory;
94  plotterFactory=0;
95 
96  delete tree;
97 
98  if ( gunParticleEnergies ) delete gunParticleEnergies;
99  gunParticleEnergies = 0;
100  if ( gunParticleTypes ) delete gunParticleTypes;
101  gunParticleTypes = 0;
102 
103  delete instance;
104 
105  G4cout << "XrayFluoAnalysisManager delete" << G4endl;
106 
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
111 XrayFluoAnalysisManager* XrayFluoAnalysisManager::getInstance()
112 
113 {
114  if (instance == 0) {instance = new XrayFluoAnalysisManager;}
115  return instance;
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 
120 void XrayFluoAnalysisManager::CreatePersistency(G4String fileName,G4String thePersistencyType,
121  G4bool readOnly, G4bool createNew)
122 {
123 
124  if (tree) {
125  delete tree;
126  tree = 0;
127  }
128  if (treeDet) {
129  delete treeDet;
130  treeDet =0;
131  }
132  if (histogramFactory) {
133  delete histogramFactory;
134  histogramFactory = 0;
135  }
136  if (tupleFactory) {
137  delete tupleFactory;
138  tupleFactory = 0;
139  }
140 
141  if (tupleDetFactory) {
142  delete tupleDetFactory;
143  tupleDetFactory = 0;
144  }
145  G4String fileNameDet = fileName+"Detector";
146 
147  AIDA::ITreeFactory* treeFactory = analysisFactory->createTreeFactory();
148  if(treeFactory) {
149  if (thePersistencyType == "hbook") {
150  fileName = fileName + ".hbk";
151  fileNameDet = fileNameDet + ".hbk";
152  }
153  else if (thePersistencyType == "xml"){
154  fileName = fileName + ".xml";
155  fileNameDet = fileNameDet + ".xml";
156  }
157 
158  if (phaseSpaceFlag) {
159 
160  tree = treeFactory->create(fileName,thePersistencyType,readOnly,createNew); // output file
161  treeDet = treeFactory->create(fileNameDet,thePersistencyType,readOnly,createNew); // output file
162  if(analysisFactory) {
163 
164  tupleDetFactory = analysisFactory->createTupleFactory(*treeDet);
165  tupleFactory = analysisFactory->createTupleFactory(*tree);
166 
167  }
168  }
169 
170  // trees to be stored in case of phase-space production
171  else {
172 
173  tree = treeFactory->create(fileName,thePersistencyType,readOnly,createNew); // output file
174 
175  if(analysisFactory) {
176 
177  histogramFactory = analysisFactory->createHistogramFactory(*tree);
178 
179  }
180  }
181 
182  delete treeFactory; // Will not delete the ITree.
183  }
184 }
185 
186 
187 void XrayFluoAnalysisManager::book()
188 {
189 
190  if (phaseSpaceFlag) {
191  /*
192  // Book clouds
193 
194  cloud_1 = histogramFactory->createCloud1D("Gamma Exting Sample","ciao!",-1);
195  cloud_2 = histogramFactory->createCloud1D("Gamma Incident on the detector","ciao!",-1);
196  cloud_3 = histogramFactory->createCloud1D("Electrons Exiting the Sample","ciao!",-1);
197  beamCloud = histogramFactory->createCloud1D("Incident Radiation Spectrum","ciao!",-1);
198  */
199  // Book output Tuple
200 
201  // Book tuple
202  std::vector<std::string> columnNames;
203  columnNames.push_back("particle");
204  columnNames.push_back("energies");
205  columnNames.push_back("momentumTheta");
206  columnNames.push_back("momentumPhi");
207  columnNames.push_back("processes");
208  columnNames.push_back("material");
209  columnNames.push_back("origin");
210  columnNames.push_back("depth");
211 
212  std::vector<std::string> columnTypes;
213  columnTypes.push_back("int");
214  columnTypes.push_back("double");
215  columnTypes.push_back("double");
216  columnTypes.push_back("double");
217  columnTypes.push_back("int"); // useful for hbk
218  columnTypes.push_back("int"); // useful for hbk
219  columnTypes.push_back("int");
220  columnTypes.push_back("double");
221 
222  tupleFluo = tupleFactory->create("10", "Total Tuple", columnNames, columnTypes, "");
223  assert(tupleFluo);
224 
225  }
226 
227  else {
228  // Book histograms
229 
230  histo_1 = histogramFactory->createHistogram1D("1","Energy Deposit", 500,0.,10.); //20eV def.
231  histo_2 = histogramFactory->createHistogram1D("2","Gamma born in the sample", 100,0.,10.);
232  histo_3 = histogramFactory->createHistogram1D("3","Electrons born in the sample", 100,0.,10.);
233  histo_4 = histogramFactory->createHistogram1D("4","Gammas leaving the sample", 300,0.,10.);
234  histo_5 = histogramFactory->createHistogram1D("5","Electrons leaving the sample ",200000 ,0.,10.0); // .05 eV def.
235  histo_6 = histogramFactory->createHistogram1D("6","Gammas reaching the detector", 100,0.,10.);
236  histo_7 = histogramFactory->createHistogram1D("7","Spectrum of the incident particles", 100,0.,10.);
237  histo_8 = histogramFactory->createHistogram1D("8","Protons reaching the detector", 100,0.,10.);
238  histo_9 = histogramFactory->createHistogram1D("9","Protons leaving the sample", 100,0.,10.);
239 
240  // Debugging-purpose Histos
241 
242  //histo_10 = histogramFactory->createHistogram1D("10","Photon Origin", 4,0.,3.);
243  //histo_11 = histogramFactory->createHistogram1D("11","Spectrum from LowEnPhotoELectric", 300,0.,10.);
244  //histo_12 = histogramFactory->createHistogram1D("12","Spectrum From the other processes (unknown)", 300,0.,10.);
245 
246  // IHistogram2D* histo_20 = histogramFactory->
247  //create2D("20","Phi, Theta",80 ,-3.14,3.14,80,0.,3.14);
248  }
249 }
250 
251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
252 
253 void XrayFluoAnalysisManager::LoadGunData(G4String fileName, G4bool raileighFlag) {
254 
255  // gunParticleEnergies = new std::vector<G4double>;
256  // gunParticleTypes = new std::vector<G4String>;
257 
258  G4String ext = fileName.substr(fileName.size()-3,fileName.size()-1);
259  G4String thePersistencyType;
260 
261  if (ext == "xml") {
262 
263  thePersistencyType = "xml";
264  }
265  else if (ext == "hbk") {
266 
267  thePersistencyType = "hbook";
268  }
269 
270  gunParticleEnergies = new std::vector<G4double>;
271  gunParticleTypes = new std::vector<G4String>;
272 
273  // dal punto di vista della programmazione e' un po' una porcata.....
274 
275  AIDA::ITreeFactory* treeDataFactory = analysisFactory->createTreeFactory();
276  AIDA::ITree* treeData = treeDataFactory->create(fileName,thePersistencyType,false, false); // input file
277  AIDA::IManagedObject* mo = treeData->find("10");
278  assert(mo);
279  G4cout <<"mo Name:" << mo->name()<< G4endl;
280  AIDA::ITuple* tupleData = dynamic_cast<AIDA::ITuple*>(mo);
281 
282  //AIDA::ITuple* tupleData = (AIDA::ITuple*) treeData->find("10");
283 
284  assert(tupleData);
285  // da riscrivere un pochino melgio. Del resto serve per pochissmo, poi non serve piu'
286 
287 
288  tupleData->start();
289  G4int processesIndex = tupleData->findColumn("processes");
290  G4int energyIndex = tupleData->findColumn("energies");
291  G4int particleIndex = tupleData->findColumn("particle");
292 
293  while (tupleData->next()) {
294  if (raileighFlag ^ (!raileighFlag && (tupleData->getInt(processesIndex) == 1 ||
295  tupleData->getInt(processesIndex) == 2)) ) {
296  gunParticleEnergies->push_back(tupleData->getDouble(energyIndex)*MeV);
297  if (tupleData->getInt(particleIndex) == 1 ) gunParticleTypes->push_back("gamma");
298  if (tupleData->getInt(particleIndex) == 0 ) gunParticleTypes->push_back("e-");
299  }
300 
301  }
302  G4cout << "Maximum mumber of events: "<< gunParticleEnergies->size() <<G4endl;
303 
304 
305 }
306 
307 std::vector<G4double>* XrayFluoAnalysisManager::GetEmittedParticleEnergies() {
308 
309  return gunParticleEnergies;
310 
311 }
312 
313 std::vector<G4String>* XrayFluoAnalysisManager::GetEmittedParticleTypes() {
314 
315  return gunParticleTypes;
316 
317 }
318 
319 
320 void XrayFluoAnalysisManager::finish()
321 {
322 
323  if (tupleFluo) {ExtractData();};
324 
325  if(plotter) {
326  // Wait for the keyboard return to avoid destroying the plotter window too quickly.
327  G4cout << "Press <ENTER> to exit" << G4endl;
328  G4cin.get();
329  plotter->hide(); //hide plotter windows, but doesn't delete plotter
330  }
331 
332  if(tree) {
333  tree->commit(); // Write histos in file.
334  tree->close();
335  }
336  if (treeDet) {
337  treeDet->commit();
338  treeDet->close();
339  }
340 
341  deletePersistencyFile = false;
342 
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346 
347 void XrayFluoAnalysisManager::InitializePlotter()
348 {
349 
350  // If speciefied (visPlotter=true),
351  // a window for the visulizaton of partial results is created
352  if(plotterFactory && visPlotter && deletePersistencyFile)
353  {
354  plotter = plotterFactory->create();
355  // Set the page title :
356  plotter->setParameter("pageTitle","XrayFluo");
357  }
358 
359  if(plotter && visPlotter) {
360  plotter->show(); // shows plotter window
361  }
362 
363 }
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
365 
366 void XrayFluoAnalysisManager::PlotCurrentResults()
367 {
368  if(plotter) {
369  if (phaseSpaceFlag){
370 
371 
372  // altra porcata di programmazione e cmq mi serve un istogramma. Da rivedere il tutto.
373  // Plotting the spectrum of Gamma exiting the sample - cloud1
374  AIDA::ICloud1D& cloud = *cloud_1;
375  AIDA::IFilter* filterGamma = tupleFactory->createFilter(
376  " Particle == std::string(\"gamma\") ");
377  AIDA::IEvaluator* evaluatorEnergy = tupleFactory->createEvaluator("Energies");
378 
379  // di nuovo, del resto serve solo per un attimo.. da correggere cmq.
380 
381  filterGamma->initialize(*tupleFluo);
382  evaluatorEnergy->initialize(*tupleFluo);
383  tupleFluo->project(cloud,*evaluatorEnergy,*filterGamma);
384 
385  plotter->currentRegion().plot( cloud, "Exiting Gammas " );
386  plotter->refresh();
387  }
388 
389  else{
390  // Plotting the Detector Energy Deposit - histo_1
391  AIDA::IHistogram1D& histo1p = *histo_1;
392  plotter->currentRegion().plot( histo1p, "Detector Energy Deposition" );
393  plotter->refresh();
394  }
395 
396  }
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 G4bool XrayFluoAnalysisManager::GetDeletePersistencyFileFlag()
401 {
402  return deletePersistencyFile;
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 
407 void XrayFluoAnalysisManager::SetPhysicFlag(G4bool val)
408 {
409  physicFlag = val;
410 }
411 
412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413 
414 
415 
416 
417 
418 void XrayFluoAnalysisManager::analyseStepping(const G4Step* aStep)
419 {
420 
421  if (phaseSpaceFlag){
422 
423 
424  G4String particleType="";
425  G4String parentProcess="";
426  G4ThreeVector momentum(0.,0.,0.);
427  G4double particleEnergy=0;
428  G4String sampleMaterial="";
429  G4double particleDepth=0;
430  G4int isBornInTheSample=0;
431  XrayFluoDetectorConstruction* detector = XrayFluoDetectorConstruction::GetInstance();
432 
433  // Select volume from which the step starts
434  if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
435  G4ThreeVector creationPos = aStep->GetTrack()->GetVertexPosition();
436  // qua serve ancora una selezione: deve essere interno al sample
437  //codice "a mano" per controllare il volume di orgine della traccia
438  /*
439  G4Double sampleXplus = detector->GetSampleXplusFaceCoord();
440  G4Double sampleXminus = detector->GetSampleXminusFaceCoord();
441  G4Double sampleYplus = detector->GetSampleYplusFaceCoord();
442  G4Double sampleYminus = detector->GetSampleYminusFaceCoord();
443  G4Double sampleZplus = detector->GetSampleYplusFaceCoord();
444  G4Double sampleZminus = detector->GetSampleYminusFaceCoord();
445 
446  G4bool ZsampleCheck = (creationPos.z()<=sampleZminus) && (creationPos.z() >= sampleZplus);
447  G4bool XsampleCheck = (creationPos.x()<=sampleZminus) && (creationPos.x() >= sampleXplus);
448  G4bool YsampleCheck = (creationPos.y()<=sampleZminus) && (creationPos.y() >= sampleYplus);
449  */
450  G4VPhysicalVolume* creationPosVolume = detector->GetGeometryNavigator()->LocateGlobalPointAndSetup(creationPos);
451 
452  // if physicflag is on, I record all the data of all the particle born in the sample.
453  // If it is off (but phase space production is on) I record data
454  // only for particles creted inside the sample and exiting it.
455  if (physicFlag ^ (!physicFlag &&
456  (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
457  //ZsampleCheck && XsampleCheck && YsampleCheck)) {
458  ) ) {
459  //
460  //
461  // G4cout << "physicFlag: "<< physicFlag << G4endl
462  // << "NextVolume: "<< aStep->GetTrack()->GetNextVolume()->GetName() << G4endl;
463 
464  // extracting information needed
465  particleType = aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetParticleName();
466  momentum = aStep->GetTrack()->GetDynamicParticle()->GetMomentum();
467  particleEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
468  if (creationPosVolume->GetName() == "Sample" ) {
469  isBornInTheSample = 1;
470  particleDepth = creationPosVolume->GetLogicalVolume()->GetSolid()->DistanceToOut(creationPos, G4ThreeVector(0,0,-1));
471  }
472  else {
473  particleDepth = -1;
474  }
475  // metodo per ottenere la profondita' "a mano"
476  // particleDepth = std::abs(creationPos.z() - detector->GetSampleIlluminatedFaceCoord());
477 
478  G4int parent;
479  if(aStep->GetTrack()->GetCreatorProcess()){
480  parentProcess = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
481  parent = 5;
482  // hack for HBK
483  if (parentProcess == "") parent = 0;
484  if (parentProcess == "IONI") parent = 1;
485  if (parentProcess == "LowEnPhotoElec") parent = 2;
486  if (parentProcess == "Transportation") parent = 3;// should never happen
487  if (parentProcess == "initStep") parent = 4;// should never happen
488  }
489  else {
490  parentProcess = "Not Known -- (primary generator + Rayleigh)";
491  parent = 5;
492  }
493  G4int sampleMat=0;
494  if(aStep->GetTrack()){
495  sampleMaterial = aStep->GetTrack()->GetMaterial()->GetName();
496  if (sampleMaterial == ("Dolorite" || "Anorthosite" || "Mars1" || "IceBasalt" || "HPGe")) sampleMat=1;
497  }
498  // filling tuple
499 
500  // G4cout<< particleType << G4endl;
501  G4int part = -1 ;
502  if (particleType == "gamma") part =1;
503  if (particleType == "e-") part = 0;
504  if (particleType == "proton") part = 2;
505 
506 
507  tupleFluo->fill(0,part);
508  tupleFluo->fill(1,particleEnergy);
509  tupleFluo->fill(2,momentum.theta());
510  tupleFluo->fill(3,momentum.phi());
511  tupleFluo->fill(4,parent); //hacked to be used with hbk
512  tupleFluo->fill(5,sampleMat); //hacked to be used with hbk
513  tupleFluo->fill(6,isBornInTheSample);
514  tupleFluo->fill(7,particleDepth);
515 
516  tupleFluo->addRow();
517  }
518  }
519  }
520 
521 
522  // Normal behaviour, without creation of phase space
523  else {
524 
525  G4double gammaAtTheDetPre=0;
526  G4double protonsAtTheDetPre=0;
527  G4double gammaLeavingSample=0;
528  //G4double gammaLeavingSamplePhi=0;
529  //G4double gammaLeavingSampleTheta=0;
530 
531  G4double eleLeavingSample=0;
532  G4double protonsLeavSam=0;
533  G4double gammaBornInSample=0;
534  G4double eleBornInSample=0;
535 
536  // Filling the histograms with data, passing thru stepping.
537 
538 
539  // Select volume from wich the step starts
540  if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
541  // Select volume from wich the step ends
542  if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" ) {
543  // Select the particle type
544  if ((aStep->GetTrack()->GetDynamicParticle()
545  ->GetDefinition()-> GetParticleName()) == "gamma" )
546 
547  {
548 
549 
550  // Control histos, used for debugging purpose.
551 
552  // if(aStep->GetTrack()->GetCreatorProcess()){
553  // G4String process = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
554  // G4cout << "Il processo di origine e' " << process << G4endl;
555  // if(histo_10) {
556  // histo_10->fill(1.);
557  // gammaLeavingSample =
558  // (aStep->GetPreStepPoint()->GetKineticEnergy());
559  // if(histo_11) {
560  // histo_11->fill(gammaLeavingSample/keV);
561 
562  // }
563  // }
564  // }
565  // else {
566  // //G4cout << "Sembra che non ci sia un processo d'origine" << G4endl;
567  // if(histo_10) {
568  // histo_10->fill(2.);
569 
570  // gammaLeavingSample =
571  // (aStep->GetPreStepPoint()->GetKineticEnergy());
572  // if(histo_12) {
573  // histo_12->fill(gammaLeavingSample/keV);
574 
575  // }
576 
577  // }
578  // }
579 
580 
581  // //
582  // Filling Histos //
583  // //
584 
585 
586  gammaLeavingSample =
587  (aStep->GetPreStepPoint()->GetKineticEnergy());
588  if(histo_4) {
589  histo_4->fill(gammaLeavingSample/keV);
590 
591  }
592 
593  // Other debugging purpose histos
594  /* gammaLeavingSamplePhi =
595  (aStep->GetPreStepPoint()->GetMomentumDirection().phi());
596  G4cout << "questo e' Phi: " << gammaLeavingSamplePhi << G4endl;
597  gammaLeavingSampleTheta =
598  (aStep->GetPreStepPoint()->GetMomentumDirection().theta());
599  G4cout << "questo e' Theta: " << gammaLeavingSampleTheta << G4endl;
600  if(histo_20) {
601  // G4cout << "histo_20 esiste" << G4endl;
602  histo_20->fill(gammaLeavingSamplePhi,gammaLeavingSampleTheta,1.);
603  } */
604 
605 
606  }
607  }
608  }
609 
610 
611 
612  if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
613 
614  if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" )
615  {
616  if ((aStep->GetTrack()->GetDynamicParticle()
617  ->GetDefinition()-> GetParticleName()) == "e-" )
618  {
619  eleLeavingSample = (aStep->GetPreStepPoint()->GetKineticEnergy());
620  if(histo_5) {
621  histo_5->fill(eleLeavingSample/keV);
622  }
623  }
624  else if ((aStep->GetTrack()->GetDynamicParticle()
625  ->GetDefinition()-> GetParticleName()) == "proton" )
626  {
627  protonsLeavSam = (aStep->GetPreStepPoint()->GetKineticEnergy());
628  if(histo_9) {
629  histo_9->fill(protonsLeavSam/keV);
630  }
631  }
632 
633  }
634  }
635 
636  // electrons from the detector -- Debugging
637 
638 // if((aStep->GetTrack()->GetDynamicParticle()
639 // ->GetDefinition()-> GetParticleName()) == "e-" ){
640 
641 // if(1== (aStep->GetTrack()->GetCurrentStepNumber())){
642 
643 // if(0 != aStep->GetTrack()->GetParentID()){
644 // if(aStep->GetTrack()->GetVolume()->GetName() == "HPGeDetector")
645 
646 // if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" ){
647 
648 // if(aStep->GetTrack()->GetCreatorProcess()){
649 // G4String process = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
650 // G4cout << "Origin Porcess Name: " << process << G4endl;
651 // if(histo_10) {
652 // histo_10->fill(1.);
653 // gammaLeavingSample =
654 // (aStep->GetPreStepPoint()->GetKineticEnergy());
655 // if(histo_11) {
656 // histo_11->fill(gammaLeavingSample/keV);
657 
658 // }
659 // }
660 // }
661 // else {
662 // G4cout << "No Origin Process Name" << G4endl;
663 // if(histo_10) {
664 // histo_10->fill(2.);
665 
666 // gammaLeavingSample =
667 // (aStep->GetPreStepPoint()->GetKineticEnergy());
668 // if(histo_12) {
669 // histo_12->fill(gammaLeavingSample/keV);
670 
671 // }
672 
673 // }
674 // }
675 // }
676 // }
677 
678 // }
679 // }
680 
681 
682 
683 
684  if((aStep->GetTrack()->GetDynamicParticle()
685  ->GetDefinition()-> GetParticleName()) == "gamma" )
686 
687  {if(1== (aStep->GetTrack()->GetCurrentStepNumber()))
688 
689  {if(0 != aStep->GetTrack()->GetParentID())
690 
691  {if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
692  {
693  gammaBornInSample = (aStep->GetPreStepPoint()->GetKineticEnergy());
694  if(histo_2) {
695  histo_2->fill(gammaBornInSample/keV);
696  }
697  }
698  }
699  }
700  }
701  if((aStep->GetTrack()->GetDynamicParticle()
702  ->GetDefinition()-> GetParticleName()) == "e-" )
703 
704  {if(1== (aStep->GetTrack()->GetCurrentStepNumber()))
705 
706  {if(0 != aStep->GetTrack()->GetParentID())
707 
708  {if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
709 {
710  eleBornInSample = (aStep->GetPreStepPoint()->GetKineticEnergy());
711  if(histo_3) {
712  histo_3->fill(eleBornInSample/keV);
713  }
714  }
715  }
716  }
717  }
718 
719  if(aStep->GetTrack()->GetNextVolume()){
720 
721  //if(aStep->GetTrack()->GetVolume()->GetName() == "World"){
722 
723  if(aStep->GetTrack()->GetNextVolume()->GetName() == "HPGeDetector")
724 
725  {
726  if ((aStep->GetTrack()->GetDynamicParticle()
727  ->GetDefinition()-> GetParticleName()) == "gamma" )
728  {
729 
730  gammaAtTheDetPre =
731  (aStep->GetPreStepPoint()->GetKineticEnergy());
732  if(histo_6) {
733  histo_6->fill( gammaAtTheDetPre/keV);
734  }
735  }
736  else if ((aStep->GetTrack()->GetDynamicParticle()
737  ->GetDefinition()-> GetParticleName()) == "proton" )
738  {
739  protonsAtTheDetPre =
740  (aStep->GetPreStepPoint()->GetKineticEnergy());
741  if(histo_8) {
742  histo_8->fill( protonsAtTheDetPre/keV);
743  }
744  }
745  }
746  //} close of if(aStep->GetTrack()->GetVolume()->GetName() == "World"){
747  }
748  }
749 }
750 
751 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
752 
753 void XrayFluoAnalysisManager::ExtractData(){
754 
755  if (tupleFluo->rows()) {
756 
757 
758  // AIDA::IFilter* filterGamma = tupleFactory->createFilter(" Particle == std::string(\"gamma\")");
759  // AIDA::IFilter* filterEminus = tupleFactory->createFilter(" Particle == std::string(\"e-\")");
760 
761 
762  AIDA::IFilter* filterAngle = tupleFactory->createFilter(
763 
764 "(momentumPhi >= (220. * (3.1415926/180.) )) && (momentumPhi <= (230. * (3.1415926/180.) )) && (momentumTheta >= (130. * (3.1415926/180.) )) && (momentumTheta <= (140. * (3.1415926/180.) ))" );
765 
766 
767  // filterGamma ->initialize(*tupleFluo);
768  // filterEminus ->initialize(*tupleFluo);
769  filterAngle->initialize(*tupleFluo);
770 
771  // Create IEvaluator and initialize it to this ITuple
772  // AIDA::IEvaluator* evaluatorEnergy = tupleFactory->createEvaluator("Energies");
773  // evaluatorEnergy->initialize(*tupleFluo);
774 
775  // tupleFluo->project(*cloud_1,*evaluatorEnergy,*filterGamma);
776  // tupleFluo->project(*cloud_2,*evaluatorEnergy,*filterAngle);
777  // tupleFluo->project(*cloud_3,*evaluatorEnergy,*filterEminus);
778 
779  tupleDetFluo = tupleDetFactory->createFiltered("1", *tupleFluo, *filterAngle);
780  assert(tupleDetFluo);
781  }
782 }
783 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
784 
785 
786 void XrayFluoAnalysisManager::analyseEnergyDep(G4double energyDep)
787 {
788 
789  // Filling of Energy Deposition, called by XrayFluoEventAction
790 
791  histo_1->fill(energyDep/keV);
792 
793 }
794 
795 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
796 
797 void XrayFluoAnalysisManager::analysePrimaryGenerator(G4double energy)
798 {
799 
800  // Filling of energy spectrum histogram of the primary generator
801 
802  if (phaseSpaceFlag){
803  // beamCloud->fill(energy/keV);
804  }
805  else {
806  histo_7->fill(energy/keV);
807  }
808 }
809 
810 
811 // not used -- Created for future development
812 
813 void XrayFluoAnalysisManager::SetOutputFileName(G4String newName)
814 {
815 
816  outputFileName = newName;
817 }
818 
819 void XrayFluoAnalysisManager::SetOutputFileType(G4String newType)
820 {
821 
822  persistencyType = newType;
823 }
824 
825 #endif
G4Track::GetParentID
G4int GetParentID() const
MeV
static const double MeV
Definition: G4SIunits.hh:193
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4StepPoint::GetStepStatus
G4StepStatus GetStepStatus() const
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:176
assert
#define assert(x)
Definition: mymalloc.cc:1309
XrayFluoAnalysisManager.hh
G4Step.hh
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
XrayFluoDetectorConstruction::GetInstance
static XrayFluoDetectorConstruction * GetInstance()
Definition: XrayFluoDetectorConstruction.cc:156
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4Track::GetNextVolume
G4VPhysicalVolume * GetNextVolume() const
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4Track::GetCreatorProcess
const G4VProcess * GetCreatorProcess() const
G4cin
#define G4cin
Definition: G4ios.hh:60
G4cout
G4GLOB_DLL std::ostream G4cout
G4Track::GetCurrentStepNumber
G4int GetCurrentStepNumber() const
G4StepPoint::GetPhysicalVolume
G4VPhysicalVolume * GetPhysicalVolume() const
G4VPhysicalVolume::GetName
const G4String & GetName() const
G4bool
bool G4bool
Definition: G4Types.hh:79
XrayFluoDetectorConstruction.hh
G4Step
Definition: G4Step.hh:76
G4VProcess::GetProcessName
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4Track::GetVertexPosition
const G4ThreeVector & GetVertexPosition() const
G4Track::GetMaterial
G4Material * GetMaterial() const
fGeomBoundary
Definition: G4StepStatus.hh:54
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4VProcess.hh
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:156
G4Step::GetPostStepPoint
G4StepPoint * GetPostStepPoint() const
G4ios.hh
G4Navigator::LocateGlobalPointAndSetup
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:118
G4Track::GetVolume
G4VPhysicalVolume * GetVolume() const
G4endl
#define G4endl
Definition: G4ios.hh:61
instance
static MCTruthManager * instance
Definition: MCTruthManager.cc:43
keV
static const double keV
Definition: G4SIunits.hh:195
G4VPhysicalVolume.hh
G4StepPoint::GetKineticEnergy
G4double GetKineticEnergy() const
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4VSolid::DistanceToOut
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4Step::GetTrack
G4Track * GetTrack() const
G4LogicalVolume::GetSolid
G4VSolid * GetSolid() const
G4DynamicParticle::GetMomentum
G4ThreeVector GetMomentum() const