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

28.10.2011 Kamil Sedlak

Browse Source 
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).
This commit is contained in:
sedlak
2011-10-28 14:04:11 +00:00
parent 89c6f27ae1
commit f93a3c070b
7 changed files with 306 additions and 165 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
include/musrScintSD.hh
View File

@ -105,7 +105,7 @@ class musrScintSD : public G4VSensitiveDetector
OPSAhistoNbin=nBins; OPSAhistoMin=min; OPSAhistoMax=max; OPSAhistoNbin=nBins; OPSAhistoMin=min; OPSAhistoMax=max;
OPSAhistoBinWidth=(max-min)/nBins; OPSAhistoBinWidth1000=OPSAhistoBinWidth*1000; OPSAhistoBinWidth=(max-min)/nBins; OPSAhistoBinWidth1000=OPSAhistoBinWidth*1000;
} }
void ProcessOpticalPhoton(G4Step* aStep, G4double APDcellsTimeVariation); void ProcessOpticalPhoton(G4Step* aStep);
void EndOfEvent_OptiacalPhotons(); void EndOfEvent_OptiacalPhotons();
void ReadInPulseShapeArray(const char* filename); void ReadInPulseShapeArray(const char* filename);
void FindCFDtime(G4double& OPSA_CFD_time, G4double& OPSA_CFD_ampl, G4double timeOfFirstPhoton); 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; if (sigma!=0) APDcellsTimeVariationRequested = true;
APDcellsTimeVariationSigma=sigma; APDcellsTimeVariationSigma=sigma;
} }
void SetAPDcrossTalk(G4double crosstalkProb) {
if (crosstalkProb>0) APDcrossTalkRequested = true;
APDcrossTalkProb = crosstalkProb;
}
private: private:
static musrScintSD* pointer; static musrScintSD* pointer;
@ -168,6 +172,12 @@ class musrScintSD : public G4VSensitiveDetector
G4double APDcell_ax, APDcell_ay, APDcell_az; // dimensions of APD cells G4double APDcell_ax, APDcell_ay, APDcell_az; // dimensions of APD cells
G4bool APDcellsTimeVariationRequested; // simmulate effect of detection time variations between differenct cells G4bool APDcellsTimeVariationRequested; // simmulate effect of detection time variations between differenct cells
G4double APDcellsTimeVariationSigma; // sigma of the 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...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

165
musrSimAna/musrAnalysis.cxx
View File

