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