From f93a3c070bcdf3a879c761246c9c064b04468cd3 Mon Sep 17 00:00:00 2001 From: Kamil Sedlak Date: Fri, 28 Oct 2011 14:04:11 +0000 Subject: [PATCH] 28.10.2011 Kamil Sedlak 1) musrSimAna significantly rewritten, but still needs to be checked. 2) changed way how some optical photon properties are simulated, + added cross-talk effects (if requested). --- include/musrScintSD.hh | 12 ++- musrSimAna/musrAnalysis.cxx | 165 +++++++++++++++++--------------- musrSimAna/musrAnalysis.hh | 29 ++++-- musrSimAna/musrCounter.cxx | 130 +++++++++++++++---------- musrSimAna/musrCounter.hh | 20 +++- src/musrDetectorConstruction.cc | 5 + src/musrScintSD.cc | 110 +++++++++++++++------ 7 files changed, 306 insertions(+), 165 deletions(-) diff --git a/include/musrScintSD.hh b/include/musrScintSD.hh index 38ed1b0..b652b58 100644 --- a/include/musrScintSD.hh +++ b/include/musrScintSD.hh @@ -105,7 +105,7 @@ class musrScintSD : public G4VSensitiveDetector OPSAhistoNbin=nBins; OPSAhistoMin=min; OPSAhistoMax=max; OPSAhistoBinWidth=(max-min)/nBins; OPSAhistoBinWidth1000=OPSAhistoBinWidth*1000; } - void ProcessOpticalPhoton(G4Step* aStep, G4double APDcellsTimeVariation); + void ProcessOpticalPhoton(G4Step* aStep); void EndOfEvent_OptiacalPhotons(); void ReadInPulseShapeArray(const char* filename); void FindCFDtime(G4double& OPSA_CFD_time, G4double& OPSA_CFD_ampl, G4double timeOfFirstPhoton); @@ -119,6 +119,10 @@ class musrScintSD : public G4VSensitiveDetector if (sigma!=0) APDcellsTimeVariationRequested = true; APDcellsTimeVariationSigma=sigma; } + void SetAPDcrossTalk(G4double crosstalkProb) { + if (crosstalkProb>0) APDcrossTalkRequested = true; + APDcrossTalkProb = crosstalkProb; + } private: static musrScintSD* pointer; @@ -168,6 +172,12 @@ class musrScintSD : public G4VSensitiveDetector 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 + G4bool APDcrossTalkRequested; + G4double APDcrossTalkProb; + + void FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time); + + }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/musrSimAna/musrAnalysis.cxx b/musrSimAna/musrAnalysis.cxx index 7efc8eb..eddf752 100644 --- a/musrSimAna/musrAnalysis.cxx +++ b/musrSimAna/musrAnalysis.cxx @@ -100,7 +100,7 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) { pcoincwin = 40; vcoincwin = 80; muonRateFactor = 1.; - rewindTimeBins = 1000000000; + // rewindTimeBins = 1000000000; dataWindowMin = -0.2; dataWindowMax = 10.0; pileupWindowMin = -10.5; @@ -176,7 +176,9 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) { } else if (strncmp(tmpString0,"REWINDTIMEBINS=",strlen("REWINDTIMEBINS="))==0) { sscanf(&line[strlen("REWINDTIMEBINS=")],"%d",&tmpivar); - rewindTimeBins = tmpivar; std::cout<<" JUHELAK: rewindTimeBins="<(ieventToDebug_tmp,iLevelToDebug_tmp)); } + else if (strcmp(tmpString0,"WRITE_OUT_DUMP_FILE")==0) { + // int clock_channelID_tmp; + // int clock_interval_tmp; + // sscanf(&line[0],"%*s %s %d %d",tmpString1,&clock_channelID_tmp,&clock_interval_tmp); + sscanf(&line[0],"%*s %s",tmpString1); + // clock_channelID = clock_channelID_tmp; + // clock_interval = clock_interval_tmp; + musrCounter::bool_WriteDataToDumpFile = true; + // std::cout<<"tmpString1="< STOP !!!"< S T O P"<(tmpChannel,iCoinc)); // tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc; } while (jjsecond)->GetCounterType(); if (DetectorType=='M') { (it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin); - (it->second)->SetCoincidenceTimeWindow_M(pileupWindowBinMin,pileupWindowBinMax); + // (it->second)->SetCoincidenceTimeWindow_M(pileupWindowBinMin,pileupWindowBinMax); (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('M',-mcoincwin+pileupWindowBinMin,-mcoincwin,mcoincwin); (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin); } else if (DetectorType=='P') { (it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin); - (it->second)->SetCoincidenceTimeWindow_P(dataWindowBinMin,dataWindowBinMax); + // (it->second)->SetCoincidenceTimeWindow_P(dataWindowBinMin,dataWindowBinMax); (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('P',-pcoincwin+dataWindowBinMin,-pcoincwin,pcoincwin); (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin); } @@ -835,7 +849,7 @@ void musrAnalysis::CreateHistograms() { hInfo->Fill(1, fieldValue); hInfo->Fill(6, runID); hInfo->Fill(7, hGeantParameters->GetBinContent(7)); - std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<RewindHitsInCounter(timeBinsToRewind); + // (*it).second->RewindHitsInCounter(timeBinsToRewind); + (*it).second->RewindHitsInCounter(); } } @@ -1019,7 +1046,8 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) { // std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM); // mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM); - mCounterHitExistsForThisEventID = MuonCounterHit(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep); + // mCounterHitExistsForThisEventID = MuonCounterHit(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep); + mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep); timeBinOfThePreviouslyProcessedHit = timeBin0; //___________________________________________________________ @@ -1225,7 +1253,8 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) { // std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM); // mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM); - mCounterHitExistsForThisEventID = MuonCounterHit(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep); + // mCounterHitExistsForThisEventID = MuonCounterHit(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep); + mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep); timeBinOfThePreviouslyProcessedHit = timeBin0; // if (mCounterHitExistsForThisEventID) std::cout<<" YYYYYYYYYYYYYYYYYYY check this : LOOOPING AGAIN"<second)->GetNextGoodPositron(evID,dataBinMinimum,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); - // std::cout<<"Debug 40 positronQuality="<0) { - // std::cout<<"Debug 50"< double hit ==> through away this event. - // std::cout<<"Debug 60"<2) {std::cout<<"DEBUGEVENT:"<1) goodPositronFound = true; - } - // std::cout<<"Debug 80"<0); - // std::cout<<"Debug 90"<second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); + if (positronQuality==3) return false; // double hit was found in the same counter + if (positronQuality==2) { + if (goodPositronFound) return false; // double hit was found in a different counter + goodPositronFound = true; + } } + + + // Long64_t dataBinMinimum = dataBinMin; + // do { + // // std::cout<<"Debug 30"<second)->GetNextGoodPositron(evID,dataBinMinimum,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); + // // std::cout<<"Debug 40 positronQuality="<0) { + // // std::cout<<"Debug 50"< double hit ==> through away this event. + // // std::cout<<"Debug 60"<2) {std::cout<<"DEBUGEVENT:"<1) goodPositronFound = true; + // } + // // std::cout<<"Debug 80"<0); + // // std::cout<<"Debug 90"< +#include #include #include #include @@ -293,12 +295,13 @@ public : virtual void SaveHistograms(char* runChar,char* v1190FileName); virtual void RemoveOldHitsFromCounters(Long64_t timeBinLimit); // virtual void RewindAllTimeInfo(Double_t timeToRewind); - virtual void RewindAllTimeInfo(Long64_t timeBinsToRewind); + // virtual void RewindAllTimeInfo(Long64_t timeBinsToRewind); + virtual void RewindAllTimeInfo(); virtual void PrintHitsInAllCounters(); virtual void InitialiseEvent(); virtual Double_t PreprocessEvent(Long64_t iEn); virtual Bool_t PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep); - virtual Bool_t MuonCounterHit(Int_t evID, Long64_t timeBinMin, Long64_t& timeBin0, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep); + // virtual Bool_t MuonCounterHit(Int_t evID, Long64_t timeBinMin, Long64_t& timeBin0, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep); void CopySubstring(char* inputChar,int iStart,int iEnd,char* outputChar); void MyPrintTree(); void MyPrintConditions(); @@ -328,17 +331,12 @@ public : Int_t pcoincwin; Int_t vcoincwin; Double_t muonRateFactor; - Long64_t rewindTimeBins; Double_t dataWindowMin; Double_t dataWindowMax; Double_t pileupWindowMin; Double_t pileupWindowMax; Double_t promptPeakWindowMin; Double_t promptPeakWindowMax; - Long64_t pileupWindowBinMin; - Long64_t pileupWindowBinMax; - Long64_t dataWindowBinMin; - Long64_t dataWindowBinMax; Int_t overallBinDelay; Bool_t boolInfinitelyLowMuonRate; @@ -378,6 +376,13 @@ public : public: static const Int_t nrConditions = 31; Bool_t condition[nrConditions]; + static Long64_t rewindTimeBins; + // static Int_t clock_channelID; + // static Long64_t clock_interval; + static Long64_t pileupWindowBinMin; + static Long64_t pileupWindowBinMax; + static Long64_t dataWindowBinMin; + static Long64_t dataWindowBinMax; private: // HISTOGRAMS: @@ -417,6 +422,16 @@ private: #endif #ifdef musrAnalysis_cxx + +Long64_t musrAnalysis::rewindTimeBins = 1000000000; +Long64_t musrAnalysis::pileupWindowBinMin; +Long64_t musrAnalysis::pileupWindowBinMax; +Long64_t musrAnalysis::dataWindowBinMin; +Long64_t musrAnalysis::dataWindowBinMax; + +//Long64_t musrAnalysis::clock_interval = 512000; +//Int_t musrAnalysis::clock_channelID = 31; + musrAnalysis::musrAnalysis(TTree *tree) { variableMap["muDecayPosX"]=&muDecayPosX; diff --git a/musrSimAna/musrCounter.cxx b/musrSimAna/musrCounter.cxx index 677d695..91dcd9b 100644 --- a/musrSimAna/musrCounter.cxx +++ b/musrSimAna/musrCounter.cxx @@ -1,6 +1,7 @@ //#include #include "musrCounter.hh" #include "TCanvas.h" +#include "musrAnalysis.hh" typedef std::map debugEventMapType; debugEventMapType debugEventMap; @@ -8,6 +9,10 @@ Bool_t bool_debugingRequired; Bool_t musrCounter::bool_ignoreUnperfectMuons = true; Bool_t musrCounter::bool_ignoreUnperfectPositrons = true; +Bool_t musrCounter::bool_WriteDataToDumpFile = false; +//Long64_t musrCounter::previousClock = -1; +//Long64_t musrCounter::CLOCK_INTERVAL = 512000; +ofstream musrCounter::dumpFile; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -68,7 +73,7 @@ void musrCounter::DrawTDChistogram() { } //================================================================ -void musrCounter::FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID){ +void musrCounter::FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID, Int_t eventNum){ //cDEL std::cout<<"FillHitInCounter: timeBin="<first; - // int tempEvnr= it->second; - // hitMap_TMP.insert( std::pair(tempBinT-timeBinsToRewind,tempEvnr) ); hitInfo* tempEvnr= it->second; - tempEvnr->RewindTimeBin2(timeBinsToRewind); - hitMap_TMP.insert( std::pair(tempBinT-timeBinsToRewind,tempEvnr) ); + tempEvnr->RewindTimeBin2(musrAnalysis::rewindTimeBins); + hitMap_TMP.insert( std::pair(tempBinT-musrAnalysis::rewindTimeBins,tempEvnr) ); } hitMap.swap(hitMap_TMP); } //================================================================ -Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, char motherCounter, Bool_t ignoreHitsAtBinZero, Long64_t timeBinMinimum, Long64_t timeBinMaximum){ +Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, char motherCounter){ // timeBin ... time bin, at which the coincidence is searched - // counterTimeShiftOfRequestingCounter ... time shift (in bin units) of the counter, for which the coincidence is searched - // ignoreHitsAtBinZero ... if "true", hits at timeBin will be ignored (needed for searching of coincidence of M counter - // with other M counters or P counters with other P counters) - // "false" should be used for coincidence detectors and vetos. + if (hitMap.empty()) return false; - // Long64_t timeBinToTest = timeBin; Long64_t timeBinMin; Long64_t timeBinMax; // If timeBinMinimum and timeBinMaximum are not specified, use internal time window of the detector (koincidence or veto detectors). // Otherwise use timeBinMinimum and timeBinMaximum (e.g.coincidence of a positron counter with other positron counters). - if (timeBinMinimum!=-123456789) timeBinMin = timeBin + timeBinMinimum; // time window requested through "timeBinMinimum" - else if (counterType == 'V') timeBinMin = timeBin + antiCoincidenceTimeWindowMin; // this is veto detector - else if (motherCounter=='M') timeBinMin = timeBin + coincidenceTimeWindowMin_M; // this is coinc. detector (or M counter) connected to M + if (counterType == 'V') timeBinMin = timeBin + antiCoincidenceTimeWindowMin; // this is veto detector + else if (motherCounter=='M') timeBinMin = timeBin + coincidenceTimeWindowMin_M; // this is coinc. detector connected to M else if (motherCounter=='P') timeBinMin = timeBin + coincidenceTimeWindowMin_P; // this is coinc. detector connected to P - if (timeBinMaximum!=-123456789) timeBinMax = timeBin + timeBinMaximum; // time window requested through "timeBinMinimum" - else if (counterType == 'V') timeBinMax = timeBin + antiCoincidenceTimeWindowMax; // this is veto detector - else if (motherCounter=='M') timeBinMax = timeBin + coincidenceTimeWindowMax_M; // this is coinc. detector (or M counter) connected to M + if (counterType == 'V') timeBinMax = timeBin + antiCoincidenceTimeWindowMax; // this is veto detector + else if (motherCounter=='M') timeBinMax = timeBin + coincidenceTimeWindowMax_M; // this is coinc. detector connected to M else if (motherCounter=='P') timeBinMax = timeBin + coincidenceTimeWindowMax_P; // this is coinc. detector connected to P - // Long64_t timeBinMin = (timeBinMinimum==-123456789) ? timeBin + coincidenceTimeWindowMin : timeBin + timeBinMinimum; - // Long64_t timeBinMax = (timeBinMaximum==-123456789) ? timeBin + coincidenceTimeWindowMax : timeBin + timeBinMaximum; for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { Long64_t timeBinOfCount_tmp = it->first; if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) { - if ((timeBin!=timeBinOfCount_tmp)||(!ignoreHitsAtBinZero)) { - return true; - } + // if ((timeBin!=timeBinOfCount_tmp)||(!ignoreHitsAtBinZero)) { + return true; + // } } else if (timeBinOfCount_tmp > timeBinMax) return false; } @@ -175,50 +176,55 @@ Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& Int_t eventNumber = (it->second)->eventIDnumber; if (eventNumber!=evtID) continue; // This trigger hit does not correspond to the currently processed event - // ==> skip it, becuase it was already proceesed or will be processed in future + // ==> skip it, because it was already proceesed or will be processed in future numberOfMuonCandidates++; - // std::cout<<"*** "<second)->eventEntry; idet = (it->second)->det_i; idetID = (it->second)->det_id; idetEdep = (it->second)->det_edep; - // timeBinOfNextHit = timeBinOfCount_tmp; return true; - - endOfThisHit: - continue; } return false; } //================================================================ -Bool_t musrCounter::CheckForPileupMuons(Long64_t timeBin0, Int_t kEntry_NN) { +Bool_t musrCounter::CheckForPileupMuons(Long64_t timeBin0) { // Check for pileup muons. If double hit in M-counter is found, return true. - Long64_t timeBinMinimum = timeBin0; for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { Long64_t timeBinOfCount_tmp = it->first; - if (timeBinOfCount_tmp < timeBinMinimum+coincidenceTimeWindowMin_M) continue; // This hit happened too long ago - if (timeBinOfCount_tmp > timeBinMinimum+coincidenceTimeWindowMax_M) break; // This hit happened too late - if ((timeBinOfCount_tmp == timeBinMinimum)&&( ((it->second)->eventEntry) == kEntry_NN)) continue; - // This is the M-counter hit for which we check the pileup -> ignore this hit in this double-hit check. - + // std::cout<<"timeBin0="< S T O P !!!\n"; exit(1);} for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { - Int_t eventNumber = (it->second)->eventIDnumber; + // Int_t eventNumber = (it->second)->eventIDnumber; Long64_t timeBinOfCount_tmp = it->first; if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) { if (bool_debugingRequired) { @@ -262,7 +268,7 @@ Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_ } // Hit candidate was found. Now check its coincidences and vetos - positronQuality=1; + Bool_t bool_coincidenceConditions = true; for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) { if (bool_debugingRequired) { if (debugEventMap[evtID]>4) { (itCounter->second)->myPrintThisCounter(evtID); } @@ -272,7 +278,9 @@ Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_ if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"< remove the candidate. - goto endOfThisHit; // no coincidence found ==> skip hit + bool_coincidenceConditions = false; + goto CoincidencesChecked; + // goto endOfThisHit; // no coincidence found ==> skip hit } } for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) { @@ -281,10 +289,17 @@ Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_ if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"< remove the candidate. - goto endOfThisHit; // coincidence with veto found ==> skip hit + bool_coincidenceConditions = false; + goto CoincidencesChecked; + // goto endOfThisHit; // coincidence with veto found ==> skip hit } } + CoincidencesChecked: + if (!bool_coincidenceConditions) continue; // This hit does not fulfill coincidence and veto criteria + if (positronQuality==2) return 3; // An electron was already found before, and now again ==> double hit + else positronQuality=2; + kEntry = (it->second)->eventEntry; idet = (it->second)->det_i; idetID = (it->second)->det_id; @@ -294,11 +309,6 @@ Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_ if (bool_debugingRequired) { if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"< (previousClock+musrAnalysis::clock_interval)) { +// previousClock += musrAnalysis::clock_interval; +// // dumpFile<<"-1\t"<=musrAnalysis::rewindTimeBins) tdcBin -= musrAnalysis::rewindTimeBins; + else if (tdcBin<0) tdcBin += musrAnalysis::rewindTimeBins; + // Int_t tdc = (tdcBin +#include #include #include #include @@ -55,17 +57,21 @@ class musrCounter { 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(); - void FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID); + void FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID, Int_t eventNum); void RemoveHitsInCounter(Long64_t timeBinLimit); - void RewindHitsInCounter(Long64_t timeBinsToRewind); - Bool_t IsInCoincidence(Long64_t timeBin, char motherCounter, Bool_t ignoreHitsAtBinZero=false, Long64_t timeBinMinimum=-123456789, Long64_t timeBinMaximum=-123456789); + // void RewindHitsInCounter(Long64_t timeBinsToRewind); + void RewindHitsInCounter(); + Bool_t IsInCoincidence(Long64_t timeBin, char motherCounter); Bool_t GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep); - Bool_t CheckForPileupMuons(Long64_t timeBin0, Int_t kEntry_NN); + Bool_t CheckForPileupMuons(Long64_t timeBin0); 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); void myPrintThisCounter(Int_t evtID, Int_t detail=2); + // void CheckClockInfo(Long64_t timeBin); Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;} Long64_t GetNumberOfMuonCandidatesAfterVK(){return numberOfMuonCandidatesAfterVK;} Long64_t GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval(){return numberOfMuonCandidatesAfterVKandDoubleHitRemoval;} + void DumpInfoToDumpFile(Int_t eventNr, Int_t detID, Long64_t tdcBin); + void WriteRewindIntoDumpFile(); private: // static musrCounter* pointerToAnalysis; @@ -93,10 +99,14 @@ class musrCounter { Long64_t numberOfMuonCandidates; Long64_t numberOfMuonCandidatesAfterVK; Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval; - + public: static Bool_t bool_ignoreUnperfectMuons; static Bool_t bool_ignoreUnperfectPositrons; + static Bool_t bool_WriteDataToDumpFile; + static ofstream dumpFile; + // static Long64_t previousClock; + // static Long64_t CLOCK_INTERVAL; }; #endif diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 046b8e4..3be7592 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -960,6 +960,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() { sscanf(&line[0],"%*s %*s %*s %lf",&sigma); myMusrScintSD ->SetAPDcellsTimeVariationSigma(sigma*nanosecond); } + else if (strcmp(varName,"SetAPDcrossTalk")==0) { + double crossTalkProb; + sscanf(&line[0],"%*s %*s %*s %lf",&crossTalkProb); + myMusrScintSD ->SetAPDcrossTalk(crossTalkProb); + } else { G4cout<<"musrDetectorConstruction.cc: ERROR: Unknown parameterName \"" <GetDefinition()->GetParticleName(); // if (particleName=="opticalphoton") {G4cout<<"UFON JE TU: edep="<GetPostStepPoint()->GetProcessDefinedStep(); // G4String processName = (process) ? process->GetProcessName() : "Unknown"; @@ -269,29 +272,28 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep, G4double APDcellsTimeVaria // 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 (APDcellsTimeVariationRequested) tmpTime += APDcellsTimeVariation; - // if (APDcellsTimeVariationRequested) tmpTime += G4RandGauss::shoot(0,APDcellsTimeVariationSigma); - - if (APDcellsEffectRequested) { - APDcellID = FindAPDcellID(aStep); - // G4cout<<"Cell ID="<begin(); it2 != optHitDetectorMap->end(); it2++ ) { - if ((it2->second)==APDcellID) { - // G4cout<<"Already fired cell ="<first)) { - optHitDetectorMap->erase(it2); - optHitDetectorMap->insert(std::pair(tmpTime,APDcellID)); - } - break; - } - } - } - if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair(tmpTime,APDcellID)); + // G4bool APDcellAlreadyFired = false; + G4int APDcellID = (APDcellsEffectRequested) ? FindAPDcellID(aStep) : 0; + FireAPDcell(optHitDetectorMap,APDcellID,tmpTime); + // if (APDcellsEffectRequested) { + // APDcellID = FindAPDcellID(aStep); + // // G4cout<<"Cell ID="<begin(); it2 != optHitDetectorMap->end(); it2++ ) { + // if ((it2->second)==APDcellID) { + // // G4cout<<"Already fired cell ="<first)) { + // optHitDetectorMap->erase(it2); + // optHitDetectorMap->insert(std::pair(tmpTime,APDcellID)); + // } + // break; + // } + // } + // } + // if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair(tmpTime,APDcellID)); + // // G4cout<<" tmpTime="<empty()) continue; + + // Do APD cell variation time if requested. Even if not requested, generate the random numbers in order to + // have reproducible simulation. + for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) { + G4double APDcellsTimeVariation = G4RandGauss::shoot(0,APDcellsTimeVariationSigma); + if (APDcellsTimeVariationRequested) { + G4int APDcellID = it2->second; + G4double tmpTime = it2->first + APDcellsTimeVariation; + optHitDetectorMap->erase(it2); + optHitDetectorMap->insert(std::pair(tmpTime,APDcellID)); + } + } + + // Simulate cross talk if requested. Even if not requested, generate the random numbers in order to + // have reproducible simulation. + for (optHitDetectorMapType::const_iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) { + G4double rand = G4UniformRand(); + if (APDcrossTalkRequested) { + if (randsecond; + G4double tmpTime = it2->first; + if (!APDcellsEffectRequested) FireAPDcell(optHitDetectorMap,APDcellID,tmpTime); + else { + if (rand<(0.5*APDcrossTalkProb)) APDcellID--; + else APDcellID++; + if ((APDcellID<0)||(APDcellID>=APDcell_nx*APDcell_ny*APDcell_nz)) continue; // the cross-talk cell is outside range + FireAPDcell(optHitDetectorMap,APDcellID,tmpTime); + } + } + } + } + optHitDetectorMapType::const_iterator it2_START = optHitDetectorMap->begin(); optHitDetectorMapType::const_iterator it2_STOP = it2_START; optHitDetectorMapType::const_iterator it2_LAST = optHitDetectorMap->end(); it2_LAST--; @@ -908,6 +942,28 @@ G4int musrScintSD::FindAPDcellID(G4Step* aStep) { //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time) { + if (!APDcellsEffectRequested) { + optHitDetectorMap->insert(std::pair(time,0)); + } + else { + G4bool APDcellAlreadyFired = false; + for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) { + if ((it2->second)==APDcellID) { // this cell already fired before ==> check times + APDcellAlreadyFired = true; + if (time<(it2->first)) { // the new cell fired before the already saved cell ==> replace it + optHitDetectorMap->erase(it2); + optHitDetectorMap->insert(std::pair(time,APDcellID)); + } + break; + } + } + if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair(time,APDcellID)); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + void musrScintSD::FindCFDtime(G4double& OPSA_CFD_time, G4double& OPSA_CFD_ampl, G4double timeOfFirstPhoton) { OPSA_CFD->Reset("M"); for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of electronic signal histogram