@ -100,7 +100,7 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
pcoincwin = 40; pcoincwin = 40;
vcoincwin = 80; vcoincwin = 80;
muonRateFactor = 1.; muonRateFactor = 1.;
rewindTimeBins = 1000000000; // rewindTimeBins = 1000000000;
dataWindowMin = -0.2; dataWindowMin = -0.2;
dataWindowMax = 10.0; dataWindowMax = 10.0;
pileupWindowMin = -10.5; pileupWindowMin = -10.5;
@ -176,7 +176,9 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
} }
else if (strncmp(tmpString0,"REWINDTIMEBINS=",strlen("REWINDTIMEBINS="))==0) { else if (strncmp(tmpString0,"REWINDTIMEBINS=",strlen("REWINDTIMEBINS="))==0) {
sscanf(&line[strlen("REWINDTIMEBINS=")],"%d",&tmpivar); sscanf(&line[strlen("REWINDTIMEBINS=")],"%d",&tmpivar);
rewindTimeBins = tmpivar; std::cout<<" JUHELAK: rewindTimeBins="<<rewindTimeBins<<std::endl; if (tmpivar<0) {rewindTimeBins = Long64_t(tmpivar)* Long64_t(tmpivar);}
else {rewindTimeBins = tmpivar;}
std::cout<<" JUHELAK: rewindTimeBins="<<rewindTimeBins<<std::endl;
} }
else if (strncmp(tmpString0,"DATAWINDOWMIN=",strlen("DATAWINDOWMIN="))==0) { else if (strncmp(tmpString0,"DATAWINDOWMIN=",strlen("DATAWINDOWMIN="))==0) {
sscanf(&line[strlen("DATAWINDOWMIN=")],"%g",&tmpfvar); sscanf(&line[strlen("DATAWINDOWMIN=")],"%g",&tmpfvar);
@ -218,6 +220,23 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
bool_debugingRequired=true; bool_debugingRequired=true;
debugEventMap.insert(std::pair<int,int>(ieventToDebug_tmp,iLevelToDebug_tmp)); debugEventMap.insert(std::pair<int,int>(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="<<tmpString1<<std::endl;
if ((strcmp(tmpString1,"")!=0)&&(strcmp(tmpString1,"Unset")!=0)) musrCounter::dumpFile.open(tmpString1);
// musrCounter::dumpFile.open("dumpFile.txt");
if (! (musrCounter::dumpFile.is_open()) ) {
std::cout<<"Writing data into a DumpFile was requested, but the dump file can not be opened!"<<std::cout;
std::cout<<" ===> STOP !!!"<<std::endl;
exit(1);
}
}
else if (strcmp(tmpString0,"USE_UNPERFECT_POSITRONS_IN_DOUBLE_HITS")==0) { else if (strcmp(tmpString0,"USE_UNPERFECT_POSITRONS_IN_DOUBLE_HITS")==0) {
musrCounter::bool_ignoreUnperfectPositrons = false; musrCounter::bool_ignoreUnperfectPositrons = false;
} }
@ -697,11 +716,6 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
jj+=strlen(tmpCoinc); jj+=strlen(tmpCoinc);
int iCoinc=atoi(tmpCoinc); int iCoinc=atoi(tmpCoinc);
if (iCoinc==987654321) continue; if (iCoinc==987654321) continue;
//ckdel 20.9.2010 if ((maxChannels<=tmpChannel)||(maxChannels<=iCoinc)) {
//ckdel 20.9.2010 std::cout<<" Channel number ("<<tmpChannel<<") or channel coincidence required ("<<iCoinc
//ckdel 20.9.2010 <<") larger than maxChannels ("<<maxChannels<<") ==> S T O P"<<std::endl;
//ckdel 20.9.2010 exit(1);
//ckdel 20.9.2010 }
tmpCoincidenceMultimap.insert(std::pair<int, int>(tmpChannel,iCoinc)); tmpCoincidenceMultimap.insert(std::pair<int, int>(tmpChannel,iCoinc));
// tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc; // tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc;
} while (jj<endOfCoincidenceArray); } while (jj<endOfCoincidenceArray);
@ -766,13 +780,13 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
char DetectorType = (it->second)->GetCounterType(); char DetectorType = (it->second)->GetCounterType();
if (DetectorType=='M') { if (DetectorType=='M') {
(it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin); (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)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('M',-mcoincwin+pileupWindowBinMin,-mcoincwin,mcoincwin);
(it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin); (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin);
} }
else if (DetectorType=='P') { else if (DetectorType=='P') {
(it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin); (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)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('P',-pcoincwin+dataWindowBinMin,-pcoincwin,pcoincwin);
(it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin); (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin);
} }
@ -835,7 +849,7 @@ void musrAnalysis::CreateHistograms() {
hInfo->Fill(1, fieldValue); hInfo->Fill(1, fieldValue);
hInfo->Fill(6, runID); hInfo->Fill(6, runID);
hInfo->Fill(7, hGeantParameters->GetBinContent(7)); hInfo->Fill(7, hGeantParameters->GetBinContent(7));
std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<<fieldValue<<"T, omega="<<omega<<std::endl; std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<<fieldValue<<"T, omega="<<omega<<std::endl;
} }
//================================================================ //================================================================
@ -884,6 +898,9 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
fHistograms->Close(); fHistograms->Close();
delete fHistograms; delete fHistograms;
if (musrCounter::bool_WriteDataToDumpFile==true) {
musrCounter::dumpFile.close();
}
//============================== //==============================
@ -956,7 +973,8 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
if (!boolInfinitelyLowMuonRate) RemoveOldHitsFromCounters(currentTimeBin-1); // Remove all hits that can not influence the next event if (!boolInfinitelyLowMuonRate) RemoveOldHitsFromCounters(currentTimeBin-1); // Remove all hits that can not influence the next event
else RemoveOldHitsFromCounters(std::numeric_limits<Long64_t>::max()/2); else RemoveOldHitsFromCounters(std::numeric_limits<Long64_t>::max()/2);
if (currentTimeBin>rewindTimeBins) { // When the time variables start to be too large, rewind the time info everywhere; if (currentTimeBin>rewindTimeBins) { // When the time variables start to be too large, rewind the time info everywhere;
RewindAllTimeInfo(rewindTimeBins); // RewindAllTimeInfo(rewindTimeBins);
RewindAllTimeInfo();
} }
} }
@ -974,14 +992,23 @@ void musrAnalysis::RemoveOldHitsFromCounters(Long64_t timeBinLimit) {
//================================================================ //================================================================
//void musrAnalysis::RewindAllTimeInfo(Double_t timeToRewind) { //void musrAnalysis::RewindAllTimeInfo(Double_t timeToRewind) {
void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) { //void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) {
void musrAnalysis::RewindAllTimeInfo() {
Long64_t timeBinsToRewind = rewindTimeBins;
if (musrCounter::bool_WriteDataToDumpFile) {
// Long64_t currentTimeBin = Long64_t(currentTime/tdcresolution);
// mCounter->CheckClockInfo(currentTimeBin);
mCounter->WriteRewindIntoDumpFile();
// musrCounter::previousClock -= timeBinsToRewind;
}
currentTime -= timeBinsToRewind*tdcresolution; currentTime -= timeBinsToRewind*tdcresolution;
nextUnfilledEventTime -= timeBinsToRewind*tdcresolution; nextUnfilledEventTime -= timeBinsToRewind*tdcresolution;
numberOfRewinds++; numberOfRewinds++;
// Loop over all timing information and rewind it by "timeToRewind"; // Loop over all timing information and rewind it by "timeToRewind";
// std::cout<<"musrAnalysis::RewindAllTimeInfo()"<<std::endl; // std::cout<<"musrAnalysis::RewindAllTimeInfo()"<<std::endl;
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) { for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
(*it).second->RewindHitsInCounter(timeBinsToRewind); // (*it).second->RewindHitsInCounter(timeBinsToRewind);
(*it).second->RewindHitsInCounter();
} }
} }
@ -1019,7 +1046,8 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl; // std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM); // mCounterHitExistsForThisEventID = mCounter->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 = 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; timeBinOfThePreviouslyProcessedHit = timeBin0;
//___________________________________________________________ //___________________________________________________________
@ -1225,7 +1253,8 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl; // std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM); // mCounterHitExistsForThisEventID = mCounter->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 = 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; timeBinOfThePreviouslyProcessedHit = timeBin0;
// if (mCounterHitExistsForThisEventID) std::cout<<" YYYYYYYYYYYYYYYYYYY check this : LOOOPING AGAIN"<<std::endl; // if (mCounterHitExistsForThisEventID) std::cout<<" YYYYYYYYYYYYYYYYYYY check this : LOOOPING AGAIN"<<std::endl;
oncePerEvent=false; oncePerEvent=false;
@ -1322,37 +1351,14 @@ Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
else { else {
// Double_t omega = 851.372*fieldNomVal[0]; // Double_t omega = 851.372*fieldNomVal[0];
Double_t dPhaseShift = (omega==0) ? 0:phaseShiftMap[det_ID[i]]/omega; Double_t dPhaseShift = (omega==0) ? 0:phaseShiftMap[det_ID[i]]/omega;
//cDEL if (det_ID[i]<20) std::cout<<"phaseShift = "<<phaseShiftMap[det_ID[i]]<<" dPhaseShift="<<dPhaseShift<<" tdcresolution="<<tdcresolution<<std::endl;
Double_t t1,t2; Double_t t1 = (globTime+det_time_start[i] )/tdcresolution;
if (bool_muDecayTimeTransformation) { Double_t t2 = (globTime+det_time_start[i]+dPhaseShift)/tdcresolution;
// THIS OPTION WAS SUPOSED TO ALLOW THE USER TO SOMEHOW (IN A TRICKY AND NOT 100% RELIABLE WAY) // }
// SHIFT THE POSITRON COUNTS TO A RESTRICTED TIME WINDOW, AS IF THE MUON DECAYED IN SOME RESTRICTED
// TIME INTERVAL. THIS, HOWEVER, DOES NOT WORK, IN THE WAY IT IS IMPLEMENTED HERE, BECAUSE THE
// MUON SPIN ROTATION AT REST IS IGNORED HERE - IT WOULD NEED SOME FURTHER WORK TO DO, AND
// BECOMES DIFFICULT TO INTERPRET LATER ==> WORK ON THIS WAS STOPPED.
Double_t diff = muDecayTime_t_max - muDecayTime_t_min;
Double_t ttt1 = det_time_start[i];
if ((ttt1 > muDecayTime_Transformation_min) && (ttt1 < muDecayTime_Transformation_max)) {
if (muDecayTime < muDecayTime_t_min) {
ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)+1) * diff + det_time_start[i];
}
else if (muDecayTime > muDecayTime_t_max) {
ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)) * diff + det_time_start[i];
}
}
t1 = (globTime+ttt1 )/tdcresolution;
t2 = (globTime+ttt1+dPhaseShift)/tdcresolution;
// std::cout<<"DEBUG: det_time_start[i]="<<det_time_start[i]<<" ttt1="<<ttt1<<" t1="<<t1<<" t2="<<t2<<std::endl;
}
else {
t1 = (globTime+det_time_start[i] )/tdcresolution;
t2 = (globTime+det_time_start[i]+dPhaseShift)/tdcresolution;
}
Long64_t timeBin = Long64_t(t1); Long64_t timeBin = Long64_t(t1);
Long64_t timeBin2 = Long64_t(t2); Long64_t timeBin2 = Long64_t(t2);
(*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i]); (*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i],eventID);
} }
} }
@ -1363,14 +1369,14 @@ Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
} }
//================================================================ //================================================================
Bool_t musrAnalysis::MuonCounterHit(Int_t evID, Long64_t timeBinMin, Long64_t& timeBin0, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep) { //Bool_t musrAnalysis::MuonCounterHit(Int_t evID, Long64_t timeBinMin, Long64_t& timeBin0, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep) {
Bool_t mCounterHitCanditateExists = mCounter->GetNextGoodMuon(evID,timeBinMin,timeBin0,kEntry,idet,idetID,idetEdep); // Bool_t mCounterHitCanditateExists = mCounter->GetNextGoodMuon(evID,timeBinMin,timeBin0,kEntry,idet,idetID,idetEdep);
if (!mCounterHitCanditateExists) return false; // if (!mCounterHitCanditateExists) return false;
// Check for other muons within the pileup window: // // Check for other muons within the pileup window:
if ( mCounter->CheckForPileupMuons(timeBin0,kEntry) ) return false; // This muon candidate is killed due to a double hit rejection. // if ( mCounter->CheckForPileupMuons(timeBin0,kEntry) ) return false; // This muon candidate is killed due to a double hit rejection.
return true; // return true;
} //}
//
//================================================================ //================================================================
Bool_t musrAnalysis::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) { Bool_t musrAnalysis::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) {
@ -1379,39 +1385,48 @@ Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_
} }
if (pCounterMap.empty()) return false; if (pCounterMap.empty()) return false;
Bool_t positronHitFound = false; // Bool_t positronHitFound = false;
Bool_t goodPositronFound = false; Bool_t goodPositronFound = false;
Int_t positronQuality = 0; Int_t positronQuality = 0;
// std::cout<<"Debug 10------------------------------"<<std::endl; // std::cout<<"Debug 10------------------------------"<<std::endl;
if (musrMode=='D') { if (musrMode=='D') {
// Loop over all positron counters // Loop over all positron counters
for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) { for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) {
// std::cout<<"Debug 20"<<std::endl; positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep);
Long64_t dataBinMinimum = dataBinMin; if (positronQuality==3) return false; // double hit was found in the same counter
do { if (positronQuality==2) {
// std::cout<<"Debug 30"<<std::endl; if (goodPositronFound) return false; // double hit was found in a different counter
positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMinimum,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); goodPositronFound = true;
// std::cout<<"Debug 40 positronQuality="<<positronQuality<<std::endl; }
dataBinMinimum=tBin1+1; // If a positron is found, one needs to search again through the remaining hits in this counter.
if (positronQuality>0) {
// std::cout<<"Debug 50"<<std::endl;
if (positronHitFound) { // Positron found now but also previously ==> double hit ==> through away this event.
// std::cout<<"Debug 60"<<std::endl;
if (bool_debugingRequired) {
if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID
<<"\"PositronCounterHit\": Coincidence with other positron candidate in other counter."<<std::endl;}
}
return false;
} // end of "positronHitFound"
// std::cout<<"Debug 70"<<std::endl;
positronHitFound = true;
if (positronQuality>1) goodPositronFound = true;
}
// std::cout<<"Debug 80"<<std::endl;
} while (positronQuality>0);
// std::cout<<"Debug 90"<<std::endl;
} }
// Long64_t dataBinMinimum = dataBinMin;
// do {
// // std::cout<<"Debug 30"<<std::endl;
// positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMinimum,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep);
// // std::cout<<"Debug 40 positronQuality="<<positronQuality<<std::endl;
// dataBinMinimum=tBin1+1; // If a positron is found, one needs to search again through the remaining hits in this counter.
// if (positronQuality>0) {
// // std::cout<<"Debug 50"<<std::endl;
// if (positronHitFound) { // Positron found now but also previously ==> double hit ==> through away this event.
// // std::cout<<"Debug 60"<<std::endl;
// if (bool_debugingRequired) {
// if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID
// <<"\"PositronCounterHit\": Coincidence with other positron candidate in other counter."<<std::endl;}
// }
// return false;
// } // end of "positronHitFound"
// // std::cout<<"Debug 70"<<std::endl;
// positronHitFound = true;
// if (positronQuality>1) goodPositronFound = true;
// }
// // std::cout<<"Debug 80"<<std::endl;
// } while (positronQuality>0);
// // std::cout<<"Debug 90"<<std::endl;
// }
// std::cout<<"Debug 100 positronQuality="<<positronQuality<<std::endl; // std::cout<<"Debug 100 positronQuality="<<positronQuality<<std::endl;
if (goodPositronFound) return true; if (goodPositronFound) return true;
} }
// std::cout<<"Debug 110"<<std::endl; // std::cout<<"Debug 110"<<std::endl;

