#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 #include // for memset #include #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*)"<entries()<<" entries"<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 "<entries()<<" entries"<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"<DarkWindow); } } else if(DarkMode == 0) { result.push_back(std::pair(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 "< >::iterator it2 = result.begin(); it2 != result.end(); it2++) { #ifdef HYPER_VERBOSITY if(detectorElement=="OD") G4cout<<"WCSimWCAddDarkNoise::AddDarkNoise ☆ adding dark noise in window "<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 "<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; lentries()<<"\n"; //Set up proper indices for tubes which have already been hit int num_hit_b4=0; for (int g=0; gGetTubeID(); //G4cout<<"totalpe "<GetTotalPe()<<"\n"; //should this be tube-1? PMTindex[tube] += (*WCHCPMT)[g]->GetTotalPe(); num_hit_b4 += (*WCHCPMT)[g]->GetTotalPe(); //G4cout<<"TotalPe "<<(*WCHCPMT)[g]->GetTotalPe()<<" "< *pmts = myDetector->Get_Pmts(); std::vector *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 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*)"<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 "<( 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 "<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; gentries(); 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(t1, t2)); } } #ifdef HYPER_VERBOSITY if(detectorElement=="mrd") G4cout<<"WCSimWCAddDarkNoise::FindDarkNoiseRanges ☆ "< >::iterator it = ranges.begin(); std::pair 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); }