#include "WCSimWCAddDarkNoise.hh" #include "WCSimWCPMT.hh" #include "WCSimWCDigi.hh" #include "WCSimWCHit.hh" #include "WCSimWCTrigger.hh" #include "G4EventManager.hh" #include "G4Event.hh" #include "G4SDManager.hh" #include "G4DigiManager.hh" #include "G4ios.hh" #include "G4RotationMatrix.hh" #include "G4ThreeVector.hh" #include "G4LogicalVolumeStore.hh" #include "WCSimDetectorConstruction.hh" #include "WCSimPmtInfo.hh" #include "WCSimDarkRateMessenger.hh" #include <vector> #include <utility> // for memset #include <cstring> #ifndef WCSIMWCADDDARKNOISE_VERBOSE //#define WCSIMWCADDDARKNOISE_VERBOSE #endif #ifndef NPMTS_VERBOSE #define NPMTS_VERBOSE 10 #endif #ifndef HYPER_VERBOSITY // #define HYPER_VERBOSITY #endif WCSimWCAddDarkNoise::WCSimWCAddDarkNoise(G4String name, WCSimDetectorConstruction* inDetector, G4String detectorElement) :G4VDigitizerModule(name), fCalledAddDarkNoise(false), myDetector(inDetector), detectorElement(detectorElement) { //Set defaults to be unphysical, so that we know if they have been overwritten by the user PMTDarkRate = -99; ConvRate = -99; //Get the user options // DarkRateMessenger = new WCSimDarkRateMessenger(this); DarkRateMessenger = WCSimDarkRateMessenger::GetInstance(); DarkRateMessenger->AddDarkRateInstance(this, detectorElement); DarkRateMessenger->Initialize(); ReInitialize(); } WCSimWCAddDarkNoise::~WCSimWCAddDarkNoise(){ // delete DarkRateMessenger; DarkRateMessenger->RemoveDarkRateInstance(detectorElement); // DarkRateMessenger singleton calls it's destructor when map is empty DarkRateMessenger = 0; } void WCSimWCAddDarkNoise::SetPMTDarkDefaults() { //Grab Dark Rate and Conversion from PMT itself // G4String WCIDCollectionName = myDetector->GetIDCollectionName(); // WCSimPMTObject * PMT = myDetector->GetPMTPointer(WCIDCollectionName); G4String WCCollectionName; if(detectorElement == "tank") WCCollectionName = myDetector->GetIDCollectionName(); else if(detectorElement == "OD") WCCollectionName = myDetector->GetODCollectionName(); else G4cout << "### detectorElement undefined ..." << G4endl; WCSimPMTObject * PMT = myDetector->GetPMTPointer(WCCollectionName); double const conversion_to_kHz = 1000000; //ToDo: remove this and treat DarkRate in CLHEP units throughout the class. double defaultPMTDarkRate = PMT->GetDarkRate() * conversion_to_kHz; double defaultConvRate = PMT->GetDarkRateConversionFactor(); //Only set the defaults if the user hasn't overwritten the unphysical defaults if(PMTDarkRate < -98) PMTDarkRate = defaultPMTDarkRate; if(ConvRate < -98) ConvRate = defaultConvRate; } void WCSimWCAddDarkNoise::AddDarkNoise(){ //Grab the PMT-specific defaults if(!fCalledAddDarkNoise) { SetPMTDarkDefaults(); fCalledAddDarkNoise = true; } //clear the result and range vectors ReInitialize(); G4DigiManager* DigiMan = G4DigiManager::GetDMpointer(); // Get the PMT collection ID // G4int WCHCID = DigiMan->GetDigiCollectionID("WCRawPMTSignalCollection"); // Get the PMT collection ID G4String thecollectionName; if(detectorElement=="tank") thecollectionName="WCRawPMTSignalCollection"; else if(detectorElement=="OD") thecollectionName="WCRawPMTSignalCollection_OD"; else G4cout << "### detectorElement undefined ..." << G4endl; G4int WCHCID = DigiMan->GetDigiCollectionID(thecollectionName); // Get the PMT Digits collection WCSimWCDigitsCollection* WCHCPMT = (WCSimWCDigitsCollection*)(DigiMan->GetDigiCollection(WCHCID)); #ifdef HYPER_VERBOSITY if(detectorElement=="OD"){G4cout<<"WCSimWCAddDarkNoise::AddDarkNoise ☆ retrieved hit collection (WCSimWCDigitsCollection*)"<<thecollectionName<<" which has "<<WCHCPMT->entries()<<" entries"<<G4endl;} #endif G4String thetriggertype=""; if(detectorElement=="OD") { // check to see if this detector element uses the tank for triggering WCSimWCTriggerBase *WCTM = (WCSimWCTriggerBase *) DigiMan->FindDigitizerModule("WCReadout_OD"); thetriggertype = WCTM->GetTriggerClassName(); } WCSimWCDigitsCollection* WCHCPMT_tank=NULL; if(thetriggertype=="TankDigits"){ //if triggertype is TankDigits, then we should add noise in the windows around the *tank* digits, not this //collection's digits, because those are the windows that will be relevant when reading out triggers. int WCHCID_tank = DigiMan->GetDigiCollectionID("WCRawPMTSignalCollection"); WCHCPMT_tank = (WCSimWCDigitsCollection*)(DigiMan->GetDigiCollection(WCHCID_tank)); (WCSimWCDigitsCollection*)(DigiMan->GetDigiCollection(WCHCID)); #ifdef HYPER_VERBOSITY if(detectorElement=="OD"){ G4cout<<"WCSimWCAddDarkNoise::AddDarkNoise ☆ retrieved hit collection (WCSimWCDigitsCollection*)WCRawPMTSignalCollection for finding dark noise windows, which has "<<WCHCPMT_tank->entries()<<" entries"<<G4endl;} #endif } // if ((WCHCPMT != NULL) && (this->PMTDarkRate > 1E-307)) { if (((WCHCPMT != NULL)||(thetriggertype=="TankDigits"&&WCHCPMT_tank!=NULL)) && (this->PMTDarkRate > 1E-307)) { //Determine ranges for adding noise // if(DarkMode == 1) // FindDarkNoiseRanges(WCHCPMT,this->DarkWindow); if(DarkMode == 1){ if(thetriggertype=="TankDigits"){ WCSimWCAddDarkNoise* WCDNM_tank = (WCSimWCAddDarkNoise*)DigiMan->FindDigitizerModule("WCDarkNoise"); int DarkWindow_tank = WCDNM_tank->GetDarkWindow(); #ifdef HYPER_VERBOSITY if(detectorElement=="OD") G4cout<<"WCSimWCAddDarkNoise::AddDarkNoise ☆ calling FindDarkNoiseRanges() with tank raw collection and darkwindow size"<<G4endl; #endif FindDarkNoiseRanges(WCHCPMT_tank,DarkWindow_tank); } else { FindDarkNoiseRanges(WCHCPMT,this->DarkWindow); } } else if(DarkMode == 0) { result.push_back(std::pair<double,double>(DarkLow,DarkHigh)); } //Call routine to add dark noise here. //loop over pairs which represent ranges. //Add noise to those ranges #ifdef HYPER_VERBOSITY if(detectorElement=="OD") G4cout<<"WCSimWCAddDarkNoise::AddDarkNoise ☆ adding dark noise hits in "<<result.size()<<" window(s) around found digits."<<G4endl; #endif int windowfordarknoise=0; for(std::vector<std::pair<double, double> >::iterator it2 = result.begin(); it2 != result.end(); it2++) { #ifdef HYPER_VERBOSITY if(detectorElement=="OD") G4cout<<"WCSimWCAddDarkNoise::AddDarkNoise ☆ adding dark noise in window "<<windowfordarknoise<<G4endl; windowfordarknoise++; #endif AddDarkNoiseBeforeDigi(WCHCPMT,it2->first,it2->second); } } } void WCSimWCAddDarkNoise::AddDarkNoiseBeforeDigi(WCSimWCDigitsCollection* WCHCPMT, double num1 ,double num2) { // Introduces dark noise into each PMT during an event window // This won't introduce noise only events, and isn't written // to handle different rates for each PMT (although this shouldn't // be too difficult to add at a later time) // // Added by: Morgan Askins (maskins@ucdavis.edu) G4int number_entries = WCHCPMT->entries(); // const G4int number_pmts = myDetector->GetTotalNumPmts(); #ifdef HYPER_VERBOSITY if(detectorElement=="OD")G4cout<<"WCSimWCAddDarkNoise::AddDarkNoiseBeforeDigi ☆ adding dark noise to collection of "<<number_entries<<" entries"<<G4endl; #endif G4int thenum_pmts; if(detectorElement=="tank") thenum_pmts = myDetector->GetTotalNumPmts(); else if(detectorElement=="OD") thenum_pmts = myDetector->GetTotalNumODPmts(); else G4cout << "### detectorElement undefined ..." << G4endl; const G4int number_pmts=thenum_pmts; int *PMTindex = new int [number_pmts+1]; //initialize PMTindex for (int l=0; l<number_pmts+1; l++){ PMTindex[l] =0; } // G4cout<<"entries before "<<WCHCPMT->entries()<<"\n"; //Set up proper indices for tubes which have already been hit int num_hit_b4=0; for (int g=0; g<number_entries; g++){ G4int tube = (*WCHCPMT)[g]->GetTubeID(); //G4cout<<"totalpe "<<tube<<" "<<(*WCHCPMT)[g]->GetTotalPe()<<"\n"; //should this be tube-1? PMTindex[tube] += (*WCHCPMT)[g]->GetTotalPe(); num_hit_b4 += (*WCHCPMT)[g]->GetTotalPe(); //G4cout<<"TotalPe "<<(*WCHCPMT)[g]->GetTotalPe()<<" "<<PMTindex[tube]<<"\n"; } // Get the info for pmt positions // std::vector<WCSimPmtInfo*> *pmts = myDetector->Get_Pmts(); std::vector<WCSimPmtInfo*> *pmts; if(detectorElement=="tank") pmts = myDetector->Get_Pmts(); else if(detectorElement=="OD")pmts = myDetector->Get_ODPmts(); else G4cout << "### detectorElement undefined ..." << G4endl; // It works out that the pmts here are ordered ! // pmts->at(i) has tubeid i+1 std::vector<int> list; list.assign( number_pmts+1, 0 ); for( int h = 0; h < number_entries; h++ ) { list[(*WCHCPMT)[h]->GetTubeID()] = h+1; } // Add noise to PMT's here, do so in the range num1 to num2 double current_time = 0; double pe = 0.0; //Calculate the time window size double windowsize = num2 - num1; G4DigiManager* DMman = G4DigiManager::GetDMpointer(); // WCSimWCPMT* WCPMT = (WCSimWCPMT*)DMman->FindDigitizerModule("WCReadoutPMT"); G4String thewcpmtname; if(detectorElement=="tank") thewcpmtname="WCReadoutPMT"; else if(detectorElement=="OD") thewcpmtname="WCReadoutPMT_OD"; else G4cout << "### detectorElement undefined ..." << G4endl; WCSimWCPMT* WCPMT = (WCSimWCPMT*)DMman->FindDigitizerModule(thewcpmtname); #ifdef HYPER_VERBOSITY if(detectorElement=="OD"){G4cout<<"WCSimWCAddDarkNoise::AddDarkNoiseBeforeDigi ☆ getting (WCSimWCPMT*)"<<thewcpmtname<<G4endl;} #endif //average number of PMTs with noise double ave=number_pmts * this->PMTDarkRate * this->ConvRate * windowsize * 1E-6; //poisson distributed noise, number of noise hits to add int nnoispmt = CLHEP::RandPoisson::shoot(ave); G4String thecollectionName; if(detectorElement=="tank") thecollectionName="WCRawPMTSignalCollection"; else if(detectorElement=="OD") thecollectionName="WCRawPMTSignalCollection_OD"; else G4cout << "### detectorElement undefined ..." << G4endl; #ifdef HYPER_VERBOSITY if(detectorElement=="OD"){G4cout<<"WCSimWCAddDarkNoise::AddDarkNoiseBeforeDigi ☆ adding "<<nnoispmt<<" dark hits to "<<thecollectionName<<G4endl;} #endif #ifdef WCSIMWCADDDARKNOISE_VERBOSE G4cout << "WCSimWCAddDarkNoise::AddDarkNoiseBeforeDigi Going to add " << nnoispmt << " dark noise hits in time window [" << num1 << "," << num2 << "] duration " << num2 - num1 << G4endl; #endif for( int i = 0; i < nnoispmt; i++ ) { //time of noise hit to be generated //A time from t=num1 to num2 current_time = num1 + G4UniformRand()*windowsize; //now a random PMT. Assuming noise levels are the same for //each PMT. int noise_pmt = static_cast<int>( G4UniformRand() * number_pmts ) + 1; //so that pmt numbers runs from 1 to Npmt if( list[ noise_pmt ] == 0 ) { //PMT has no hits yet. Create a new WCSimWCDigi WCSimWCDigi* ahit = new WCSimWCDigi(); ahit->SetTubeID( noise_pmt); //G4cout<<"setting new noise pmt "<<noise_pmt<<" "<<ahit->GetTubeID()<<"\n"; // This Logical volume is GlassFaceWCPMT // ahit->SetLogicalVolume(G4LogicalVolumeStore::GetInstance()->GetVolume(myDetector->GetDetectorName()+"-glassFaceWCPMT")); if(detectorElement=="tank") ahit->SetLogicalVolume(G4LogicalVolumeStore::GetInstance()->GetVolume(myDetector->GetDetectorName()+"-glassFaceWCPMT")); else if(detectorElement=="OD")ahit->SetLogicalVolume(G4LogicalVolumeStore::GetInstance()->GetVolume(myDetector->GetDetectorName()+"-glassFaceWCPMT_OD")); else G4cout << "### detectorElement undefined ..." << G4endl; //G4cout<<"1 "<<(G4LogicalVolumeStore::GetInstance()->GetVolume("glassFaceWCPMT"))->GetName()<<"\n"; //G4cout<<"2 "<<(*WCHCPMT)[0]->GetLogicalVolume()->GetName()<<"\n"; ahit->SetTrackID(-1); ahit->SetParentID(PMTindex[noise_pmt], -1); // Set the position and rotation of the pmt Double_t hit_pos[3]; Double_t hit_rot[3]; // TODO: need to change the format of hit_pos to G4ThreeVector // and change hit_rot to G4RotationMatrix WCSimPmtInfo* pmtinfo = (WCSimPmtInfo*)pmts->at( noise_pmt -1 ); hit_pos[0] = 10*pmtinfo->Get_transx(); hit_pos[1] = 10*pmtinfo->Get_transy(); hit_pos[2] = 10*pmtinfo->Get_transz(); hit_rot[0] = pmtinfo->Get_orienx(); hit_rot[1] = pmtinfo->Get_orieny(); hit_rot[2] = pmtinfo->Get_orienz(); G4RotationMatrix pmt_rotation(hit_rot[0], hit_rot[1], hit_rot[2]); G4ThreeVector pmt_position(hit_pos[0], hit_pos[1], hit_pos[2]); ahit->SetRot(pmt_rotation); ahit->SetPos(pmt_position); ahit->SetTime(PMTindex[noise_pmt],current_time); ahit->SetPreSmearTime(PMTindex[noise_pmt],current_time); //presmear==postsmear for dark noise pe = WCPMT->rn1pe(); ahit->SetPe(PMTindex[noise_pmt],pe); //Added this line to increase the totalPe by 1 ahit->AddPe(current_time); WCHCPMT->insert(ahit); PMTindex[noise_pmt]++; number_entries ++; list[ noise_pmt ] = number_entries; // Add this PMT to the end of the list #ifdef WCSIMWCADDDARKNOISE_VERBOSE if(noise_pmt < NPMTS_VERBOSE) G4cout << "WCSimWCAddDarkNoise::AddDarkNoiseBeforeDigi Added NEW DIGI with dark noise hit at time " << current_time << " to PMT " << noise_pmt << G4endl; #endif } else { (*WCHCPMT)[ list[noise_pmt]-1 ]->AddPe(current_time); pe = WCPMT->rn1pe(); (*WCHCPMT)[ list[noise_pmt]-1 ]->SetPe(PMTindex[noise_pmt],pe); (*WCHCPMT)[ list[noise_pmt]-1 ]->SetTime(PMTindex[noise_pmt],current_time); (*WCHCPMT)[ list[noise_pmt]-1 ]->SetPreSmearTime(PMTindex[noise_pmt],current_time); //presmear==postsmear for dark noise (*WCHCPMT)[ list[noise_pmt]-1 ]->SetParentID(PMTindex[noise_pmt],-1); PMTindex[noise_pmt]++; #ifdef WCSIMWCADDDARKNOISE_VERBOSE if(noise_pmt < NPMTS_VERBOSE) G4cout << "WCSimWCAddDarkNoise::AddDarkNoiseBeforeDigi Added to exisiting digi a dark noise hit at time " << current_time << " to PMT " << noise_pmt << G4endl; #endif } }//i (number of noise hits to add) delete [] PMTindex; return; } void WCSimWCAddDarkNoise::FindDarkNoiseRanges(WCSimWCDigitsCollection* WCHCPMT, double width) { //Loop over all Hits and assign a time window around each hit //store these in the ranges vector as pairs for (int g=0; g<WCHCPMT->entries(); g++){ for (int gp=0; gp<(*WCHCPMT)[g]->GetTotalPe(); gp++){ double time = (*WCHCPMT)[g]->GetTime(gp); //lets assume a 5us window. So we centre this on the hit time. //t1 is the lower limit of the window. double t1=time - width/2.; double t2=time + width/2.; ranges.push_back(std::pair<double, double>(t1, t2)); } } #ifdef HYPER_VERBOSITY if(detectorElement=="mrd") G4cout<<"WCSimWCAddDarkNoise::FindDarkNoiseRanges ☆ "<<ranges.size()<<" windows around digits found before pruning."<<G4endl; #endif //we need to ensure that the ranges found above are sorted first //for the algorithm below to work sort(ranges.begin(),ranges.end()); //check if the vector range has any entries //If no entries this indicates that no digits were //found in the WCSimWCDigitsCollection, which can cause //segmentation faults in the next part of the code in this method //which searches for overlapping ranges. //Set range vector to have one element from 0 to 0 (so no noise digits will be added) if(ranges.size() == 0) { //push back a range of 0 and 0 and return result.push_back(std::make_pair(0.,0.)); return; } //the ranges vector contains overlapping ranges //this loop removes overlaps //output are pairs stored in the result vector std::vector<std::pair<double, double> >::iterator it = ranges.begin(); std::pair<double, double> current = *(it)++; for( ; it != ranges.end(); it++) { if (current.second >= it->first){ current.second = std::max(current.second, it->second); } else { result.push_back(current); current = *(it); } } result.push_back(current); //now we should have a vector of non-overlapping range pairs to pass to the //dark noise routine } void WCSimWCAddDarkNoise::SaveOptionsToOutput(WCSimRootOptions * wcopt, string tag) { wcopt->SetPMTDarkRate(tag, PMTDarkRate); wcopt->SetConvRate(tag, ConvRate); wcopt->SetDarkHigh(tag, DarkHigh); wcopt->SetDarkLow(tag, DarkLow); wcopt->SetDarkWindow(tag, DarkWindow); wcopt->SetDarkMode(tag, DarkMode); }