289 lines
11 KiB
C++
289 lines
11 KiB
C++
// Geant4 simulation for MuSR
|
|
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
|
|
// DATE : 2008-05
|
|
//
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
#ifndef lem4RootOutput_h
|
|
#define lem4RootOutput_h 1
|
|
|
|
//#include "G4UserRunAction.hh"
|
|
#include "globals.hh"
|
|
#include "G4ThreeVector.hh"
|
|
// ROOT
|
|
#include "TFile.h"
|
|
#include "TTree.h"
|
|
#include "TH1.h"
|
|
#include "TH2.h"
|
|
#include "TVectorD.h"
|
|
//
|
|
#include <map>
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
class lem4RootOutput {
|
|
public:
|
|
lem4RootOutput();
|
|
~lem4RootOutput();
|
|
static lem4RootOutput* GetRootInstance();
|
|
|
|
public:
|
|
void BeginOfRunAction();
|
|
void EndOfRunAction();
|
|
void FillEvent();
|
|
void ClearAllRootVariables();
|
|
// void SetSensitiveDetectorMapping(std::string logivol, int volumeID);
|
|
void SetVolumeIDMapping(std::string logivol, int volumeID);
|
|
G4int ConvertVolumeToID(std::string logivol);
|
|
void SetSpecialSaveVolumeDefined() {boolIsAnySpecialSaveVolumeDefined=true;};
|
|
|
|
|
|
// Getting variables (just for debugging)
|
|
G4double GetDecayPositionZ() {return muDecayPosZ_t;};
|
|
G4double GetDecayTime() {return muDecayTime_t*microsecond;};
|
|
G4double GetTimeInTarget() {return muTargetTime_t*microsecond;};
|
|
|
|
// Setting variables common to the whole event:
|
|
void SetRunID (G4int id) {runID_t = id;};
|
|
void SetEventID (G4int id) {eventID_t = id;};
|
|
void SetDecayDetectorID (std::string detectorName) {muDecayDetID_t = SensDetectorMapping[detectorName];};
|
|
void SetBField (G4double F[6]) {B_t[0]=F[0]/tesla; B_t[1]=F[1]/tesla; B_t[2]=F[2]/tesla;
|
|
B_t[3]=F[3]/tesla; B_t[4]=F[4]/tesla; B_t[5]=F[5]/tesla;};
|
|
void SetDecayPolarisation (G4ThreeVector pol) {muDecayPolX_t=pol.x(); muDecayPolY_t=pol.y(); muDecayPolZ_t=pol.z();};
|
|
void SetDecayPosition (G4ThreeVector pos) {muDecayPosX_t=pos.x()/mm; muDecayPosY_t=pos.y()/mm;
|
|
muDecayPosZ_t=pos.z()/mm;};
|
|
void SetGlobTime (G4double gt) {globalTime_t = gt/microsecond;};
|
|
|
|
void SetDetectorInfo (G4int nDetectors, G4int ID, G4double edep, G4double edep_el, G4double edep_pos,
|
|
G4double edep_gam, G4double edep_mup,G4int nsteps, G4double length, G4double t1,
|
|
G4double t2, G4double x, G4double y, G4double z) {
|
|
if ((nDetectors<0)||(nDetectors>=(det_nMax-1))) {
|
|
G4cout<<"SERIOUS ERROR! lem4RootOutput.hh: nDetectors="<<nDetectors<<
|
|
" is larger than det_nMax="<<det_nMax<<G4endl;
|
|
return;
|
|
}
|
|
else {
|
|
det_n=nDetectors+1; det_ID[nDetectors]=ID; det_edep[nDetectors]=edep/MeV;
|
|
det_edep_el[nDetectors]=edep_el/MeV; det_edep_pos[nDetectors]=edep_pos/MeV;
|
|
det_edep_gam[nDetectors]=edep_gam/MeV; det_edep_mup[nDetectors]=edep_mup/MeV;
|
|
det_nsteps[nDetectors]=nsteps; det_length[nDetectors]=length/mm;
|
|
det_time_start[nDetectors]=t1/microsecond; det_time_end[nDetectors]=t2/microsecond;
|
|
det_x[nDetectors]=x/mm;det_y[nDetectors]=y/mm;det_z[nDetectors]=z/mm;
|
|
// G4cout << "Saved to root: eventID="<<eventID_t<<"; det_ID="<<det_ID[nDetectors]
|
|
// <<"; det_edep="<<det_edep[nDetectors]<<G4endl;
|
|
}
|
|
//G4cout << "eventID="<<eventID_t<<"; det_n="<<det_n<<"; det_ID="<<det_ID[nDetectors]<<"; "
|
|
// << "Deposited E (el, pos, gam, mu)="<<det_edep_el[nDetectors]<<", "
|
|
// << det_edep_pos[nDetectors]<<", "<< det_edep_gam[nDetectors]<<", "
|
|
// << det_edep_mup[nDetectors]<<G4endl;
|
|
//G4cout << " det_time_start="<<det_time_start[nDetectors]
|
|
// << " det_time_end="<<det_time_end[nDetectors] << G4endl;
|
|
};
|
|
|
|
void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z,
|
|
G4double px, G4double py, G4double pz) ;
|
|
|
|
void SetInitialMuonParameters(G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz,
|
|
G4double xpolaris, G4double ypolaris, G4double zpolaris) {
|
|
muIniPosX_t=x; muIniPosY_t=y; muIniPosZ_t=z;
|
|
muIniMomX_t=px; muIniMomY_t=py; muIniMomZ_t=pz;
|
|
muIniPolX_t=xpolaris; muIniPolY_t=ypolaris; muIniPolZ_t=zpolaris;
|
|
}
|
|
|
|
|
|
// void SetSubtrackInfo (G4int iSubtrack, G4int partID, G4int trID, G4int volID, G4double E,
|
|
// G4double xx, G4double yy, G4double zz, G4double xxPost, G4double yyPost, G4double zzPost,
|
|
// G4double Eki, G4double stepL, G4int nrOfSubst) {
|
|
// if ((iSubtrack<0)||(iSubtrack>=maxNSubTracks)) {
|
|
// G4cout<<"SERIOUS ERROR! lem4RootOutput.hh: iSubtrack="<<iSubtrack<<
|
|
// " is larger than maxNSubTracks="<<maxNSubTracks<<G4endl;
|
|
// return;
|
|
// }
|
|
// else {
|
|
// nSubTracks=iSubtrack+1;
|
|
// particleID_t[iSubtrack]=partID;
|
|
// trackID_t[iSubtrack]=trID;
|
|
// logicalVolumeID_t[iSubtrack]=volID;
|
|
// edep_t[iSubtrack]=E/MeV;
|
|
// x_t[iSubtrack]=xx/mm;
|
|
// y_t[iSubtrack]=yy/mm;
|
|
// z_t[iSubtrack]=zz/mm;
|
|
// post_x_t[iSubtrack]=xxPost/mm;
|
|
// post_y_t[iSubtrack]=yyPost/mm;
|
|
// post_z_t[iSubtrack]=zzPost/mm;
|
|
// kineticEnergy_t[iSubtrack]=Eki/MeV;
|
|
// stepLength_t[iSubtrack]=stepL/mm;
|
|
// nrOfSubsteps_t[iSubtrack]=nrOfSubst;
|
|
// }
|
|
// };
|
|
void SetPolInTarget(G4ThreeVector pol) {muTargetPolX_t=pol.x(); muTargetPolY_t=pol.y(); muTargetPolZ_t=pol.z();};
|
|
void SetTimeInTarget(G4double time) {muTargetTime_t = time/microsecond;};
|
|
void SetInitialPositronMomentum(G4ThreeVector mom) {posIniMomx_t=mom.x(); posIniMomy_t=mom.y(); posIniMomz_t=mom.z();};
|
|
void SetDecayTime(G4double time) {muDecayTime_t=time/microsecond;};
|
|
void SetFieldValue(G4double value) {fieldValue_t=value/tesla;};
|
|
void SetNrFieldNomVal(G4int n) {nFieldNomVal = n;}
|
|
void SetFieldNomVal(G4int i, G4double value);
|
|
G4int GetNrOfVolumes() {return det_nMax;}
|
|
|
|
void StoreGeantParameter(Int_t i, Double_t value) {
|
|
if (i<maxNGeantParameters) { GeantParametersD[i]=value; }
|
|
else {G4cout<<"lem4RootOutput.hh::StoreGeantParameter: index="<<i<<" out of range"
|
|
<<" (maxNGeantParameters=" <<maxNGeantParameters<<")"<<G4endl;}
|
|
};
|
|
|
|
TH2F *htest1, *htest2;
|
|
TH1F *htest3, *htest4, *htest5, *htest6, *htest7, *htest8;
|
|
|
|
public:
|
|
static G4bool store_runID;
|
|
static G4bool store_eventID;
|
|
static G4bool store_BFieldAtDecay;
|
|
static G4bool store_muIniPosX;
|
|
static G4bool store_muIniPosY;
|
|
static G4bool store_muIniPosZ;
|
|
static G4bool store_muIniMomX;
|
|
static G4bool store_muIniMomY;
|
|
static G4bool store_muIniMomZ;
|
|
static G4bool store_muIniPolX;
|
|
static G4bool store_muIniPolY;
|
|
static G4bool store_muIniPolZ;
|
|
static G4bool store_muDecayDetID;
|
|
static G4bool store_muDecayPosX;
|
|
static G4bool store_muDecayPosY;
|
|
static G4bool store_muDecayPosZ;
|
|
static G4bool store_muDecayTime;
|
|
static G4bool store_muDecayPolX;
|
|
static G4bool store_muDecayPolY;
|
|
static G4bool store_muDecayPolZ;
|
|
static G4bool store_muTargetTime;
|
|
static G4bool store_muTargetPolX;
|
|
static G4bool store_muTargetPolY;
|
|
static G4bool store_muTargetPolZ;
|
|
static G4bool store_posIniMomX;
|
|
static G4bool store_posIniMomY;
|
|
static G4bool store_posIniMomZ;
|
|
static G4bool store_globalTime;
|
|
static G4bool store_fieldValue;
|
|
static G4bool store_det_ID;
|
|
static G4bool store_det_edep;
|
|
static G4bool store_det_edep_el;
|
|
static G4bool store_det_edep_pos;
|
|
static G4bool store_det_edep_gam;
|
|
static G4bool store_det_edep_mup;
|
|
static G4bool store_det_nsteps;
|
|
static G4bool store_det_length;
|
|
static G4bool store_det_start;
|
|
static G4bool store_det_end;
|
|
static G4bool store_det_x;
|
|
static G4bool store_det_y;
|
|
static G4bool store_det_z;
|
|
static G4bool store_fieldNomVal;
|
|
|
|
|
|
private:
|
|
TFile* rootFile;
|
|
TTree* rootTree;
|
|
static lem4RootOutput* pointerToRoot;
|
|
static const Int_t maxNGeantParameters=30;
|
|
Double_t GeantParametersD[maxNGeantParameters]; // parameters transfered from GEANT to Root
|
|
// 0 ... fieldOption: 0 ... no field, 1 ... uniform, 2 ... gaussian, 3 ... from table
|
|
// 1 ... fieldValue: intensity of the magnetic field
|
|
// 2 ... minimum of the generated decay time of the muon (in microsecond)
|
|
// 3 ... maximum of the generated decay time of the muon (in microsecond)
|
|
// 4 ... muon mean life time (in microsecond)
|
|
|
|
// Variables common to the whole event:
|
|
Int_t runID_t;
|
|
Int_t eventID_t;
|
|
Double_t B_t[6];
|
|
Double_t muIniPosX_t, muIniPosY_t, muIniPosZ_t;
|
|
Double_t muIniMomX_t, muIniMomY_t, muIniMomZ_t;
|
|
Double_t muIniPolX_t, muIniPolY_t, muIniPolZ_t;
|
|
Int_t muDecayDetID_t;
|
|
Double_t muDecayPolX_t, muDecayPolY_t, muDecayPolZ_t;
|
|
Double_t muTargetPolX_t, muTargetPolY_t, muTargetPolZ_t;
|
|
Double_t muDecayPosX_t, muDecayPosY_t, muDecayPosZ_t;
|
|
Double_t muDecayTime_t;
|
|
Double_t posIniMomx_t, posIniMomy_t, posIniMomz_t;
|
|
Double_t globalTime_t, muTargetTime_t;
|
|
Double_t fieldValue_t;
|
|
|
|
public:
|
|
static const Int_t maxNFieldnNominalValues=30;
|
|
private:
|
|
Int_t nFieldNomVal;
|
|
Double_t fieldNomVal[maxNFieldnNominalValues];
|
|
|
|
// Variables for a particle in a given detector within the event
|
|
public:
|
|
static const Int_t maxNSubTracks=30;
|
|
private:
|
|
// Int_t nSubTracks;
|
|
// // Char_t particleName_t[6];
|
|
// Int_t particleID_t[maxNSubTracks];
|
|
// G4int trackID_t[maxNSubTracks];
|
|
// Double_t edep_t[maxNSubTracks];
|
|
// Double_t x_t[maxNSubTracks];
|
|
// Double_t y_t[maxNSubTracks];
|
|
// Double_t z_t[maxNSubTracks];
|
|
// Double_t post_x_t[maxNSubTracks];
|
|
// Double_t post_y_t[maxNSubTracks];
|
|
// Double_t post_z_t[maxNSubTracks];
|
|
// // Char_t logicalVolume_t[13];
|
|
// Int_t logicalVolumeID_t[maxNSubTracks];
|
|
// Double_t kineticEnergy_t[maxNSubTracks];
|
|
// Double_t stepLength_t[maxNSubTracks];
|
|
// Int_t nrOfSubsteps_t[maxNSubTracks];
|
|
|
|
|
|
// Variables for the activity inside a given detector
|
|
public:
|
|
static const Int_t det_nMax=100; // must be by 1 higher than the real number of detector "hits", because
|
|
// else the detector nr. 0 is counted (0 is used if no
|
|
// SensDetectorMapping correspond to a given logical volume).
|
|
private:
|
|
G4int det_n;
|
|
G4int det_ID[det_nMax];
|
|
G4double det_edep[det_nMax];
|
|
G4int det_nsteps[det_nMax];
|
|
G4double det_length[det_nMax];
|
|
G4double det_edep_el[det_nMax];
|
|
G4double det_edep_pos[det_nMax];
|
|
G4double det_edep_gam[det_nMax];
|
|
G4double det_edep_mup[det_nMax];
|
|
G4double det_time_start[det_nMax];
|
|
G4double det_time_end[det_nMax];
|
|
G4double det_x[det_nMax];
|
|
G4double det_y[det_nMax];
|
|
G4double det_z[det_nMax];
|
|
|
|
// G4bool firstStepInVolume_t;
|
|
// G4bool lastStepInVolume_t;
|
|
|
|
public:
|
|
static const Int_t save_nMax=1000;
|
|
|
|
private:
|
|
G4int save_n;
|
|
G4int save_detID[save_nMax];
|
|
G4int save_particleID[save_nMax];
|
|
G4double save_ke[save_nMax];
|
|
G4double save_x[save_nMax];
|
|
G4double save_y[save_nMax];
|
|
G4double save_z[save_nMax];
|
|
G4double save_px[save_nMax];
|
|
G4double save_py[save_nMax];
|
|
G4double save_pz[save_nMax];
|
|
|
|
G4bool boolIsAnySpecialSaveVolumeDefined;
|
|
|
|
std::map<std::string,int> SensDetectorMapping;
|
|
};
|
|
|
|
//G4int lem4RootOutput::maxNSubTracks=30;
|
|
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
#endif
|