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

29.3.2011 Kamil Sedlak

Browse Source 
1) Small change in optical photon simulation
2) Some other small changes, probably not very important
This commit is contained in:
sedlak 2011-03-29 15:34:13 +00:00
parent 2a1c05b09b
commit 9bc0d53074
9 changed files with 199 additions and 62 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.

2
doc/musrSim.tex
View File

@ -356,7 +356,7 @@ is just a dead material influencing the penetrating particles but not detecting
some volume or whether it is read in from an external file as a field-map.
\item \emph{X}, \emph{Y}, \emph{Z} (floats) -- position of the centre of the field
in the {\bf global} coordinate system. IMPORTANT: For some technical internal
Geant4 reasons, this \bf{POSITION HAS TO LAY WITHIN THE \emph{logicalVolume}!}
Geant4 reasons, this {\bf POSITION HAS TO LAY WITHIN THE \emph{logicalVolume}!}
Note that the logical volume may be positioned somewhere deep in a volume
structure, not directly within the ``World'' volume, and
therefore the (local) coordinates in the definition of the logical volume

16
include/musrScintSD.hh
View File

@ -95,6 +95,16 @@ class musrScintSD : public G4VSensitiveDetector
void ProcessOpticalPhoton(G4Step*);
void EndOfEvent_OptiacalPhotons();
void ReadInPulseShapeArray(const char* filename);
G4int FindAPDcellID(G4Step* aStep);
void SetAPDcellSizes(G4int nx, G4int ny, G4int nz, G4double half_ax, G4double half_ay, G4double half_az) {
APDcellsEffectRequested = true;
APDcell_nx = nx; APDcell_ny = ny; APDcell_nz = nz;
APDcell_ax = 2*half_ax/ nx; APDcell_ay = 2*half_ay/ ny; APDcell_az = 2*half_az/ nz;
}
void SetAPDcellsTimeVariationSigma(G4double sigma) {
if (sigma!=0) APDcellsTimeVariationRequested = true;
APDcellsTimeVariationSigma=sigma;
}
private:
static musrScintSD* pointer;
@ -139,7 +149,11 @@ class musrScintSD : public G4VSensitiveDetector
G4double pulseShapeArray[10001];
TH1D* OPSA_CFD;
G4double OPSA_CFD_a1,OPSA_CFD_delay, OPSA_CFD_timeShiftOffset;
G4bool APDcellsEffectRequested; // simulate effects of finite cell number
G4int APDcell_nx, APDcell_ny, APDcell_nz; // number of cells in APD along x, y and z directions
G4double APDcell_ax, APDcell_ay, APDcell_az; // dimensions of APD cells
G4bool APDcellsTimeVariationRequested; // simmulate effect of detection time variations between differenct cells
G4double APDcellsTimeVariationSigma; // sigma of the detection time variations between differenct cells
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

54
musrSimAna/musrAnalysis.cxx
View File