29
musrSimAna/musrAnalysis.hh
View File

@ -8,6 +8,8 @@
#ifndef musrAnalysis_h #ifndef musrAnalysis_h
#define musrAnalysis_h #define musrAnalysis_h
#include <iostream>
#include <fstream>
#include <list> #include <list>
#include <map> #include <map>
#include <TROOT.h> #include <TROOT.h>
@ -293,12 +295,13 @@ public :
virtual void SaveHistograms(char* runChar,char* v1190FileName); virtual void SaveHistograms(char* runChar,char* v1190FileName);
virtual void RemoveOldHitsFromCounters(Long64_t timeBinLimit); virtual void RemoveOldHitsFromCounters(Long64_t timeBinLimit);
// virtual void RewindAllTimeInfo(Double_t timeToRewind); // 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 PrintHitsInAllCounters();
virtual void InitialiseEvent(); virtual void InitialiseEvent();
virtual Double_t PreprocessEvent(Long64_t iEn); 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 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 CopySubstring(char* inputChar,int iStart,int iEnd,char* outputChar);
void MyPrintTree(); void MyPrintTree();
void MyPrintConditions(); void MyPrintConditions();
@ -328,17 +331,12 @@ public :
Int_t pcoincwin; Int_t pcoincwin;
Int_t vcoincwin; Int_t vcoincwin;
Double_t muonRateFactor; Double_t muonRateFactor;
Long64_t rewindTimeBins;
Double_t dataWindowMin; Double_t dataWindowMin;
Double_t dataWindowMax; Double_t dataWindowMax;
Double_t pileupWindowMin; Double_t pileupWindowMin;
Double_t pileupWindowMax; Double_t pileupWindowMax;
Double_t promptPeakWindowMin; Double_t promptPeakWindowMin;
Double_t promptPeakWindowMax; Double_t promptPeakWindowMax;
Long64_t pileupWindowBinMin;
Long64_t pileupWindowBinMax;
Long64_t dataWindowBinMin;
Long64_t dataWindowBinMax;
Int_t overallBinDelay; Int_t overallBinDelay;
Bool_t boolInfinitelyLowMuonRate; Bool_t boolInfinitelyLowMuonRate;
@ -378,6 +376,13 @@ public :
public: public:
static const Int_t nrConditions = 31; static const Int_t nrConditions = 31;
Bool_t condition[nrConditions]; 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: private:
// HISTOGRAMS: // HISTOGRAMS:
@ -417,6 +422,16 @@ private:
#endif #endif
#ifdef musrAnalysis_cxx #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) musrAnalysis::musrAnalysis(TTree *tree)
{ {
variableMap["muDecayPosX"]=&muDecayPosX; variableMap["muDecayPosX"]=&muDecayPosX;

130
musrSimAna/musrCounter.cxx
View File

@ -1,6 +1,7 @@
//#include <iostream> //#include <iostream>
#include "musrCounter.hh" #include "musrCounter.hh"
#include "TCanvas.h" #include "TCanvas.h"
#include "musrAnalysis.hh"
typedef std::map<int,int> debugEventMapType; typedef std::map<int,int> debugEventMapType;
debugEventMapType debugEventMap; debugEventMapType debugEventMap;
@ -8,6 +9,10 @@ Bool_t bool_debugingRequired;
Bool_t musrCounter::bool_ignoreUnperfectMuons = true; Bool_t musrCounter::bool_ignoreUnperfectMuons = true;
Bool_t musrCounter::bool_ignoreUnperfectPositrons = 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...... //....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="<<timeBin<<" timeBin2="<<timeBin2<<" counterTimeShift="<< counterTimeShift<<std::endl; //cDEL std::cout<<"FillHitInCounter: timeBin="<<timeBin<<" timeBin2="<<timeBin2<<" counterTimeShift="<< counterTimeShift<<std::endl;
//cDEL std::cout<<" FillHitInCounter I: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl; //cDEL std::cout<<" FillHitInCounter I: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl;
// std::cout<<"musrCounter::FillHitInCounter:"<<counterNr<<std::endl; // std::cout<<"musrCounter::FillHitInCounter:"<<counterNr<<std::endl;
@ -77,6 +82,11 @@ void musrCounter::FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t tim
hitInfo* hInfo = new hitInfo(kEntry,eveID,iDet,detectorID,edep,timeBin2-counterTimeShift); hitInfo* hInfo = new hitInfo(kEntry,eveID,iDet,detectorID,edep,timeBin2-counterTimeShift);
//cDEL std::cout<<detectorID<<" FillHitInCounter II: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl; //cDEL std::cout<<detectorID<<" FillHitInCounter II: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl;
hitMap.insert( std::pair<Long64_t,hitInfo*>(timeBin-counterTimeShift,hInfo) ); hitMap.insert( std::pair<Long64_t,hitInfo*>(timeBin-counterTimeShift,hInfo) );
if (bool_WriteDataToDumpFile) { // write data into an output file
// CheckClockInfo(timeBin);
// dumpFile<<eventNum<<"\t"<<detectorID<<"\t"<<timeBin<<"\n";
DumpInfoToDumpFile(eventNum,detectorID,timeBin);
}
} }
} }
@ -99,55 +109,46 @@ void musrCounter::RemoveHitsInCounter(Long64_t timeBinLimit) {
} }
//================================================================ //================================================================
void musrCounter::RewindHitsInCounter(Long64_t timeBinsToRewind) { //void musrCounter::RewindHitsInCounter(Long64_t timeBinsToRewind) {
void musrCounter::RewindHitsInCounter() {
// Reset time in hits from the Counter class // Reset time in hits from the Counter class
if (hitMap.empty()) return; if (hitMap.empty()) return;
// Long64_t timeBinsToRewind = musrAnalysis::rewindTimeBins;
hitMap_TYPE hitMap_TMP; hitMap_TYPE hitMap_TMP;
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
// std::cout<<"musrCounter::RewindHitsInCounter: "<<std::endl;
Long64_t tempBinT = it->first; Long64_t tempBinT = it->first;
// int tempEvnr= it->second;
// hitMap_TMP.insert( std::pair<Long64_t,int>(tempBinT-timeBinsToRewind,tempEvnr) );
hitInfo* tempEvnr= it->second; hitInfo* tempEvnr= it->second;
tempEvnr->RewindTimeBin2(timeBinsToRewind); tempEvnr->RewindTimeBin2(musrAnalysis::rewindTimeBins);
hitMap_TMP.insert( std::pair<Long64_t,hitInfo*>(tempBinT-timeBinsToRewind,tempEvnr) ); hitMap_TMP.insert( std::pair<Long64_t,hitInfo*>(tempBinT-musrAnalysis::rewindTimeBins,tempEvnr) );
} }
hitMap.swap(hitMap_TMP); 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 // 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; if (hitMap.empty()) return false;
// Long64_t timeBinToTest = timeBin;
Long64_t timeBinMin; Long64_t timeBinMin;
Long64_t timeBinMax; Long64_t timeBinMax;
// If timeBinMinimum and timeBinMaximum are not specified, use internal time window of the detector (koincidence or veto detectors). // 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). // 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" if (counterType == 'V') timeBinMin = timeBin + antiCoincidenceTimeWindowMin; // this is veto detector
else 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=='M') timeBinMin = timeBin + coincidenceTimeWindowMin_M; // this is coinc. detector (or M counter) connected to M
else if (motherCounter=='P') timeBinMin = timeBin + coincidenceTimeWindowMin_P; // this is coinc. detector connected to P 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" if (counterType == 'V') timeBinMax = timeBin + antiCoincidenceTimeWindowMax; // this is veto detector
else 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=='M') timeBinMax = timeBin + coincidenceTimeWindowMax_M; // this is coinc. detector (or M counter) connected to M
else if (motherCounter=='P') timeBinMax = timeBin + coincidenceTimeWindowMax_P; // this is coinc. detector connected to P 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) { for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
Long64_t timeBinOfCount_tmp = it->first; Long64_t timeBinOfCount_tmp = it->first;
if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) { if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) {
if ((timeBin!=timeBinOfCount_tmp)||(!ignoreHitsAtBinZero)) { // if ((timeBin!=timeBinOfCount_tmp)||(!ignoreHitsAtBinZero)) {
return true; return true;
} // }
} }
else if (timeBinOfCount_tmp > timeBinMax) return false; 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; Int_t eventNumber = (it->second)->eventIDnumber;
if (eventNumber!=evtID) continue; // This trigger hit does not correspond to the currently processed event 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++; numberOfMuonCandidates++;
// std::cout<<"*** "<<evtID<<" eventNumber="<<eventNumber<<" canditas="<<numberOfMuonCandidates<<std::endl;
// Hit candidate was found. Now check its coincidences and vetos // Hit candidate was found. Now check its coincidences and vetos
Bool_t bool_coincidenceConditions = true;
for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) { for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) {
if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') )) { // no coincidence found ==> skip hit if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') )) { // no coincidence found ==> skip hit
if (bool_ignoreUnperfectMuons) hitMap.erase(it); if (bool_ignoreUnperfectMuons) hitMap.erase(it);
goto endOfThisHit; bool_coincidenceConditions = false;
goto MuonCoincidencesChecked;
} }
} }
for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) { for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) {
if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') ) { // coincidence with veto found ==> skip hit if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') ) { // coincidence with veto found ==> skip hit
if (bool_ignoreUnperfectMuons) hitMap.erase(it); if (bool_ignoreUnperfectMuons) hitMap.erase(it);
goto endOfThisHit; bool_coincidenceConditions = false;
goto MuonCoincidencesChecked;
} }
} }
MuonCoincidencesChecked:
if (!bool_coincidenceConditions) continue; // This hit does not fulfill coincidence and veto criteria
numberOfMuonCandidatesAfterVK++; numberOfMuonCandidatesAfterVK++;
if ( CheckForPileupMuons(timeBinOfNextHit) ) {
// std::cout<<"CheckForPileupMuons=true"<<std::endl;
continue; // This muon candidate is killed due to a double hit rejection.
}
kEntry = (it->second)->eventEntry; kEntry = (it->second)->eventEntry;
idet = (it->second)->det_i; idet = (it->second)->det_i;
idetID = (it->second)->det_id; idetID = (it->second)->det_id;
idetEdep = (it->second)->det_edep; idetEdep = (it->second)->det_edep;
// timeBinOfNextHit = timeBinOfCount_tmp;
return true; return true;
endOfThisHit:
continue;
} }
return false; 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. // Check for pileup muons. If double hit in M-counter is found, return true.
Long64_t timeBinMinimum = timeBin0; Long64_t timeBinMinimum = timeBin0;
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
Long64_t timeBinOfCount_tmp = it->first; Long64_t timeBinOfCount_tmp = it->first;
if (timeBinOfCount_tmp < timeBinMinimum+coincidenceTimeWindowMin_M) continue; // This hit happened too long ago // std::cout<<"timeBin0="<<timeBin0<<" timeBinOfCount_tmp="<<timeBinOfCount_tmp<<std::endl;
if (timeBinOfCount_tmp > timeBinMinimum+coincidenceTimeWindowMax_M) break; // This hit happened too late if (timeBinOfCount_tmp < timeBinMinimum+musrAnalysis::pileupWindowBinMin) continue; // This hit happened too long ago
if ((timeBinOfCount_tmp == timeBinMinimum)&&( ((it->second)->eventEntry) == kEntry_NN)) continue; if (timeBinOfCount_tmp > timeBinMinimum+musrAnalysis::pileupWindowBinMax) break; // This hit happened too late
// This is the M-counter hit for which we check the pileup -> ignore this hit in this double-hit check. // if ((timeBinOfCount_tmp == timeBinMinimum)&&( ((it->second)->eventEntry) == kEntry_NN)) continue;
if (timeBinOfCount_tmp == timeBinMinimum) continue;
// This is the M-counter hit for which we check the pileup -> ignore this hit in this double-hit check.
// We have found the hit, which could be the double hit. If it is required, we now have // We have found the hit, which could be the double hit. If it is required, we now have
// to check the coincidences and anticoincidences of this double-hit muon: // to check the coincidences and anticoincidences of this double-hit muon:
if (!bool_ignoreUnperfectMuons) { if (!bool_ignoreUnperfectMuons) {
@ -244,15 +250,15 @@ Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_
// INPUT PARAMETERS: evtID, timeBinMin // INPUT PARAMETERS: evtID, timeBinMin
// OUTPUT PARAMETERS: timeBinOfNextGoodHit // OUTPUT PARAMETERS: timeBinOfNextGoodHit
// positronQuality = 0 ... no positron candidate found // positronQuality = 0 ... no positron candidate found
// = 1 ... positron candidate found, but vetoed or not in a required coincidence with other detector
// = 2 ... good positron found // = 2 ... good positron found
// = 3 ... double hit
// Loop over the hits in the counter // Loop over the hits in the counter
Int_t positronQuality=0; Int_t positronQuality=0;
if (bool_debugingRequired) {if (debugEventMap[evtID]>4) myPrintThisCounter(evtID);} if (bool_debugingRequired) {if (debugEventMap[evtID]>4) myPrintThisCounter(evtID);}
if (hitMap.empty()) return 0; if (hitMap.empty()) return 0;
if (counterType!='P') {std::cout<<"\n!!! FATAL ERROR !!! musrCounter::GetNextGoodPositron: not the positron counter! ==> S T O P !!!\n"; exit(1);} if (counterType!='P') {std::cout<<"\n!!! FATAL ERROR !!! musrCounter::GetNextGoodPositron: not the positron counter! ==> S T O P !!!\n"; exit(1);}
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { 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; Long64_t timeBinOfCount_tmp = it->first;
if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) { if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) {
if (bool_debugingRequired) { 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 // 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) { for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) {
if (bool_debugingRequired) { if (bool_debugingRequired) {
if (debugEventMap[evtID]>4) { (itCounter->second)->myPrintThisCounter(evtID); } 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:"<<evtID<<"\"GetNextGoodPositron\": Coincidence required but not found"<<std::endl;} if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Coincidence required but not found"<<std::endl;}
} }
if (bool_ignoreUnperfectPositrons) hitMap.erase(it); // no coincidence found ==> remove the candidate. if (bool_ignoreUnperfectPositrons) hitMap.erase(it); // no coincidence found ==> 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) { 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:"<<evtID<<"\"GetNextGoodPositron\": Coincidence vith veto detector found"<<std::endl;} if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Coincidence vith veto detector found"<<std::endl;}
} }
if (bool_ignoreUnperfectPositrons) hitMap.erase(it); // coincidence with veto found ==> remove the candidate. if (bool_ignoreUnperfectPositrons) hitMap.erase(it); // coincidence with veto found ==> 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; kEntry = (it->second)->eventEntry;
idet = (it->second)->det_i; idet = (it->second)->det_i;
idetID = (it->second)->det_id; idetID = (it->second)->det_id;
@ -294,11 +309,6 @@ Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_
if (bool_debugingRequired) { if (bool_debugingRequired) {
if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Good positron candidate found in this counter."<<std::endl;} if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Good positron candidate found in this counter."<<std::endl;}
} }
return 2;
endOfThisHit:
;
// continue;
} }
return positronQuality; return positronQuality;
} }
@ -348,3 +358,23 @@ void musrCounter::myPrintThisCounter(Int_t evtID, Int_t detail) {
} }
//================================================================ //================================================================
//void musrCounter::CheckClockInfo(Long64_t timeBin) {
// Int_t clock_detectorID=musrAnalysis::clock_channelID;
// if (timeBin > (previousClock+musrAnalysis::clock_interval)) {
// previousClock += musrAnalysis::clock_interval;
// // dumpFile<<"-1\t"<<clock_detectorID<<"\t"<<previousClock<<"\n";
// DumpInfoToDumpFile(-1,clock_detectorID,previousClock);
// }
//}
//================================================================
void musrCounter::DumpInfoToDumpFile(Int_t eventNr, Int_t detID, Long64_t tdcBin) {
if (tdcBin>=musrAnalysis::rewindTimeBins) tdcBin -= musrAnalysis::rewindTimeBins;
else if (tdcBin<0) tdcBin += musrAnalysis::rewindTimeBins;
// Int_t tdc = (tdcBin<musrAnalysis::rewindTimeBins) ? tdcBin : tdcBin-musrAnalysis::rewindTimeBins;
dumpFile<<eventNr<<"\t"<<detID<<"\t"<<tdcBin<<"\n";
}
//================================================================
void musrCounter::WriteRewindIntoDumpFile() {
DumpInfoToDumpFile(-2,1000,0);
}
//================================================================

