LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

7.2.2011 Kamil Sedlak

Browse Source 
1) Constant fraction discriminator added to the analysis of the 
     optical photon signal. 
  2) Bug fixes in the optical photon simulation.
  3) A few other small improvements.
This commit is contained in:
sedlak 2011-02-07 16:07:38 +00:00
parent d11e19a7d8
commit a5e99ed164
10 changed files with 10389 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/musrSim.pdf
View File

Binary file not shown.

73
doc/musrSim.tex
View File

@ -126,7 +126,8 @@ in the macro file:
\end{itemize}
By default, the output of the simulation is written out in the subdirectory ``data'' with
the name ``musr\_RUNNUMBER.root''.
the name ``musr\_RUNNUMBER.root''. The default ``data'' directory can be changed by the
command ``/musr/command rootOutputDirectoryName \emph{dirName}''.
(Note that the execution of the simulation can be terminated gently by creating a file ``RUNNUMBER.stop'' in the working directory.)
\section{Conventions}
@ -747,10 +748,72 @@ in order to simulate optical photons:
E.g.\ REFLECTIVITY or EFFICIENCY for a given surface may be assigned this way.
\end{description}
One has to assign a non-zero EFFICIENCY and a REFLECTIVITY smaller than 1 to a boundary surface
between the scintillator and sensitive device (e.g.\ an APD).
\item {\bf /musr/command OPSA \emph{subcommand} \emph{parameters ...}}\\
Different commands related to Optical Photon Signal Analysis (OPSA):\\
\begin{description}
\item{\bf signalSeparationTime \emph{timeSeparation}} \\
If two subsequent optical photons arrive to the same active detector
in time smaller than \emph{timeSeparation}, they will form
the same ``optical photon signal''. Otherwise (i.e.\ if there is no photon
detected for time larger than \emph{timeSeparation}) the next arriving photon
will start a new ``optical photon signal''. See the similarity
to the command ``/musr/command signalSeparationTime \emph{timeSeparation}''
for the deposited energy signals.
\item{\bf OPSAhist \emph{nBins} \emph{min} \emph{max}} \\
Defines ``OPSA'' histograms -- i.e.\ histograms that
contain the time distribution of the arival of optical photons.\\
\emph{nBins} ... number of bins of the histogram\\
\emph{min} ... minimum of the x-axis of the histogram\\
\emph{max}... maximum of the x-axis the histogram\\
In fact three kinds of histograms are created -- see command
``/musr/command OPSA eventsForOPSAhistos \emph{eventID} \emph{detID}''.
\item{\bf eventsForOPSAhistos \emph{eventID} \emph{detID}}\\
Defines event number and detector ID, for which histograms of the
OPSA will be stored. If \emph{detID} is set to zero, all detectors will be
stored for the given event.
The naming convention for the histograms is the following:
OPSAhist\_\emph{eventID}\_\emph{detID}\_\emph{n} ,\\
where \emph{n} is the number of the ``optical photon signal'' in the
event \emph{eventID} in the detector \emph{detID}.\\
There are other two histograms, namely
``OPSAshape\_\emph{eventID}\_\emph{detID}\_\emph{n}'', which shows
the signal from OPSAhist convoluted with the responce function
of the optical detection device as e.g.\ G-APD, and
``OPSA\_CFD\_\emph{eventID}\_\emph{detID}\_\emph{n}'', which shows
the signal from OPSAshape after the constant fraction discriminator.
\item{\bf pulseShapeArray \emph{pulseShapeFileName}}\\
\emph{pulseShapeFileName} is the name of the file that contains responce
function (pulse shape) of a single cell (single photon) detected by
the photosensitive detector. The data format is very strict:\\
\% ... comments \\
first column ... time in picosecond (it has to be: 0, 1, 2, ... iPSmax
where iPSmax is smaller than 10000)\\
second column ... amplitude of the APD response function\\
{\bf Internally in musrSim, the data read from this file are interpolated to the
centers of the bins of the histograms defined by
``/musr/command OPSA OPSAhist ''}.
\item{\bf CFD \emph{a1} \emph{delay} \emph{timeShiftOffset}}\\
\item{\bf minNrOfDetectedPhotons}\\
\item{\bf photonFractions}\\
\end{description}
\end{description}
\subsection*{Tips and tricks}
\begin{itemize}
\item One has to assign a non-zero EFFICIENCY and a REFLECTIVITY smaller than 1 to a boundary surface
between the scintillator and sensitive device (e.g.\ an APD).
\item To get a perfectly reflecting surface, one can use REFLECTIVITY=1; EFFICIENCY=0;
\emph{surfaceType}=dielectric\_metal, \emph{surfaceFinish}=polished, \emph{surfaceModel}=unified.
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Some other parameters}
%
@ -760,6 +823,10 @@ between the scintillator and sensitive device (e.g.\ an APD).
\emph{nrOfEvents} (int) -- number of events to be simulated.\\
(This is a default Geant4 command, which has to be specified in any simulation run).
\item{\bf /musr/command rootOutputDirectoryName \emph{dirName}} \\
(default: /musr/command rootOutputDirectoryName data )\\
Specify the (sub)directory to which the output file will be written out.
\item{\bf /musr/command logicalVolumeToBeReweighted mu \emph{logicalVolume} \emph{weight} }\\
(default: not defined; no reweighting is done unless explicitly requested by this command.) \\
Events can be reweighted by this command. If muon {\bf stops and decays} in the