@ -455,6 +455,39 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
nscan = sscanf(pch,"%d",&N1);
} while (nscan==1);
}
else if (strcmp(tmpString0,"setSpecialAnticoincidenceTimeWindow")==0) {
int channelNumber;
double min, max;
char unit[100]="Undefined";
Long64_t lmin, lmax;
sscanf(&line[0],"%*s %d %lf %lf %s",&channelNumber,&min,&max,unit);
if ((strcmp(unit,"bin")==0)||(strcmp(unit,"bins")==0)) {lmin = Long64_t (min); lmax = Long64_t (max);}
else if ((strcmp(unit,"ns") ==0)||(strcmp(unit,"nanosecond") ==0)) {
lmin = Long64_t (min*nanosecond/tdcresolution);
lmax = Long64_t (max*nanosecond/tdcresolution);
}
else if ((strcmp(unit,"mus") ==0)||(strcmp(unit,"microsecond") ==0)) {
lmin = Long64_t (min*microsecond/tdcresolution);
lmax = Long64_t (max*microsecond/tdcresolution);
}
else {
std::cout<<" ERROR! musrAnalysis: unknown unit in setSpecialAnticoincidenceTimeWindow !!! (unit="<<unit<<")"<<std::endl;
std::cout<<" ==> S T O P "<<std::endl;
exit(1);
}
std::cout<<" lmin="<<lmin<<" bins"<<std::endl;
std::cout<<" lmax="<<lmax<<" bins"<<std::endl;
if (allCounterMap.find(channelNumber)==allCounterMap.end()) {
std::cout<<" ERROR! musrAnalysis: unknown counter ("<<channelNumber<<") in setSpecialAnticoincidenceTimeWindow"<<std::endl;
std::cout<<" ==> S T O P "<<std::endl;
exit(1);
}
musrCounter* counter = (allCounterMap.find(channelNumber))->second;
counter->SetAntiCoincidenceTimeWindowMin(lmin);
counter->SetAntiCoincidenceTimeWindowMax(lmax);
}
else if (strcmp(tmpString0,"artificiallyChangeMuDecayTime")==0) {
float min, max, mmmin, mmmax;
sscanf(&line[0],"%*s %g %g %g %g %g %g %g",&min,&max,&mmmin,&mmmax);
@ -759,8 +792,11 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
}
Long64_t numberOfMuonCandidates = mCounter->GetNumberOfMuonCandidates();
Long64_t numberOfMuonCandidatesAfterVK = mCounter->GetNumberOfMuonCandidatesAfterVK();
Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000;
hInfo->Fill(5,(Double_t) numberOfMuonCandidates); // number of "triggers"
hInfo->Fill(4,(Double_t) numberOfMuonCandidatesAfterVK); // number of "triggers"
hInfo->Fill(3,durationOfExperiment);
//==============================
// Write out the histograms
std::cout<<"musrAnalysis::SaveHistograms()"<<std::endl;
@ -780,14 +816,18 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
std::cout<<"==================================================================="<<std::endl;
std::cout<<"Number of raw trigger counts (ignoring pileup gate, vetos and coincidences): "<<numberOfMuonCandidates<<std::endl;
std::cout<<"Number of trigger counts (after vetos and coincidences but not pile-up rejected): "<<numberOfMuonCandidatesAfterVK<<std::endl;
std::cout<<"Number of triggered events (i.e. only good \"muons\"): "<<numberOfGoodMuons<<std::endl;
Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000;
std::cout<<"Duration of the \"experiment\": "<<durationOfExperiment<<" second"<<std::endl;
// std::cout<<" (numberOfRewinds="<<numberOfRewinds<<")"<<std::endl;
std::cout<<"In this run, MUONRATEFACTOR was set to "<<muonRateFactor<<", and the raw trigger rate was "
<<numberOfMuonCandidates/durationOfExperiment<<" muons/second"<<std::endl;
std::cout<<" To get event rate of 30000 events/second, set MUONRATEFACTOR to "
<< muonRateFactor*( (numberOfMuonCandidates/durationOfExperiment)/30000 )<<std::endl;
std::cout<<"MUONRATEFACTOR was set to "<<muonRateFactor<<std::endl;
std::cout<<" Raw trigger rate = "<<numberOfMuonCandidates/durationOfExperiment<<" muons/second";
std::cout<<" (to get 30000 events/second, set MUONRATEFACTOR = "
<< muonRateFactor*( (numberOfMuonCandidates/durationOfExperiment)/30000 )<<")"<<std::endl;
std::cout<<" Trigger rate after V & K = "<<numberOfMuonCandidatesAfterVK/durationOfExperiment<<" muons/second";
std::cout<<" (to get 30000 events/second, set MUONRATEFACTOR = "
<< muonRateFactor*( (numberOfMuonCandidatesAfterVK/durationOfExperiment)/30000 )<<")"<<std::endl;
std::cout<<"========================== E N D =================================="<<std::endl;
}
@ -985,7 +1025,7 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
goodEvent_R_det = goodEvent_det && ( (find(R_posCounterList.begin(), R_posCounterList.end(), idetP_ID)) != R_posCounterList.end() );
// std::cout<<"goodEvent_F_det="<<goodEvent_F_det<<std::endl;
if (pileupEvent&&goodEvent_F_det) {
std::cout<<" DEBUG: Pileup Event: eventID = "<<eventID<<" pileup_eventID = "<<pileup_eventID<<" det_time10 = "<<det_time10<<std::endl;
// std::cout<<" DEBUG: Pileup Event: eventID = "<<eventID<<" pileup_eventID = "<<pileup_eventID<<" det_time10 = "<<det_time10<<std::endl;
// debugEventMap.insert(std::pair<int,int>(eventID,10));
}

15
musrSimAna/musrCounter.cxx
View File

