diff --git a/doc/musrSim.pdf b/doc/musrSim.pdf index 8a4c3b9..915785f 100644 Binary files a/doc/musrSim.pdf and b/doc/musrSim.pdf differ diff --git a/doc/musrSim.tex b/doc/musrSim.tex index 5a4f289..926ae3b 100644 --- a/doc/musrSim.tex +++ b/doc/musrSim.tex @@ -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 diff --git a/include/musrScintSD.hh b/include/musrScintSD.hh index 8445b88..05dbbef 100644 --- a/include/musrScintSD.hh +++ b/include/musrScintSD.hh @@ -95,51 +95,65 @@ 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; - musrScintHitsCollection* scintCollection; - G4bool myStoreOnlyEventsWithHits; - G4int myStoreOnlyEventsWithHitInDetID; - G4double mySignalSeparationTime; - G4bool myStoreOnlyTheFirstTimeHit; - G4bool boolIsVvvInfoRequested; - musrRootOutput* myRootOutput; + static musrScintSD* pointer; + musrScintHitsCollection* scintCollection; + G4bool myStoreOnlyEventsWithHits; + G4int myStoreOnlyEventsWithHitInDetID; + G4double mySignalSeparationTime; + G4bool myStoreOnlyTheFirstTimeHit; + G4bool boolIsVvvInfoRequested; + musrRootOutput* myRootOutput; - // for optical photon counting - typedef std::multimap optHitDetectorMapType; - typedef std::map optHitMapType; - optHitMapType optHitMap; - // for optical photon signal analysis (OPSA) - G4int OPSA_minNrOfDetectedPhotons; - G4double OPSA_signalSeparationTime; - G4double OPSA_fracA; - G4double OPSA_fracB; - G4double OPSA_C_threshold; - G4double OPSA_D_threshold; - typedef std::multimap OPSA_signal_MapType; - OPSA_signal_MapType OPSA_signal_Map; + // for optical photon counting + typedef std::multimap optHitDetectorMapType; + typedef std::map optHitMapType; + optHitMapType optHitMap; + // for optical photon signal analysis (OPSA) + G4int OPSA_minNrOfDetectedPhotons; + G4double OPSA_signalSeparationTime; + G4double OPSA_fracA; + G4double OPSA_fracB; + G4double OPSA_C_threshold; + G4double OPSA_D_threshold; + typedef std::multimap OPSA_signal_MapType; + OPSA_signal_MapType OPSA_signal_Map; - G4bool bool_multimapOfEventIDsForOPSAhistosEXISTS; - G4bool bool_StoreThisOPSAhistSUMMED; - G4bool bool_StoreThisOPSAhistALL; - typedef std::multimap multimapOfEventIDsForOPSAhistos_Type; - multimapOfEventIDsForOPSAhistos_Type multimapOfEventIDsForOPSAhistos; - TH1D* OPSAhisto; - Int_t OPSAhistoNbin; - Double_t OPSAhistoMin, OPSAhistoMax, OPSAhistoBinWidth, OPSAhistoBinWidth1000; - typedef std::map 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; - + G4bool bool_multimapOfEventIDsForOPSAhistosEXISTS; + G4bool bool_StoreThisOPSAhistSUMMED; + G4bool bool_StoreThisOPSAhistALL; + typedef std::multimap multimapOfEventIDsForOPSAhistos_Type; + multimapOfEventIDsForOPSAhistos_Type multimapOfEventIDsForOPSAhistos; + TH1D* OPSAhisto; + Int_t OPSAhistoNbin; + Double_t OPSAhistoMin, OPSAhistoMax, OPSAhistoBinWidth, OPSAhistoBinWidth1000; + typedef std::map 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; + 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...... diff --git a/musrSimAna/musrAnalysis.cxx b/musrSimAna/musrAnalysis.cxx index 9970a10..aec0808 100644 --- a/musrSimAna/musrAnalysis.cxx +++ b/musrSimAna/musrAnalysis.cxx @@ -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="< S T O P "< S T O P "<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()"<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); } } diff --git a/musrSimAna/musrCounter.hh b/musrSimAna/musrCounter.hh index 89e4a84..058622f 100644 --- a/musrSimAna/musrCounter.hh +++ b/musrSimAna/musrCounter.hh @@ -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 diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 5a1bafe..381b795 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -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 \"" + <insert(std::pair(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="<begin(); it2 != optHitDetectorMap->end(); it2++ ) { + if ((it2->second)==APDcellID) { + APDcellAlreadyFired = true; + // G4cout<<"HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH already fired cell ="<insert(std::pair(tmpTime,APDcellID)); // G4cout<<" tmpTime="<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 = "<GetPostStepPoint(); + G4ThreeVector postStepPoint = postPoint->GetPosition(); + const G4AffineTransform* trans = postPoint->GetTouchable()->GetHistory()->GetPtrTopTransform(); + G4ThreeVector postStepPointLocal = trans->TransformPoint(postStepPoint); + // G4cout<<"musrScintSD::FindAPDcellID: postStepPointLocal="<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 = "<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...... diff --git a/src/musrStackingAction.cc b/src/musrStackingAction.cc index 252d34a..92187a3 100644 --- a/src/musrStackingAction.cc +++ b/src/musrStackingAction.cc @@ -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; }