10
include/musrRootOutput.hh
View File

@ -80,7 +80,8 @@ class musrRootOutput {
G4int idVolVertex, G4int idProcVertex, G4int idTrackVertex, G4int particleID) ;
void SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeA,
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast);
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast,
G4double timeCFD, G4double amplCFD);
void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z, G4double time,
G4double px, G4double py, G4double pz, G4double polx, G4double poly, G4double polz) ;
@ -124,7 +125,7 @@ class musrRootOutput {
else {G4cout<<"musrRootOutput.hh::StoreGeantParameter: index="<<i<<" out of range"
<<" (maxNGeantParameters=" <<maxNGeantParameters<<")"<<G4endl;}
};
void setRootOutputDirectoryName(char dirName[1000]);
TH2F *htest1, *htest2;
TH1F *htest3, *htest4, *htest5, *htest6, *htest7, *htest8;
@ -222,6 +223,8 @@ class musrRootOutput {
static G4bool store_odet_timeD;
static G4bool store_odet_timeE;
static G4bool store_odet_timeLast;
static G4bool store_odet_timeCFD;
static G4bool store_odet_amplCFD;
static G4int oldEventNumberInG4EqEMFieldWithSpinFunction;
@ -230,6 +233,7 @@ class musrRootOutput {
TTree* rootTree;
static musrRootOutput* pointerToRoot;
static const Int_t maxNGeantParameters=30;
char rootOutputDirectoryName[1000];
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
@ -327,6 +331,8 @@ class musrRootOutput {
G4double odet_timeD[odet_nMax];
G4double odet_timeE[odet_nMax];
G4double odet_timeLast[odet_nMax];
G4double odet_timeCFD[odet_nMax];
G4double odet_amplCFD[odet_nMax];
public:
static const Int_t save_nMax=1000;

39
include/musrScintSD.hh
View File

@ -34,12 +34,14 @@ class G4HCofThisEvent;
class signalInfo {
public:
signalInfo(G4int id, G4int nP, G4double tFirst,
G4double tA, G4double tB, G4double tC, G4double tD, G4double tE, G4double tLast) {
G4double tA, G4double tB, G4double tC, G4double tD, G4double tE, G4double tLast,
G4double tCFD, G4double aCFD) {
detID=id; nPhot=nP; timeFirst=tFirst;
timeA=tA; timeB=tB; timeC=tC; timeD=tD; timeE=tE; timeLast=tLast;}
timeA=tA; timeB=tB; timeC=tC; timeD=tD; timeE=tE; timeLast=tLast;
timeCFD=tCFD; amplCFD=aCFD;}
~signalInfo() {}
void transferDataToRoot(musrRootOutput* myRootOut, G4int nn) {
myRootOut->SetOPSAinfo(nn,detID,nPhot,timeFirst,timeA,timeB,timeC,timeD,timeE,timeLast);
myRootOut->SetOPSAinfo(nn,detID,nPhot,timeFirst,timeA,timeB,timeC,timeD,timeE,timeLast,timeCFD,amplCFD);
}
private:
@ -55,6 +57,8 @@ class signalInfo {
G4double timeD;
G4double timeE;
G4double timeLast;
G4double timeCFD;
G4double amplCFD;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -69,17 +73,26 @@ class musrScintSD : public G4VSensitiveDetector
void Initialize(G4HCofThisEvent*);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
void EndOfEvent(G4HCofThisEvent*);
void EndOfRun();
// Optical Photon Signal Analysis (OPSA)
void Set_OPSA_minNrOfDetectedPhotons(G4int val) {OPSA_minNrOfDetectedPhotons=val;}
void Set_OPSA_SignalSeparationTime(G4double val) {OPSA_signalSeparationTime=val;}
void Set_OPSA_frac(G4double a, G4double b, G4double c, G4double d, G4double e)
{OPSA_fracA=a; OPSA_fracB=b; OPSA_fracC=c; OPSA_fracD=d; OPSA_fracE=e;}
void AddEventIDToMultimapOfEventIDsForOPSAhistos (G4int ev_ID, G4int detector_ID)
{bool_multimapOfEventIDsForOPSAhistosEXISTS=true; multimapOfEventIDsForOPSAhistos.insert(std::pair<G4int,G4int>(ev_ID,detector_ID));}
void SetOPSAhistoBinning(Int_t nBins, Double_t min, Double_t max) {OPSAhistoNbin=nBins; OPSAhistoMin=min; OPSAhistoMax=max;}
void Set_OPSA_CFD(G4double a1, G4double delay, G4double timeShiftOffset)
{OPSA_CFD_a1=a1; OPSA_CFD_delay=delay; OPSA_CFD_timeShiftOffset = timeShiftOffset;}
void AddEventIDToMultimapOfEventIDsForOPSAhistos (G4int ev_ID, G4int detector_ID) {
bool_multimapOfEventIDsForOPSAhistosEXISTS=true;
multimapOfEventIDsForOPSAhistos.insert(std::pair<G4int,G4int>(ev_ID,detector_ID));
}
void SetOPSAhistoBinning(Int_t nBins, Double_t min, Double_t max) {
OPSAhistoNbin=nBins; OPSAhistoMin=min; OPSAhistoMax=max;
OPSAhistoBinWidth=(max-min)/nBins; OPSAhistoBinWidth1000=OPSAhistoBinWidth*1000;
}
void ProcessOpticalPhoton(G4Step*);
void EndOfEvent_OptiacalPhotons();
void ReadInPulseShapeArray(const char* filename);
private:
static musrScintSD* pointer;
@ -111,7 +124,19 @@ class musrScintSD : public G4VSensitiveDetector
multimapOfEventIDsForOPSAhistos_Type multimapOfEventIDsForOPSAhistos;
TH1D* OPSAhisto;
Int_t OPSAhistoNbin;
Double_t OPSAhistoMin, OPSAhistoMax;
Double_t OPSAhistoMin, OPSAhistoMax, OPSAhistoBinWidth, OPSAhistoBinWidth1000;
typedef std::map<std::string,TH1D*> mapOfOPSAsumHistograms_Type;
mapOfOPSAsumHistograms_Type mapOfOPSAsumHistograms;
mapOfOPSAsumHistograms_Type mapOfOPSAsum0Histograms;
TH1D* OPSAhistoSUM;
TH1D* OPSAhistoSUM0;
TH1D* OPSAshape;
G4bool bool_pulseShapeExists;
G4int iPSmax;
G4double pulseShapeArray[10001];
TH1D* OPSA_CFD;
G4double OPSA_CFD_a1,OPSA_CFD_delay, OPSA_CFD_timeShiftOffset;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

10014
run/APDpulseShapeFile.txt Normal file
View File

File diff suppressed because it is too large Load Diff

24
src/musrDetectorConstruction.cc
View File

@ -645,21 +645,21 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical type \""<<type<<"\" not found!"<<G4endl;
G4cout<<" (Maybe it exists in Geant4 but was not added to musrSim ==> update musrDetectorConstruction.cc file!)"<< G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
G4cout<<" "<<line<<G4endl;
exit(1);
}
if ((OpticalModel==0)&&(strcmp(model,"glisur")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical surface model \""<<model<<"\" not found!"<<G4endl;
G4cout<<" (Maybe it exists in Geant4 but was not added to musrSim ==> update musrDetectorConstruction.cc file!)"<< G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
G4cout<<" "<<line<<G4endl;
exit(1);
}
if ((OpticalFinish==0)&&(strcmp(finish,"polished")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical surface finish \""<<finish<<"\" not found!"<<G4endl;
G4cout<<" (Maybe it exists in Geant4 but was not added to musrSim ==> update musrDetectorConstruction.cc file!)"<< G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
G4cout<<" "<<line<<G4endl;
exit(1);
}
@ -726,6 +726,16 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
sscanf(&line[0],"%*s %*s %*s %d %lf %lf",&nBins,&min,&max);
myMusrScintSD -> SetOPSAhistoBinning(nBins,min*nanosecond,max*nanosecond);
}
else if (strcmp(varName,"pulseShapeArray")==0) {
char fileName[500];
sscanf(&line[0],"%*s %*s %*s %s",fileName);
myMusrScintSD -> ReadInPulseShapeArray(fileName);
}
else if (strcmp(varName,"CFD")==0) {
double a1, delay, timeShiftOffset;
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf",&a1,&delay,&timeShiftOffset);
myMusrScintSD -> Set_OPSA_CFD(a1,delay,timeShiftOffset);
}
}
}
@ -922,6 +932,14 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if (strcmp(tmpString1,"rootOutputDirectoryName")==0){
char rootOutDirName[1000];
sscanf(&line[0],"%*s %*s %s",rootOutDirName);
// Cut out the character '/' at the end of the directory name, if it was specified in the name
if (rootOutDirName[strlen(rootOutDirName)-1]=='/') rootOutDirName[strlen(rootOutDirName)-1]='\0';
myRootOutput->setRootOutputDirectoryName(rootOutDirName);
}
else if (strcmp(tmpString1,"storeOnlyEventsWithHitInDetID")==0){
G4int variable;
sscanf(&line[0],"%*s %*s %d",&variable);

26
src/musrRootOutput.cc
View File

@ -33,6 +33,7 @@ musrRootOutput::musrRootOutput() {
pointerToRoot=this;
boolIsAnySpecialSaveVolumeDefined=false;
nFieldNomVal=0;
strcpy(rootOutputDirectoryName,"data");
ProcessIDMapping["DecayWithSpin"]=1;
ProcessIDMapping["eIoni"]=2;
@ -157,6 +158,8 @@ G4bool musrRootOutput::store_odet_timeC = true;
G4bool musrRootOutput::store_odet_timeD = true;
G4bool musrRootOutput::store_odet_timeE = true;
G4bool musrRootOutput::store_odet_timeLast = true;
G4bool musrRootOutput::store_odet_timeCFD = true;
G4bool musrRootOutput::store_odet_amplCFD = true;
G4int musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction=-1;
@ -167,7 +170,8 @@ void musrRootOutput::BeginOfRunAction() {
G4cout << "musrRootOutput::BeginOfRunAction() Defining the Root tree and branches:"<<G4endl;
G4int tmpRunNr=(G4RunManager::GetRunManager())->GetCurrentRun()->GetRunID();
char RootOutputFileName[200];
sprintf(RootOutputFileName, "data/musr_%i.root", tmpRunNr);
// sprintf(RootOutputFileName, "data/musr_%i.root", tmpRunNr);
sprintf(RootOutputFileName, "%s/musr_%i.root",rootOutputDirectoryName,tmpRunNr);
rootFile=new TFile(RootOutputFileName,"recreate");
rootTree=new TTree("t1","a simple Tree with simple variables");
if (store_runID) {rootTree->Branch("runID",&runID,"runID/I");}
@ -279,7 +283,7 @@ void musrRootOutput::BeginOfRunAction() {
}
if (store_odet_ID || store_odet_nPhot || store_odet_timeFirst || store_odet_timeA || store_odet_timeB ||
store_odet_timeC || store_odet_timeD || store_odet_timeE || store_odet_timeLast)
store_odet_timeC || store_odet_timeD || store_odet_timeE || store_odet_timeLast || store_odet_timeCFD || store_odet_amplCFD)
{rootTree->Branch("odet_n",&odet_n,"odet_n/I");}
if (store_odet_ID) {rootTree->Branch("odet_ID",&odet_ID,"odet_ID[odet_n]/I");}
if (store_odet_nPhot) {rootTree->Branch("odet_nPhot",&odet_nPhot,"odet_nPhot[odet_n]/I");}
@ -289,7 +293,10 @@ void musrRootOutput::BeginOfRunAction() {
if (store_odet_timeC) {rootTree->Branch("odet_timeC",&odet_timeC,"odet_timeC[odet_n]/D");}
if (store_odet_timeD) {rootTree->Branch("odet_timeD",&odet_timeD,"odet_timeD[odet_n]/D");}
if (store_odet_timeE) {rootTree->Branch("odet_timeE",&odet_timeE,"odet_timeE[odet_n]/D");}
if (store_odet_timeLast) {rootTree->Branch("odet_timeLast",&odet_timeLast,"odet_timeLast[odet_n]/D");}
if (store_odet_timeLast) {rootTree->Branch("odet_timeLast",&odet_timeLast,"odet_timeLast[odet_n]/D");}
if (store_odet_timeCFD) {rootTree->Branch("odet_timeCFD",&odet_timeCFD,"odet_timeCFD[odet_n]/D");}
if (store_odet_amplCFD) {rootTree->Branch("odet_amplCFD",&odet_amplCFD,"odet_amplCFD[odet_n]/D");}
// htest1 = new TH1F("htest1","The debugging histogram 1",50,-4.,4.);
// htest2 = new TH1F("htest2","The debugging histogram 2",50,0.,3.142);
@ -527,7 +534,7 @@ void musrRootOutput::SetDetectorInfoVvv (G4int nDetectors,
void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeA,
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast)
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast, G4double timeCFD, G4double amplCFD)
{
if ((nDetectors<0)||(nDetectors>=(odet_nMax-1))) {
char message[200];
@ -546,5 +553,16 @@ void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4
odet_timeD[nDetectors]=timeD/microsecond;
odet_timeE[nDetectors]=timeE/microsecond;
odet_timeLast[nDetectors]=timeLast/microsecond;
odet_timeCFD[nDetectors]=timeCFD/microsecond;
odet_amplCFD[nDetectors]=amplCFD;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::setRootOutputDirectoryName(char dirName[1000]) {
strcpy(rootOutputDirectoryName,dirName);
char message[200];
sprintf(message,"musrRootOutput.cc::setRootOutputDirectoryName: Root output file will be stored in directory %s",dirName);
musrErrorMessage::GetInstance()->musrError(INFO,message,false);
}

3
src/musrRunAction.cc
View File

@ -26,6 +26,7 @@
#include "musrRunAction.hh"
#include "musrEventAction.hh"
#include "musrSteppingAction.hh"
#include "musrScintSD.hh"
#include "G4Run.hh"
#include "musrErrorMessage.hh"
#include "F04GlobalField.hh"
@ -92,6 +93,8 @@ void musrRunAction::BeginOfRunAction(const G4Run* aRun) {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRunAction::EndOfRunAction(const G4Run* aRun) {
musrScintSD::GetInstance()->EndOfRun();
musrRootOutput::GetRootInstance()->StoreGeantParameter(5,aRun->GetNumberOfEvent());
musrRootOutput::GetRootInstance()->EndOfRunAction();

267
src/musrScintSD.cc
View File

@ -54,9 +54,9 @@
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Double_t poissonf(Double_t* x, Double_t* par) {
return par[0]*TMath::Poisson(x[0],par[1]);
}
//Double_t poissonf(Double_t* x, Double_t* par) {
// return par[0]*TMath::Poisson(x[0],par[1]);
//}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -77,6 +77,12 @@ musrScintSD::musrScintSD(G4String name)
OPSAhistoNbin = 100;
OPSAhistoMin =0;
OPSAhistoMax = 10.;
OPSAhistoBinWidth=0.1;
OPSAhistoBinWidth1000=100.;
bool_pulseShapeExists=false;
OPSA_CFD_a1 = -0.2;
OPSA_CFD_delay = 2.;
OPSA_CFD_timeShiftOffset = 0.;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -89,6 +95,8 @@ musrScintSD* musrScintSD::GetInstance() {return pointer;}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
// This "Initialize" method is called at the beginning of each event
// --> perhaps some this code could be executed once per run?
if (verboseLevel>1) G4cout<<"VERBOSE 2: musrScintSD::Initialize\n";
scintCollection = new musrScintHitsCollection
(SensitiveDetectorName,collectionName[0]);
@ -105,6 +113,14 @@ void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
boolIsVvvInfoRequested = myMusrSteppingAction->IsVvvInfoRequested();
myRootOutput = musrRootOutput::GetRootInstance();
// In case of optical photons, delete all optHitDetectorMap* from the previous event (if they exist).
if (musrParameters::boolG4OpticalPhotons) {
for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
delete (it->second);
}
optHitMap.clear();
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -224,7 +240,7 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
}
optHitDetectorMapType* optHitDetectorMap = (*iter).second;
// optHitDetectorMapType optHitDetectorMap = optHitMap[detID];
G4double tmpTime = aStep->GetPreStepPoint()->GetGlobalTime();
G4double tmpTime = aStep->GetPostStepPoint()->GetGlobalTime();
// G4cout<<" tmpTime="<<tmpTime;
// optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,boundaryStatus));
@ -243,6 +259,7 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
//is the photocathode because it is the only one with
//non-zero efficiency
optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,boundaryStatus));
// G4cout<<" tmpTime="<<tmpTime<<G4endl;
// G4cout<<" Detection"<<G4endl;
// break;
// case FresnelReflection:
@ -513,6 +530,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
G4bool boolStoreThisOPSAhist = false;
G4bool boolStoreThisOPSAhistSUMMED = false;
G4int OPSA_detID= it->first;
optHitDetectorMapType* optHitDetectorMap = it->second;
@ -529,6 +547,10 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
if ( (itOPSAhist->second == 0) || (itOPSAhist->second == OPSA_detID) ) boolStoreThisOPSAhist = true;
}
}
if (multimapOfEventIDsForOPSAhistos.find(-1)!=multimapOfEventIDsForOPSAhistos.end()) {
// The user requires to store OPSA timing histograms summed up for all events together
boolStoreThisOPSAhistSUMMED = true;
}
}
if (optHitDetectorMap->empty()) continue;
@ -537,16 +559,20 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
optHitDetectorMapType::const_iterator it2_LAST = optHitDetectorMap->end(); it2_LAST--;
Double_t time = -1000, lastTime = -1000;
G4int OPSA_nPhot = 0;
G4double OPSA_timeFirst = -1000;
G4double OPSA_timeA = -1000;
G4double OPSA_timeB = -1000;
G4double OPSA_timeC = -1000;
G4double OPSA_timeD = -1000;
G4double OPSA_timeE = -1000;
G4double OPSA_timeLast = -1000;
G4double OPSA_timeFirst = -1000000;
G4double OPSA_timeA = -1000000;
G4double OPSA_timeB = -1000000;
G4double OPSA_timeC = -1000000;
G4double OPSA_timeD = -1000000;
G4double OPSA_timeE = -1000000;
G4double OPSA_timeLast = -1000000;
G4double OPSA_CFD_time = -1000000;
G4double OPSA_CFD_ampl = -1000;
G4int iHistNr = -1;
TF1* poiss = NULL;
G4int iHistNrSUM = -1;
// TF1* poiss = NULL;
for (optHitDetectorMapType::const_iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
//delete G4cout<<"KAMIL: Nr of optHitDetectorMap->size()="<<optHitDetectorMap->size()<<G4endl;
OPSA_nPhot++;
lastTime = time;
time = it2->first;
@ -563,62 +589,147 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
it2_START = it2_STOP;
if (OPSA_nPhot >= OPSA_minNrOfDetectedPhotons) { // ignore hits with too low number of detected photons
G4double OPSA_f_nPhot = OPSA_nPhot;
G4int NA = int (OPSA_fracA * OPSA_f_nPhot + 0.5);
G4int NB = int (OPSA_fracB * OPSA_f_nPhot + 0.5);
G4int NC = int (OPSA_fracC * OPSA_f_nPhot + 0.5);
G4int ND = int (OPSA_fracD * OPSA_f_nPhot + 0.5);
G4int NE = int (OPSA_fracE * OPSA_f_nPhot + 0.5);
G4int NA = int (OPSA_fracA * OPSA_f_nPhot + 0.5); if (NA<=0) NA=1;
G4int NB = int (OPSA_fracB * OPSA_f_nPhot + 0.5); if (NB<=0) NB=1;
G4int NC = int (OPSA_fracC * OPSA_f_nPhot + 0.5); if (NC<=0) NC=1;
G4int ND = int (OPSA_fracD * OPSA_f_nPhot + 0.5); if (ND<=0) ND=1;
G4int NE = int (OPSA_fracE * OPSA_f_nPhot + 0.5); if (NE<=0) NE=1;
Int_t nP=0;
// Define OPSA histograms if required for this event
if (boolStoreThisOPSAhist) {
if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) {
iHistNr++;
char nameHist[100]; sprintf(nameHist,"OPSAhist_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
char nameHistTitle[100]; sprintf(nameHistTitle,"OPSAhist_%d_%d_%d;time (ns);Nr of photons",eeeventID,OPSA_detID,iHistNr);
char nameHist[200]; sprintf(nameHist,"OPSAhist_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
char nameHistTitle[200]; sprintf(nameHistTitle,"OPSAhist_%d_%d_%d;time (ns);Nr of photons",eeeventID,OPSA_detID,iHistNr);
OPSAhisto = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
poiss = new TF1("poiss",poissonf,0.,.5,2); // x in [0;300], 2
poiss->SetParameter(0,1);
poiss->SetParameter(1,1);
// poiss = new TF1("poiss",poissonf,0.,.5,2); // x in [0;300], 2
// poiss->SetParameter(0,1);
// poiss->SetParameter(1,1);
if (bool_pulseShapeExists) {
sprintf(nameHist,"OPSAshape_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
sprintf(nameHistTitle,"OPSAshape_%d_%d_%d;time (ns);Pulse signal",eeeventID,OPSA_detID,iHistNr);
OPSAshape = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
sprintf(nameHist,"OPSA_CFD_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
sprintf(nameHistTitle,"OPSA_CFD_%d_%d_%d;time (ns);CFD signal",eeeventID,OPSA_detID,iHistNr);
OPSA_CFD = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
}
}
if (boolStoreThisOPSAhistSUMMED) {
iHistNrSUM++;
char nameHist[200]; sprintf(nameHist,"OPSAhistSUM_%d_%d",OPSA_detID,iHistNrSUM);
char nameHist0[200]; sprintf(nameHist0,"OPSAhistSUM0_%d_%d",OPSA_detID,iHistNrSUM);
if (mapOfOPSAsumHistograms.find(nameHist) != mapOfOPSAsumHistograms.end()) {
OPSAhistoSUM = mapOfOPSAsumHistograms[nameHist];
OPSAhistoSUM0 = mapOfOPSAsum0Histograms[nameHist0];
// G4cout<<" OPSAhistoSUM histogram found:"<<OPSAhistoSUM<<G4endl;
}
else { // create histogram because it does not exist
char nameHistTitle[200]; sprintf(nameHistTitle,"OPSAhistSUM_%d_%d;time (ns);Nr of photons",OPSA_detID,iHistNrSUM);
char nameHistTitle0[200]; sprintf(nameHistTitle0,"OPSAhistSUM0_%d_%d;time (ns);Nr of photons",OPSA_detID,iHistNrSUM);
OPSAhistoSUM = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
OPSAhistoSUM0= new TH1D(nameHist0,nameHistTitle0,OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
mapOfOPSAsumHistograms.insert(std::pair<std::string,TH1D*>(nameHist,OPSAhistoSUM));
mapOfOPSAsum0Histograms.insert(std::pair<std::string,TH1D*>(nameHist0,OPSAhistoSUM0));
// G4cout<<" New OPSAhistoSUM histogram created:"<<OPSAhistoSUM<<G4endl;
}
}
for (optHitDetectorMapType::const_iterator it3 = optHitDetectorSUBmap.begin(); it3 != optHitDetectorSUBmap.end(); it3++) {
nP++;
if (nP==1) OPSA_timeFirst = it3->first;
if (nP==NA) OPSA_timeA = it3->first;
if (nP==NB) OPSA_timeB = it3->first;
if (nP==NC) OPSA_timeC = it3->first;
if (nP==ND) OPSA_timeD = it3->first;
if (nP==NE) OPSA_timeE = it3->first;
if (nP==OPSA_nPhot) OPSA_timeLast = it3->first;
if (boolStoreThisOPSAhist) OPSAhisto->Fill((it3->first)-OPSA_timeFirst+0.00000000001);
G4double timePhot = it3->first;
if (nP==1) OPSA_timeFirst = timePhot;
if (nP==NA) OPSA_timeA = timePhot;
if (nP==NB) OPSA_timeB = timePhot;
if (nP==NC) OPSA_timeC = timePhot;
if (nP==ND) OPSA_timeD = timePhot;
if (nP==NE) OPSA_timeE = timePhot;
if (nP==OPSA_nPhot) OPSA_timeLast = timePhot;
if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) OPSAhisto->Fill(timePhot-OPSA_timeFirst+0.00000000001);
if (boolStoreThisOPSAhistSUMMED) {
OPSAhistoSUM ->Fill(timePhot-OPSA_timeFirst+0.00000000001);
OPSAhistoSUM0->Fill(timePhot+0.00000000001);
}
}
// OPSAhisto->Fit(poiss,"Q");
if (boolStoreThisOPSAhist) OPSAhisto->Write();
// if required, convert the histogram with photon times to histogram of electronic pulse shape
if (bool_pulseShapeExists) {
for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of photon histogram
Double_t photonCount = OPSAhisto ->GetBinContent(iBin);
if (photonCount==0.) break; //nothing to be done for bins without any photons
for (Int_t iBinNew=0; iBinNew<(OPSAhistoNbin-iBin); iBinNew++) { // loop over all bins from the actual one to the end
Int_t iPS = int( (double(iBinNew)+0.5)*OPSAhistoBinWidth1000 );
if (iPS>=iPSmax-1) break; // out of the pulse shape info ==> leave this loop
// OPSAshape->AddBinContent(iBin+iBinNew,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth); // iPS corresponds to time in picoseconds
OPSAshape->Fill((iBin+iBinNew-0.5)*OPSAhistoBinWidth,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth);
}
}
if (boolStoreThisOPSAhist) OPSAshape -> Write();
// Now fill the histogram with the CFD signal
for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of electronic signal histogram
Double_t signalValue = OPSAshape->GetBinContent(iBin);
Double_t xValue = (iBin-0.5)*OPSAhistoBinWidth;
OPSA_CFD->Fill(xValue,signalValue*OPSA_CFD_a1);
Double_t xValueShifted = xValue + OPSA_CFD_delay;
if (xValueShifted>OPSAhistoMax) break; //return if out of range;
OPSA_CFD->Fill(xValueShifted,signalValue);
}
if (boolStoreThisOPSAhist) OPSA_CFD -> Write();
// Find the timeCFD from CFD signal
Double_t oldYvalue;
Double_t yValue=-1000;
G4bool OPSA_CFD_time_found = false;
for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of CFD histogram
Double_t xValue = (iBin-0.5)*OPSAhistoBinWidth;
oldYvalue = yValue;
yValue = OPSA_CFD->GetBinContent(iBin);
if (yValue>OPSA_CFD_ampl) OPSA_CFD_ampl=yValue; // find the CFD signal amplitude;
if (xValue < 0.6*OPSA_CFD_delay) continue; // it is believed the CFD signal can not cross y=0 much before the delay time
if ((yValue >= 0)&&(!OPSA_CFD_time_found)) { // signal just crossed the y=0 axis for the first time ==> record CFD time
OPSA_CFD_time_found = true;
OPSA_CFD_time = xValue - yValue/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
OPSA_CFD_time += OPSA_timeFirst - OPSA_CFD_delay + OPSA_CFD_timeShiftOffset; // correct for
}
}
// G4cout<<" OPSA_CFD_time = "<< OPSA_CFD_time<<" OPSA_CFD_ampl="<<OPSA_CFD_ampl<<" OPSA_timeFirst"<<OPSA_timeFirst<<G4endl;
} // end if (pulseShapeExists)
// Delete the histograms from memory if they were created
if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) {
delete OPSAhisto;
if (bool_pulseShapeExists) {
delete OPSAshape;
delete OPSA_CFD;
}
}
signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_timeFirst,OPSA_timeA,OPSA_timeB,OPSA_timeC,
OPSA_timeD,OPSA_timeE,OPSA_timeLast);
OPSA_timeD,OPSA_timeE,OPSA_timeLast,OPSA_CFD_time,OPSA_CFD_ampl);
OPSA_signal_Map.insert(std::pair<G4int,signalInfo*>(OPSA_nPhot,mySignalInfo) );
if (boolStoreThisOPSAhist) {
OPSAhisto->Fit(poiss,"Q");
OPSAhisto->Write();
// TCanvas* cc = new TCanvas();
// OPSAhisto->Draw();
// cc->Update();
}
}
OPSA_nPhot = 0;
OPSA_timeFirst = -1000;
OPSA_timeA = -1000;
OPSA_timeB = -1000;
OPSA_timeC = -1000;
OPSA_timeD = -1000;
OPSA_timeE = -1000;
OPSA_timeLast = -1000;
OPSA_timeFirst = -1000000;
OPSA_timeA = -1000000;
OPSA_timeB = -1000000;
OPSA_timeC = -1000000;
OPSA_timeD = -1000000;
OPSA_timeE = -1000000;
OPSA_timeLast = -1000000;
OPSA_CFD_time = -1000000;
OPSA_CFD_ampl = -1000;
}
}
}
// Now delete all optHitDetectorMap*
for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
delete (it->second);
}
optHitMap.clear();
// // Now delete all optHitDetectorMap*
// for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
// delete (it->second);
// }
// optHitMap.clear();
// Now store the information of all mySignalInfo* to rootOutputFile
@ -634,4 +745,58 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
OPSA_signal_Map.clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::EndOfRun() {
for (mapOfOPSAsumHistograms_Type::const_iterator it = mapOfOPSAsumHistograms.begin(); it!=mapOfOPSAsumHistograms.end(); it++) {
(it->second) -> Write();
}
for (mapOfOPSAsumHistograms_Type::const_iterator it = mapOfOPSAsum0Histograms.begin(); it!=mapOfOPSAsum0Histograms.end(); it++) {
(it->second) -> Write();
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::ReadInPulseShapeArray(const char* filename) {
bool_pulseShapeExists = true;
std::ifstream file( filename );
// FILE* file = fopen(filename,"r");
if (!(file.is_open())) {
// if (file==NULL) {
G4cout << "musrScintSD::ReadInPulseShapeArray: requested pulse shape data file \""<< filename << "\" not opened (found) !!!"<<G4endl;
musrErrorMessage::GetInstance()->musrError(FATAL,"musrScintSD: pulse shape data file not found !",false);
}
do {
file.ignore(256, '\n');
} while (file.peek() == '%');
G4int i = -1;
G4int idummy;
G4double ampl;
while (!file.eof()) {
// while (!feof(fSteeringFile)) {
// if (file.eof()) break;
file >> idummy >> ampl;
if (file.eof()) break;
i++;
if (i!=idummy) {
G4cout<<"musrScintSD::ReadInPulseShapeArray: pulse shape data file corrupted! i<>idummy (i="<<i<<", idummy="<<idummy<<")"<<G4endl;
musrErrorMessage::GetInstance()->musrError(FATAL,"musrScintSD: pulse shape data file corrupted !",false);
}
if (i>10000) {
G4cout<<"musrScintSD::ReadInPulseShapeArray: pulse shape data file can not handle such big array ! i>10000 (i="<<i<<")"<<G4endl;
musrErrorMessage::GetInstance()->musrError(FATAL,"musrScintSD: pulse shape data file problem !",false);
}
pulseShapeArray[i] = ampl;
}
iPSmax = i+1;
file.close( );
// for (int j=0; j<iPSmax; j++) {G4cout<<j<<" "<< pulseShapeArray[j]<<G4endl;}
// G4cout<<G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

4
src/musrStackingAction.cc
View File

@ -57,8 +57,8 @@ G4ClassificationOfNewTrack musrStackingAction::ClassifyNewTrack(const G4Track *
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrStackingAction::NewStage() {
// G4cout << "Number of optical photons produces in this event : "
// << opticalPhotonCounter << G4endl;
// G4cout << "Number of optical photons produces in this event : "
// << opticalPhotonCounter << G4endl;
musrRootOutput::GetRootInstance()->SetNOptPhot(opticalPhotonCounter);
}