@ -31,6 +31,7 @@ musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TY
strcpy(TDC_histoNameAdd,"Unset");
doubleHitN=0;
numberOfMuonCandidates=0;
numberOfMuonCandidatesAfterVK=0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -182,6 +183,9 @@ Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t&
for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) {
if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') ) goto endOfThisHit; // coincidence with veto found ==> skip hit
}
numberOfMuonCandidatesAfterVK++;
// Check coincidences with other hits in the M counter
// it is expected that there is only one M counter
if (this->IsInCoincidence(timeBinOfCount_tmp,'M',true) ) { // coincidence with another M-counter hit ==> skip hit
@ -300,10 +304,13 @@ void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherC
//================================================================
void musrCounter::SetCoincidenceTimeWindowOfAllVetoDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max) {
for (counterMapType::const_iterator it = vetoCounterMap.begin(); it!=vetoCounterMap.end(); ++it) {
Long64_t maxCoinc_AlreadySet = ((it->second)->GetMaxCoincidenceTimeWindow());
if (maxCoinc < maxCoinc_AlreadySet) (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
(it->second)->SetAntiCoincidenceTimeWindow(min,max);
musrCounter* counter = it->second;
Long64_t maxCoinc_AlreadySet = counter->GetMaxCoincidenceTimeWindow();
Long64_t min_AlreadySet = counter->GetAntiCoincidenceTimeWindowMin();
Long64_t max_AlreadySet = counter->GetAntiCoincidenceTimeWindowMax();
if (maxCoinc < maxCoinc_AlreadySet) counter->SetMaxCoincidenceTimeWindow(maxCoinc);
if (min < min_AlreadySet) counter->SetAntiCoincidenceTimeWindowMin(min);
if (max > max_AlreadySet) counter->SetAntiCoincidenceTimeWindowMax(max);
}
}

7
musrSimAna/musrCounter.hh
View File

@ -48,7 +48,10 @@ class musrCounter {
void SetCoincidenceTimeWindowOfAllVetoDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max);
void SetCoincidenceTimeWindow_M(Long64_t min, Long64_t max) {coincidenceTimeWindowMin_M=min; coincidenceTimeWindowMax_M=max;}
void SetCoincidenceTimeWindow_P(Long64_t min, Long64_t max) {coincidenceTimeWindowMin_P=min; coincidenceTimeWindowMax_P=max;}
void SetAntiCoincidenceTimeWindow(Long64_t min, Long64_t max) {antiCoincidenceTimeWindowMin=min; antiCoincidenceTimeWindowMax=max;}
void SetAntiCoincidenceTimeWindowMin(Long64_t min) {antiCoincidenceTimeWindowMin=min;}
void SetAntiCoincidenceTimeWindowMax(Long64_t max) {antiCoincidenceTimeWindowMax=max;}
Long64_t GetAntiCoincidenceTimeWindowMin() {return antiCoincidenceTimeWindowMin;}
Long64_t GetAntiCoincidenceTimeWindowMax() {return antiCoincidenceTimeWindowMax;}
void SetTDChistogram(char hName[200],int t0,int t1,int t2,int hNr,char hNameAdd[200]);
void FillTDChistogram(Double_t variable, Double_t vaha);
void DrawTDChistogram();
@ -60,6 +63,7 @@ class musrCounter {
Int_t GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound);
void myPrintThisCounter(Int_t evtID, Int_t detail=2);
Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;}
Long64_t GetNumberOfMuonCandidatesAfterVK(){return numberOfMuonCandidatesAfterVK;}
private:
@ -86,6 +90,7 @@ class musrCounter {
Int_t doubleHitN;
Long64_t numberOfMuonCandidates;
Long64_t numberOfMuonCandidatesAfterVK;
};
#endif

16
src/musrDetectorConstruction.cc
View File

@ -736,6 +736,22 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf",&a1,&delay,&timeShiftOffset);
myMusrScintSD -> Set_OPSA_CFD(a1,delay,timeShiftOffset);
}
else if (strcmp(varName,"APDcells")==0) {
int nx, ny, nz;
double halfLenght_x, halfLenght_y, halfLenght_z;
sscanf(&line[0],"%*s %*s %*s %d %d %d %lf %lf %lf",&nx,&ny,&nz,&halfLenght_x,&halfLenght_y,&halfLenght_z);
myMusrScintSD -> SetAPDcellSizes(nx, ny, nz, halfLenght_x, halfLenght_y, halfLenght_z);
}
else if (strcmp(varName,"SetAPDcellsTimeVariationSigma")==0) {
double sigma;
sscanf(&line[0],"%*s %*s %*s %lf",&sigma);
myMusrScintSD ->SetAPDcellsTimeVariationSigma(sigma*nanosecond);
}
else {
G4cout<<"musrDetectorConstruction.cc: ERROR: Unknown parameterName \""
<<varName<<"\" ."<<G4endl;
ReportGeometryProblem(line);
}
}
}

68
src/musrScintSD.cc
View File