20
musrSimAna/musrCounter.hh
View File

@ -1,5 +1,7 @@
#ifndef musrCounter_h #ifndef musrCounter_h
#define musrCounter_h 1 #define musrCounter_h 1
#include <iostream>
#include <fstream>
#include <TApplication.h> #include <TApplication.h>
#include <TSystem.h> #include <TSystem.h>
#include <TH1.h> #include <TH1.h>
@ -55,17 +57,21 @@ class musrCounter {
void SetTDChistogram(char hName[200],int t0,int t1,int t2,int hNr,char hNameAdd[200]); 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 FillTDChistogram(Double_t variable, Double_t vaha);
void DrawTDChistogram(); 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 RemoveHitsInCounter(Long64_t timeBinLimit);
void RewindHitsInCounter(Long64_t timeBinsToRewind); // 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();
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 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); 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 myPrintThisCounter(Int_t evtID, Int_t detail=2);
// void CheckClockInfo(Long64_t timeBin);
Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;} Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;}
Long64_t GetNumberOfMuonCandidatesAfterVK(){return numberOfMuonCandidatesAfterVK;} Long64_t GetNumberOfMuonCandidatesAfterVK(){return numberOfMuonCandidatesAfterVK;}
Long64_t GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval(){return numberOfMuonCandidatesAfterVKandDoubleHitRemoval;} Long64_t GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval(){return numberOfMuonCandidatesAfterVKandDoubleHitRemoval;}
void DumpInfoToDumpFile(Int_t eventNr, Int_t detID, Long64_t tdcBin);
void WriteRewindIntoDumpFile();
private: private:
// static musrCounter* pointerToAnalysis; // static musrCounter* pointerToAnalysis;
@ -93,10 +99,14 @@ class musrCounter {
Long64_t numberOfMuonCandidates; Long64_t numberOfMuonCandidates;
Long64_t numberOfMuonCandidatesAfterVK; Long64_t numberOfMuonCandidatesAfterVK;
Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval; Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval;
public: public:
static Bool_t bool_ignoreUnperfectMuons; static Bool_t bool_ignoreUnperfectMuons;
static Bool_t bool_ignoreUnperfectPositrons; static Bool_t bool_ignoreUnperfectPositrons;
static Bool_t bool_WriteDataToDumpFile;
static ofstream dumpFile;
// static Long64_t previousClock;
// static Long64_t CLOCK_INTERVAL;
}; };
#endif #endif

