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