@ -84,6 +84,10 @@ musrScintSD::musrScintSD(G4String name)
OPSA_CFD_a1 = -0.2;
OPSA_CFD_delay = 2.;
OPSA_CFD_timeShiftOffset = 0.;
APDcell_nx =1; APDcell_ny=10; APDcell_nz=10;
APDcell_ax =0.3; APDcell_ay=0.3; APDcell_az=0.3;
APDcellsEffectRequested=false;
APDcellsTimeVariationRequested=false;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -259,7 +263,25 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
// case Detection: //Note, this assumes that the volume causing detection
//is the photocathode because it is the only one with
//non-zero efficiency
optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,boundaryStatus));
// If the simulation of the effect of the finite number of cells in APD is requested, find out
// whether a given cell already fired - if so, do not store the photon.
G4bool APDcellAlreadyFired = false;
G4int APDcellID = 0;
if (APDcellsEffectRequested) {
APDcellID = FindAPDcellID(aStep);
// G4cout<<"Cell ID="<<APDcellID<<G4endl;
for (optHitDetectorMapType::const_iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
if ((it2->second)==APDcellID) {
APDcellAlreadyFired = true;
// G4cout<<"HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH already fired cell ="<<APDcellID<<G4endl;
break;
}
}
}
if (APDcellsTimeVariationRequested) tmpTime += G4RandGauss::shoot(0,APDcellsTimeVariationSigma);
if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
// G4cout<<" tmpTime="<<tmpTime<<G4endl;
// G4cout<<" Detection"<<G4endl;
// break;
@ -636,6 +658,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
for (optHitDetectorMapType::const_iterator it3 = optHitDetectorSUBmap.begin(); it3 != optHitDetectorSUBmap.end(); it3++) {
nP++;
G4double timePhot = it3->first;
// if (APDcellsTimeVariationRequested) timePhot += G4RandGauss::shoot(0,APDcellsTimeVariationSigma); // Shifted above
if (nP==1) OPSA_timeFirst = timePhot;
if (nP==NA) OPSA_timeA = timePhot;
if (nP==NB) OPSA_timeB = timePhot;
@ -694,17 +717,24 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
// 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;
// G4bool OPSA_CFD_time_found = false;
OPSA_CFD_ampl = OPSA_CFD->GetMaximum();
int binmax = OPSA_CFD->GetMaximumBin();
int binmim = OPSA_CFD->GetMinimumBin();
for (Int_t iBin=binmim; iBin<=binmax; iBin++) { // loop over bins between min and max of CFD histogram
// Double_t xValue = (iBin-0.5)*OPSAhistoBinWidth;
Double_t xValue = OPSA_CFD->GetXaxis()->GetBinCenter(iBin);
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;
// 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
if (yValue >= 0) { // signal just crossed the y=0 axis
// 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<<G4endl;
break;
}
}
@ -839,3 +869,25 @@ void musrScintSD::ReadInPulseShapeArray(const char* filename) {
// G4cout<<G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4int musrScintSD::FindAPDcellID(G4Step* aStep) {
// TEST:
G4StepPoint* postPoint = aStep->GetPostStepPoint();
G4ThreeVector postStepPoint = postPoint->GetPosition();
const G4AffineTransform* trans = postPoint->GetTouchable()->GetHistory()->GetPtrTopTransform();
G4ThreeVector postStepPointLocal = trans->TransformPoint(postStepPoint);
// G4cout<<"musrScintSD::FindAPDcellID: postStepPointLocal="<<postStepPointLocal<<G4endl;
G4int ix = (APDcell_nx>1) ? int ((postStepPointLocal.x() / APDcell_ax + APDcell_nx/2.)) : 0;
G4int iy = (APDcell_ny>1) ? int ((postStepPointLocal.y() / APDcell_ay + APDcell_ny/2.)) : 0;
G4int iz = (APDcell_nz>1) ? int ((postStepPointLocal.z() / APDcell_az + APDcell_nz/2.)) : 0;
G4int in = ix * APDcell_ny * APDcell_nz + iy * APDcell_nz + iz;
// G4cout<<"musrScintSD::FindAPDcellID: in = "<<in<<", ix = "<<ix<<", iy = "<<iy<<", iz = "<<iz<<G4endl;
if ((ix<0)||(ix>APDcell_nx)) musrErrorMessage::GetInstance()->musrError(FATAL,
"musrScintSD::FindAPDcellID: APD cell out of range (ix). Wrong dimensions of APD in \"/musr/command OPSA APDcells\" command?",false);
if ((iy<0)||(iy>APDcell_ny)) musrErrorMessage::GetInstance()->musrError(FATAL,
"musrScintSD::FindAPDcellID: APD cell out of range (iy). Wrong dimensions of APD in \"/musr/command OPSA APDcells\" command?",false);
if ((iz<0)||(iz>APDcell_nz)) musrErrorMessage::GetInstance()->musrError(FATAL,
"musrScintSD::FindAPDcellID: APD cell out of range (iz). Wrong dimensions of APD in \"/musr/command OPSA APDcells\" command?",false);
return in;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

3
src/musrStackingAction.cc
View File

@ -50,6 +50,9 @@ musrStackingAction::~musrStackingAction()
G4ClassificationOfNewTrack musrStackingAction::ClassifyNewTrack(const G4Track * aTrack) {
if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) {
opticalPhotonCounter++;
return fWaiting; // Optical photons will be treated after all other particles in order
// to have the reproducibility of the other particle tracks if a parameter
// for the optical photons is changed in the next run.
}
return fUrgent;
}