5
src/musrDetectorConstruction.cc
View File

@ -960,6 +960,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
sscanf(&line[0],"%*s %*s %*s %lf",&sigma); sscanf(&line[0],"%*s %*s %*s %lf",&sigma);
myMusrScintSD ->SetAPDcellsTimeVariationSigma(sigma*nanosecond); myMusrScintSD ->SetAPDcellsTimeVariationSigma(sigma*nanosecond);
} }
else if (strcmp(varName,"SetAPDcrossTalk")==0) {
double crossTalkProb;
sscanf(&line[0],"%*s %*s %*s %lf",&crossTalkProb);
myMusrScintSD ->SetAPDcrossTalk(crossTalkProb);
}
else { else {
G4cout<<"musrDetectorConstruction.cc: ERROR: Unknown parameterName \"" G4cout<<"musrDetectorConstruction.cc: ERROR: Unknown parameterName \""
<<varName<<"\" ."<<G4endl; <<varName<<"\" ."<<G4endl;

110
src/musrScintSD.cc
View File

@ -88,6 +88,8 @@ musrScintSD::musrScintSD(G4String name)
APDcell_ax =0.3; APDcell_ay=0.3; APDcell_az=0.3; APDcell_ax =0.3; APDcell_ay=0.3; APDcell_az=0.3;
APDcellsEffectRequested=false; APDcellsEffectRequested=false;
APDcellsTimeVariationRequested=false; APDcellsTimeVariationRequested=false;
APDcrossTalkRequested=false;
APDcrossTalkProb=0.;
} }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -147,8 +149,9 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
G4String particleName=aTrack->GetDefinition()->GetParticleName(); G4String particleName=aTrack->GetDefinition()->GetParticleName();
// if (particleName=="opticalphoton") {G4cout<<"UFON JE TU: edep="<<edep<<G4endl; return false;} // if (particleName=="opticalphoton") {G4cout<<"UFON JE TU: edep="<<edep<<G4endl; return false;}
if (particleName=="opticalphoton") { if (particleName=="opticalphoton") {
G4double APDcellsTimeVariation = G4RandGauss::shoot(0,APDcellsTimeVariationSigma); // G4double APDcellsTimeVariation = G4RandGauss::shoot(0,APDcellsTimeVariationSigma);
if (!musrParameters::boolG4OpticalPhotonsUnprocess) ProcessOpticalPhoton(aStep,APDcellsTimeVariation); // if (!musrParameters::boolG4OpticalPhotonsUnprocess) ProcessOpticalPhoton(aStep,APDcellsTimeVariation);
if (!musrParameters::boolG4OpticalPhotonsUnprocess) ProcessOpticalPhoton(aStep);
return false; return false;
} }
@ -203,8 +206,8 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
} }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
void musrScintSD::ProcessOpticalPhoton(G4Step* aStep, G4double APDcellsTimeVariation) { //void musrScintSD::ProcessOpticalPhoton(G4Step* aStep, G4double APDcellsTimeVariation) {
// //Was the photon absorbed by the absorption process ? // //Was the photon absorbed by the absorption process ?
// const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep(); // const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
// G4String processName = (process) ? process->GetProcessName() : "Unknown"; // 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 // 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. // whether a given cell already fired - if so, do not store the photon.
G4bool APDcellAlreadyFired = false; // G4bool APDcellAlreadyFired = false;
G4int APDcellID = 0; G4int APDcellID = (APDcellsEffectRequested) ? FindAPDcellID(aStep) : 0;
if (APDcellsTimeVariationRequested) tmpTime += APDcellsTimeVariation; FireAPDcell(optHitDetectorMap,APDcellID,tmpTime);
// if (APDcellsTimeVariationRequested) tmpTime += G4RandGauss::shoot(0,APDcellsTimeVariationSigma); // if (APDcellsEffectRequested) {
// APDcellID = FindAPDcellID(aStep);
if (APDcellsEffectRequested) { // // G4cout<<"Cell ID="<<APDcellID<<G4endl;
APDcellID = FindAPDcellID(aStep); // for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
// G4cout<<"Cell ID="<<APDcellID<<G4endl; // if ((it2->second)==APDcellID) {
for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) { // // G4cout<<"Already fired cell ="<<APDcellID<<G4endl;
if ((it2->second)==APDcellID) { // APDcellAlreadyFired = true;
// G4cout<<"Already fired cell ="<<APDcellID<<G4endl; // // Now lets check whether the new photon came earlier than the already saved one.
APDcellAlreadyFired = true; // // If so, delete the previous photon and store the time of the new one instead:
// Now lets check whether the new photon came earlier than the already saved one. // if (tmpTime<(it2->first)) {
// If so, delete the previous photon and store the time of the new one instead: // optHitDetectorMap->erase(it2);
if (tmpTime<(it2->first)) { // optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
optHitDetectorMap->erase(it2); // }
optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID)); // break;
} // }
break; // }
} // }
} // if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
} //
if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
// G4cout<<" tmpTime="<<tmpTime<<G4endl; // G4cout<<" tmpTime="<<tmpTime<<G4endl;
// G4cout<<" Detection"<<G4endl; // G4cout<<" Detection"<<G4endl;
// break; // break;
@ -590,6 +592,38 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
} }
if (optHitDetectorMap->empty()) continue; if (optHitDetectorMap->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<G4double,G4int>(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 (rand<APDcrossTalkProb) {
G4int APDcellID = it2->second;
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_START = optHitDetectorMap->begin();
optHitDetectorMapType::const_iterator it2_STOP = it2_START; optHitDetectorMapType::const_iterator it2_STOP = it2_START;
optHitDetectorMapType::const_iterator it2_LAST = optHitDetectorMap->end(); it2_LAST--; optHitDetectorMapType::const_iterator it2_LAST = optHitDetectorMap->end(); it2_LAST--;
@ -908,6 +942,28 @@ G4int musrScintSD::FindAPDcellID(G4Step* aStep) {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time) {
if (!APDcellsEffectRequested) {
optHitDetectorMap->insert(std::pair<G4double,G4int>(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<G4double,G4int>(time,APDcellID));
}
break;
}
}
if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair<G4double,G4int>(time,APDcellID));
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::FindCFDtime(G4double& OPSA_CFD_time, G4double& OPSA_CFD_ampl, G4double timeOfFirstPhoton) { void musrScintSD::FindCFDtime(G4double& OPSA_CFD_time, G4double& OPSA_CFD_ampl, G4double timeOfFirstPhoton) {
OPSA_CFD->Reset("M"); OPSA_CFD->Reset("M");
for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of electronic signal histogram for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of electronic signal histogram