xray_fluorescence/src/XrayFluoAnalysisManager.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // Author: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
0029 //
0030 // History:
0031 // -----------
0032 // 20 Jul 2007 A.Mantero, code cleaning and update. Replaced histos with tuple
0033 // 11 Jul 2003 A.Mantero, code cleaning / Plotter-XML addiction
0034 //    Sep 2002 A.Mantero, AIDA3.0 Migration
0035 // 06 Dec 2001 A.Pfeiffer updated for singleton
0036 // 30 Nov 2001 Guy Barrand : migrate to AIDA-2.2.
0037 // 28 Nov 2001 Elena Guardincerri     Created
0038 //
0039 // -------------------------------------------------------------------
0040 #include "G4AnalysisManager.hh"
0041 #include "G4RootAnalysisReader.hh"
0042
0043 #include "G4VProcess.hh"
0044 #include "XrayFluoAnalysisManager.hh"
0045 #include "G4Step.hh"
0046 #include "XrayFluoDetectorConstruction.hh"
0047 #include "G4VPhysicalVolume.hh"
0048 #include "G4Gamma.hh"
0049 #include "G4Electron.hh"
0050 #include "G4Proton.hh"
0051 #include "G4SystemOfUnits.hh"
0052
0053 using G4AnalysisReader = G4RootAnalysisReader;
0054
0055 XrayFluoAnalysisManager* XrayFluoAnalysisManager::instance = 0;
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0058
0059 namespace {
0060   //Mutex to acquire access to singleton instance check/creation
0061   G4Mutex instanceMutex = G4MUTEX_INITIALIZER;
0062   //Mutex to acquire accss to histograms creation/access
0063   //It is also used to control all operations related to histos
0064   //File writing and check analysis
0065   G4Mutex dataManipulationMutex = G4MUTEX_INITIALIZER;
0066 }
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0069
0070 XrayFluoAnalysisManager::XrayFluoAnalysisManager()
0071   :outputFileName("xrayfluo"), phaseSpaceFlag(false), physicFlag (false),
0072    gunParticleEnergies(0), gunParticleTypes(0)
0073 {
0074   dataLoaded = false;
0075   fParticleEnergyAndTypeIndex = 0;
0076
0077   //creating the messenger
0078   analisysMessenger = new XrayFluoAnalysisMessenger(this);
0079
0080   //Instantiate the analysis manager
0081   G4AnalysisManager::Instance();
0082
0083   G4cout << "XrayFluoAnalysisManager created" << G4endl;
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0087
0088 XrayFluoAnalysisManager::~XrayFluoAnalysisManager()
0089 {
0090   if ( gunParticleEnergies ) delete gunParticleEnergies;
0091   gunParticleEnergies = 0;
0092   if ( gunParticleTypes ) delete gunParticleTypes;
0093   gunParticleTypes = 0;
0094
0095   delete instance;
0096   instance = 0;
0097 }
0098
0099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0100
0101 XrayFluoAnalysisManager* XrayFluoAnalysisManager::getInstance()
0102
0103 {
0104   G4AutoLock l(&instanceMutex);
0105   if (instance == 0) {instance = new XrayFluoAnalysisManager;}
0106   return instance;
0107 }
0108
0109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0110
0111 void XrayFluoAnalysisManager::book()
0112 {
0113   G4AutoLock l(&dataManipulationMutex);
0114   // Get analysis manager
0115   G4AnalysisManager* man = G4AnalysisManager::Instance();
0116   man->SetDefaultFileType("root");
0117   // Open an output file
0118   man->OpenFile(outputFileName);
0119   man->SetVerboseLevel(1);
0120   man->SetFirstHistoId(1);
0121   man->SetFirstNtupleId(1);
0122
0123   G4cout << "Open output file: " << outputFileName << G4endl;
0124
0125   if (phaseSpaceFlag)
0126     {
0127       // Book output Tuple (ID = 1)
0128       man->CreateNtuple("101","OutputNTuple");
0129       man->CreateNtupleIColumn("particle");
0130       man->CreateNtupleDColumn("energies");
0131       man->CreateNtupleDColumn("momentumTheta");
0132       man->CreateNtupleDColumn("momentumPhi");
0133       man->CreateNtupleIColumn("processes");
0134       man->CreateNtupleIColumn("material");
0135       man->CreateNtupleIColumn("origin");
0136       man->CreateNtupleDColumn("depth");
0137       man->FinishNtuple();
0138       G4cout << "Created ntuple for phase space" << G4endl;
0139     }
0140   else {
0141     // Book histograms
0142     man->CreateH1("h1","Energy Deposit", 500,0.,10.); //20eV def.
0143     man->CreateH1("h2","Gamma born in the sample", 100,0.,10.);
0144     man->CreateH1("h3","Electrons  born in the sample", 100,0.,10.);
0145     man->CreateH1("h4","Gammas leaving the sample", 300,0.,10.);
0146     man->CreateH1("h5","Electrons leaving the sample ",200000 ,0.,10.0); // .05 eV def.
0147     man->CreateH1("h6","Gammas reaching the detector", 100,0.,10.);
0148     man->CreateH1("h7","Spectrum of the incident particles", 100,0.,10.);
0149     man->CreateH1("h8","Protons reaching the detector", 100,0.,10.);
0150     man->CreateH1("h9","Protons leaving the sample", 100,0.,10.);
0151     G4cout << "Created histos" << G4endl;
0152   }
0153 }
0154
0155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0156
0157 void XrayFluoAnalysisManager::LoadGunData(G4String fileName, G4bool raileighFlag)
0158 {
0159   G4AutoLock l(&dataManipulationMutex);
0160
0161   //Do not allow more than one thread to trigger the file reading
0162   if (dataLoaded)
0163     return;
0164
0165   // Get analysis reader manager
0166   G4AnalysisReader* analysisReader = G4AnalysisReader::Instance();
0167   analysisReader->SetVerboseLevel(1);
0168
0169   //This is for testing purposes
0170   G4int ntupleId = analysisReader->GetNtuple("101",fileName);
0171   G4cout << "Got ntupleId: " << ntupleId << G4endl;
0172
0173   gunParticleEnergies = new std::vector<G4double>;
0174   gunParticleTypes = new std::vector<G4String>;
0175
0176   G4int particleID, processID;
0177   G4double energy;
0178   analysisReader->SetNtupleIColumn("processes", processID);
0179   analysisReader->SetNtupleDColumn("energies", energy);
0180   analysisReader->SetNtupleIColumn("particles", particleID);
0181
0182   while (analysisReader->GetNtupleRow() )
0183     {
0184       if (raileighFlag ^ (!raileighFlag && (processID == 1 ||
0185                         processID == 2)) )
0186     {
0187       gunParticleEnergies->push_back(energy*MeV);
0188       if (particleID == 1)
0189         gunParticleTypes->push_back("gamma");
0190       else if (particleID == 0)
0191         gunParticleTypes->push_back("e-");
0192     }
0193     }
0194
0195   G4cout << "Maximum mumber of events: "<< gunParticleEnergies->size() <<G4endl;
0196
0197   dataLoaded= true;
0198 }
0199
0200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0201
0202 const std::pair<G4double,G4String> XrayFluoAnalysisManager::GetEmittedParticleEnergyAndType()
0203 {
0204   G4AutoLock l(&dataManipulationMutex);
0205   std::pair<G4double,G4String> result;
0206
0207   if(fParticleEnergyAndTypeIndex < (G4int) gunParticleEnergies->size())
0208     {
0209       G4double energy = gunParticleEnergies->at(fParticleEnergyAndTypeIndex);
0210       G4String name = gunParticleTypes->at(fParticleEnergyAndTypeIndex);
0211       result.first = energy;
0212       result.second = name;
0213     }
0214
0215   fParticleEnergyAndTypeIndex++;
0216   return result;
0217 }
0218
0219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0220
0221
0222 void XrayFluoAnalysisManager::finish()
0223 {
0224   G4AutoLock l(&dataManipulationMutex);
0225   G4cout << "Going to save histograms" << G4endl;
0226   // Save histograms
0227   G4AnalysisManager* man = G4AnalysisManager::Instance();
0228   man->Write();
0229   man->CloseFile();
0230 }
0231
0232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0233
0234
0235 void XrayFluoAnalysisManager::SetPhysicFlag(G4bool val)
0236 {
0237   physicFlag = val;
0238 }
0239
0240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0241
0242 void XrayFluoAnalysisManager::analyseStepping(const G4Step* aStep)
0243 {
0244   G4AutoLock l(&dataManipulationMutex);
0245   G4AnalysisManager* man = G4AnalysisManager::Instance();
0246
0247   if (phaseSpaceFlag){
0248     G4ParticleDefinition* particleType= 0;
0249     G4String parentProcess="";
0250     G4ThreeVector momentum(0.,0.,0.);
0251     G4double particleEnergy=0;
0252     G4String sampleMaterial="";
0253     G4double particleDepth=0;
0254     G4int isBornInTheSample=0;
0255     XrayFluoDetectorConstruction* detector = XrayFluoDetectorConstruction::GetInstance();
0256
0257     // Select volume from which the step starts
0258     if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
0259       G4ThreeVector creationPos = aStep->GetTrack()->GetVertexPosition();
0260       // qua serve ancora una selezione: deve essere interno al sample
0261       //codice "a mano" per controllare il volume di orgine della traccia
0262
0263       G4VPhysicalVolume* creationPosVolume = detector->GetGeometryNavigator()->LocateGlobalPointAndSetup(creationPos);
0264
0265       // if physicflag is on, I record all the data of all the particle born in the sample.
0266       // If it is off (but phase space production is on) I record data
0267       // only for particles creted inside the sample and exiting it.
0268       if (physicFlag ^ (!physicFlag &&
0269             (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
0270             ))
0271     {
0272       // extracting information needed
0273       particleType = aStep->GetTrack()->GetDynamicParticle()->GetDefinition();
0274       momentum = aStep->GetTrack()->GetDynamicParticle()->GetMomentum();
0275       particleEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
0276       if (creationPosVolume->GetName() == "Sample" ) {
0277         isBornInTheSample = 1;
0278         particleDepth = creationPosVolume->GetLogicalVolume()->GetSolid()
0279           ->DistanceToOut(creationPos, G4ThreeVector(0,0,-1));
0280       }
0281       else {
0282         particleDepth = -1;
0283       }
0284       // metodo per ottenere la profondita' "a mano"
0285       //    particleDepth = std::abs(creationPos.z() - detector->GetSampleIlluminatedFaceCoord());
0286
0287       G4int parent=0;
0288       if(aStep->GetTrack()->GetCreatorProcess())
0289         {
0290           parentProcess = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
0291           parent = 5;
0292           // hack for HBK
0293           if (parentProcess == "") parent = 0;
0294           if (parentProcess == "ioni") parent = 1;
0295           if (parentProcess == "LowEnPhotoElec") parent = 2;
0296           if (parentProcess == "Transportation") parent = 3;// should never happen
0297           if (parentProcess == "initStep") parent = 4;// should never happen
0298         }
0299       else {
0300         parentProcess = "Not Known -- (primary generator + Rayleigh)";
0301         parent = 5;
0302       }
0303       G4int sampleMat=0;
0304       if(aStep->GetTrack()){
0305         sampleMaterial = aStep->GetTrack()->GetMaterial()->GetName();
0306         if (sampleMaterial == "Dolorite"
0307           || sampleMaterial == "Anorthosite"
0308           || sampleMaterial == "Mars1"
0309           || sampleMaterial == "IceBasalt"
0310           || sampleMaterial == "HPGe") sampleMat=1;
0311       }
0312
0313       G4int part = -1 ;
0314       if (particleType == G4Gamma::Definition()) part =1;
0315       if (particleType == G4Electron::Definition()) part = 0;
0316       if (particleType == G4Proton::Definition()) part = 2;
0317
0318       //Fill ntuple
0319       man->FillNtupleIColumn(1,0, part);
0320       man->FillNtupleDColumn(1,1,particleEnergy);
0321       man->FillNtupleDColumn(1,2,momentum.theta());
0322       man->FillNtupleDColumn(1,3,momentum.phi());
0323       man->FillNtupleIColumn(1,4,parent);
0324       man->FillNtupleIColumn(1,5,sampleMat);
0325       man->FillNtupleIColumn(1,6,isBornInTheSample);
0326       man->FillNtupleDColumn(1,7,particleDepth);
0327       man->AddNtupleRow(1);
0328
0329     }
0330     }
0331   }
0332
0333   // Normal behaviour, without creation of phase space
0334   else
0335     {
0336
0337       // Filling the histograms with data, passing thru stepping.
0338
0339       // Select volume from wich the step starts
0340       if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
0341     // Select volume from wich the step ends
0342     if(aStep->GetTrack()->GetNextVolume()->GetName() == "World" ) {
0343       // Select the particle type
0344       if ((aStep->GetTrack()->GetDynamicParticle()
0345            ->GetDefinition()) == G4Gamma::Definition() )
0346
0347         {
0348           G4double gammaLeavingSample =
0349         aStep->GetPreStepPoint()->GetKineticEnergy();
0350           man->FillH1(4,gammaLeavingSample/keV);
0351
0352         }
0353       else if ((aStep->GetTrack()->GetDynamicParticle()
0354             ->GetDefinition()) == G4Electron::Definition() )
0355         {
0356           G4double eleLeavingSample =
0357         aStep->GetPreStepPoint()->GetKineticEnergy();
0358           man->FillH1(5,eleLeavingSample/keV);
0359         }
0360       else if ((aStep->GetTrack()->GetDynamicParticle()
0361             ->GetDefinition()) == G4Proton::Definition() )
0362         {
0363           G4double
0364         protonsLeavSam = aStep->GetPreStepPoint()->GetKineticEnergy();
0365           man->FillH1(9,protonsLeavSam/keV);
0366         }
0367     }
0368       }
0369
0370       if((aStep->GetTrack()->GetDynamicParticle()
0371       ->GetDefinition()) == G4Gamma::Definition() )
0372         {
0373       if(aStep->GetTrack()->GetCurrentStepNumber() == 1)
0374         {
0375           if(aStep->GetTrack()->GetParentID() != 0)
0376         {
0377           if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
0378             {
0379               G4double gammaBornInSample =
0380             aStep->GetPreStepPoint()->GetKineticEnergy();
0381               man->FillH1(2,gammaBornInSample/keV);
0382             }
0383         }
0384         }
0385     }
0386       else if ((aStep->GetTrack()->GetDynamicParticle()
0387         ->GetDefinition() ) == G4Electron::Definition() )
0388     {
0389       if(aStep->GetTrack()->GetCurrentStepNumber() == 1)
0390         {
0391           if(aStep->GetTrack()->GetParentID() != 0)
0392         {
0393           if(aStep->GetTrack()->GetVolume()->GetName() == "Sample")
0394             {
0395               G4double eleBornInSample =
0396             aStep->GetPreStepPoint()->GetKineticEnergy();
0397               man->FillH1(3,eleBornInSample/keV);
0398             }
0399         }
0400         }
0401     }
0402
0403
0404       if(aStep->GetTrack()->GetNextVolume())
0405     {
0406
0407
0408       if(aStep->GetTrack()->GetNextVolume()->GetName() == "HPGeDetector")
0409
0410         {
0411           if ((aStep->GetTrack()->GetDynamicParticle()
0412            ->GetDefinition()) == G4Gamma::Definition() )
0413         {
0414
0415           G4double gammaAtTheDetPre =
0416             aStep->GetPreStepPoint()->GetKineticEnergy();
0417           man->FillH1(6,gammaAtTheDetPre/keV);
0418         }
0419           else if ((aStep->GetTrack()->GetDynamicParticle()
0420             ->GetDefinition() ) == G4Proton::Definition() )
0421         {
0422           G4double protonsAtTheDetPre =
0423             aStep->GetPreStepPoint()->GetKineticEnergy();
0424           man->FillH1(8,protonsAtTheDetPre);
0425         }
0426         }
0427     }
0428     }
0429 }
0430
0431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0432
0433 void XrayFluoAnalysisManager::analyseEnergyDep(G4double energyDep)
0434 {
0435   G4AutoLock l(&dataManipulationMutex);
0436   // Filling of Energy Deposition, called by XrayFluoEventAction
0437   G4AnalysisManager* man = G4AnalysisManager::Instance();
0438
0439   if (!phaseSpaceFlag)
0440     man->FillH1(1,energyDep/keV);
0441
0442 }
0443
0444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0445
0446 void XrayFluoAnalysisManager::analysePrimaryGenerator(G4double energy)
0447 {
0448   G4AutoLock l(&dataManipulationMutex);
0449   // Filling of energy spectrum histogram of the primary generator
0450   G4AnalysisManager* man = G4AnalysisManager::Instance();
0451   if (!phaseSpaceFlag)
0452     man->FillH1(7,energy/keV);
0453
0454 }
0455
0456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0457
0458 void XrayFluoAnalysisManager::SetOutputFileName(G4String newName)
0459 {
0460
0461   outputFileName = newName;
0462 }
0463
0464