/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus * * MAUS is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * MAUS is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MAUS. If not, see . * */ #include "src/map/MapCppEMRMCDigitization/MapCppEMRMCDigitization.hh" namespace MAUS { PyMODINIT_FUNC init_MapCppEMRMCDigitization(void) { PyWrapMapBase::PyWrapMapBaseModInit ("MapCppEMRMCDigitization", "", "", "", ""); } MapCppEMRMCDigitization::MapCppEMRMCDigitization() : MapBase("MapCppEMRMCDigitization") { } void MapCppEMRMCDigitization::_birth(const std::string& argJsonConfigDocument) { _classname = "MapCppEMRMCDigitization"; char* pMAUS_ROOT_DIR = getenv("MAUS_ROOT_DIR"); if (!pMAUS_ROOT_DIR) { throw MAUS::Exception(Exception::recoverable, "Could not resolve ${MAUS_ROOT_DIR} environment variable", "MapCppEMRMCDigitization::birth"); } // JsonCpp setup Json::Value configJSON; configJSON = JsonWrapper::StringToJson(argJsonConfigDocument); // Fetch variables _number_of_planes = configJSON["EMRnumberOfPlanes"].asInt(); _number_of_bars = configJSON["EMRnumberOfBars"].asInt(); _secondary_n_low = configJSON["EMRsecondaryNLow"].asInt(); _secondary_tot_low = configJSON["EMRsecondaryTotLow"].asInt(); _secondary_bunching_width = configJSON["EMRsecondaryBunchingWidth"].asInt(); _do_sampling = configJSON["EMRdoSampling"].asInt(); _nph_per_MeV = configJSON["EMRnphPerMeV"].asInt(); _seed = configJSON["EMRseed"].asInt(); _trap_eff = configJSON["EMRtrapEff"].asDouble(); _QE_MAPMT = configJSON["EMRqeMAPMT"].asDouble(); _QE_SAPMT = configJSON["EMRqeSAPMT"].asDouble(); _nADC_per_pe_MAPMT = configJSON["EMRnadcPerPeMAPMT"].asDouble(); _nADC_per_pe_SAPMT = configJSON["EMRnadcPerPeSAPMT"].asDouble(); _electronics_response_spread_MAPMT = configJSON["EMRelectronicsResponseSpreadMAPMT"].asDouble(); _electronics_response_spread_SAPMT = configJSON["EMRelectronicsResponseSpreadSAPMT"].asDouble(); _dbb_count = configJSON["EMRdbbCount"].asDouble(); _fadc_count = configJSON["EMRfadcCount"].asDouble(); _time_response_spread = configJSON["EMRtimeResponseSpread"].asInt(); _tot_func_p1 = configJSON["EMRtotFuncP1"].asDouble(); _tot_func_p2 = configJSON["EMRtotFuncP2"].asDouble(); _tot_func_p3 = configJSON["EMRtotFuncP3"].asDouble(); _acquisition_window = configJSON["EMRacquisitionWindow"].asInt(); _signal_integration_window = configJSON["EMRsignalIntegrationWindow"].asInt(); _arrival_time_shift = configJSON["EMRarrivalTimeShift"].asInt(); _arrival_time_gaus_width = configJSON["EMRarrivalTimeGausWidth"].asDouble(); _arrival_time_uniform_width = configJSON["EMRarrivalTimeUniformWidth"].asDouble(); _pulse_shape_landau_width = configJSON["EMRpulseShapeLandauWidth"].asDouble(); _fom = configJSON["EMRfom"].asString(); _birks_constant = configJSON["EMRbirksConstant"].asDouble(); _signal_energy_threshold = configJSON["EMRsignalEnergyThreshold"].asDouble(); _baseline_spread = configJSON["EMRbaselineSpread"].asInt(); _maximum_noise_level = configJSON["EMRmaximumNoiseLevel"].asInt(); _deltat_shift = configJSON["EMRdeltatShift"].asDouble(); _baseline_position = configJSON["EMRbaselinePosition"].asDouble(); _arrival_time_spread = configJSON["EMRarrivalTimeSpread"].asDouble(); // Generate random pedestals and noise in the SAPMTs _rand = new TRandom3(_seed); for (size_t iPlane = 0; iPlane < _number_of_planes; iPlane++) { _baseline[iPlane] = _baseline_position + static_cast(_rand->Uniform(-_baseline_spread, _baseline_spread)); _noise_level[iPlane] = static_cast(_rand->Uniform(0, _maximum_noise_level)); _noise_position[iPlane] = static_cast(_rand->Uniform(0, 1.9999999)); } // Load the EMR calibration map bool loaded = _calibMap.InitializeFromCards(configJSON); if (!loaded) throw(Exception(Exception::recoverable, "Could not find EMR calibration map", "MapCppEMRMCDigitizer::birth")); // Load the EMR attenuation map loaded = _attenMap.InitializeFromCards(configJSON); if (!loaded) throw(Exception(Exception::recoverable, "Could not find EMR attenuation map", "MapCppEMRMCDigitizer::birth")); } void MapCppEMRMCDigitization::_death() { } void MapCppEMRMCDigitization::_process(Data *data) const { // Routine data checks before processing it if ( !data ) throw Exception(Exception::recoverable, "Data was NULL", "MapCppEMRMCDigitization::_process"); Spill* spill = data->GetSpill(); if ( !spill ) throw Exception(Exception::recoverable, "Spill was NULL", "MapCppEMRMCDigitization::_process"); if ( spill->GetDaqEventType() != "physics_event" ) return; MCEventPArray* mcEvts = spill->GetMCEvents(); if ( !mcEvts ) throw Exception(Exception::recoverable, "MCEventPArray was NULL", "MapCppEMRMCDigitization::_process"); int nPartEvents = spill->GetMCEventSize(); if ( !nPartEvents ) return; // Check the Recon event array and EMR Spill Data are initialised. If not make it so if ( !spill->GetReconEvents() ) spill->SetReconEvents(new ReconEventPArray()); if ( !spill->GetEMRSpillData() ) spill->SetEMRSpillData(new EMRSpillData()); // Set primary time window std::vector delta_t_array; delta_t_array.resize(0); for (int iPe = 0; iPe < nPartEvents; iPe++) { EMRHitArray *hits = spill->GetAnMCEvent(iPe).GetEMRHits(); if ( hits ) { Primary *primary = spill->GetAnMCEvent(iPe).GetPrimary(); double pTime = primary->GetTime(); // ns for (size_t iHit = 0; iHit < hits->size(); iHit++) { EMRHit hit = hits->at(iHit); double delta_t = hit.GetTime() - pTime; if (delta_t > 0) delta_t_array.push_back(delta_t); } } } // Procede if there's no EMR data if ( !delta_t_array.size() ) return; // Set the window from the leading time to 20 ADC counts later int lt = static_cast (*std::min_element(delta_t_array.begin(), delta_t_array.end())/_dbb_count); int deltat_min = lt; // ns int deltat_max = lt + 50; // ns, +50 ns // Create a temporary array of n+1 BarHitTmp arrays (1 per trigger + decays) EMREventPlaneHitVector mc_events_tmp = get_emr_events_tmp(nPartEvents+1); // EMREvents // Fill the temporary array with G4 spill information processMC(mcEvts, mc_events_tmp, deltat_min, deltat_max); // Create a temporary array of n plane hit arrays (1 per trigger) and a temporary array of bars EMREventPlaneHitVector emr_events_tmp = get_emr_events_tmp(nPartEvents); // EMREvents EMRBarHitArray emr_bar_hits_tmp; // EMRSpillData // Digitize the G4 spill digits digitize(mc_events_tmp, emr_events_tmp, emr_bar_hits_tmp, nPartEvents); // Create a temporary array to harbour n' plane hit arrays (1 per daughter candiate) EMREventPlaneHitVector emr_candidates_tmp; // Process the secondary hits and merge them timewise into daughter candidates process_candidates(emr_bar_hits_tmp, emr_candidates_tmp); emr_events_tmp.insert(emr_events_tmp.end(), emr_candidates_tmp.begin(), emr_candidates_tmp.end()); // Fill the Recon event array with spill information (1 per trigger + daughter candidates) fill(spill, emr_events_tmp, nPartEvents); } void MapCppEMRMCDigitization::processMC(MAUS::MCEventPArray *mcEvts, EMREventPlaneHitVector& mc_events_tmp, int deltat_min, int deltat_max) const { size_t nPartEvents = mcEvts->size(); for (size_t iPe = 0; iPe < nPartEvents; iPe++) { Primary *primary = mcEvts->at(iPe)->GetPrimary(); double pTime = primary->GetTime(); // ns EMRHitArray *EMRhits = mcEvts->at(iPe)->GetEMRHits(); for (size_t iHit = 0; iHit < EMRhits->size(); iHit++) { EMRHit hit = EMRhits->at(iHit); // Channel ID EMRChannelId *xCh = hit.GetChannelId(); int g4barid = xCh->GetBar(); int xPlane = g4barid / 59; int xBar = g4barid % 59 + 1; // Position ThreeVector xPos = hit.GetPosition(); // (mm, mm, mm) double xPathLength = hit.GetPathLength(); // mm // Time double time = hit.GetTime(); // ns int xTime = static_cast(time); // ns, leading time int xDeltaT = xTime - pTime; // ns // Energy deposition double edep = hit.GetEnergyDeposited(); // MeV int xTot = static_cast(edep*pow(10, 6)); // eV, non digitized time-over-threshold // Fill the bar hit data (from the EMR sensitive detector) EMRBarHitTmp bHitTmp; EMRBarHit bHit; bHit.SetChannel(xPlane*_number_of_bars+xBar); bHit.SetTot(xTot); bHit.SetTime(xTime); bHit.SetDeltaT(xDeltaT); bHitTmp._barhit = bHit; bHitTmp._pos = xPos; bHitTmp._path_length = xPathLength; if ( xDeltaT >= deltat_min && xDeltaT < deltat_max ) { mc_events_tmp[iPe][xPlane]._barhits.push_back(bHitTmp); mc_events_tmp[iPe][xPlane]._orientation = xPlane % 2; } else { mc_events_tmp[nPartEvents][xPlane]._barhits.push_back(bHitTmp); mc_events_tmp[nPartEvents][xPlane]._orientation = xPlane % 2; } } } } void MapCppEMRMCDigitization::digitize(EMREventPlaneHitVector mc_events_tmp, EMREventPlaneHitVector& emr_events_tmp, EMRBarHitArray& emr_bar_hits_tmp, size_t nPartEvents) const { for (size_t iPe = 0; iPe < mc_events_tmp.size(); iPe++) { for (size_t iPlane = 0; iPlane < _number_of_planes; iPlane++) { if ( !mc_events_tmp[iPe][iPlane]._barhits.size() ) continue; // Find PMT gain correction coefficient EMRChannelKey xAverageKey(iPlane, iPlane%2, -1, "emr"); double epsilon_MA = _calibMap.Eps(xAverageKey, "MA"); double epsilon_SA = _calibMap.Eps(xAverageKey, "SA"); int naph_SAPMT(0); // Number of attenuated photons reaching the SAPMT for (size_t iHit = 0; iHit < mc_events_tmp[iPe][iPlane]._barhits.size(); iHit++) { EMRBarHitTmp bHitTmp = mc_events_tmp[iPe][iPlane]._barhits[iHit]; EMRBarHit bHit = bHitTmp._barhit; int xBar = bHit.GetChannel()%_number_of_bars; EMRChannelKey xKey(iPlane, iPlane%2, xBar, "emr"); int xTot = bHit.GetTot(); // eV int xTime = bHit.GetTime(); // ns int xDeltaT = bHit.GetDeltaT(); // ns double x = bHitTmp._pos.x(); // mm double y = bHitTmp._pos.y(); // mm double xPathLength = bHitTmp._path_length; // mm //--------- MAPMT signal simulation -------// // Simulate electronic noise and MAPMT dark current // !!! TODO !!! // Geant4 information: edep, time and path length double energy = static_cast(xTot)/pow(10, 6); // MeV int time = xDeltaT; // ns if ( energy < _signal_energy_threshold ) continue; // Convert energy to the number of scintillation photons (nsph) according to Birks's law int nsph = static_cast ((_nph_per_MeV*energy)/(1.0+_birks_constant*energy/xPathLength)); if ( _do_sampling ) nsph = _rand->Poisson(nsph); // Convert nsph to the number of trapped photons (ntph) int ntph = static_cast(nsph*_trap_eff); if ( _do_sampling ) ntph = _rand->Poisson(ntph); // Attenuate ntph according to fibre length (naph) and loss in the connectors // MAPMT (half the PEs) int naph_MAPMT = static_cast (static_cast(ntph)/2 *_attenMap.fibreAtten(xKey, x, y, "MA") *_attenMap.connectorAtten(xKey, "MA")); if ( _do_sampling ) naph_MAPMT = _rand->Poisson(naph_MAPMT); // SAPMT (other half of the PEs) int naph_SAPMT_hit = static_cast (static_cast(ntph)/2 *_attenMap.fibreAtten(xKey, x, y, "SA") *_attenMap.connectorAtten(xKey, "SA")); naph_SAPMT += naph_SAPMT_hit; // Simulate crosstalk and misalignment // !!! TODO !!! // Convert naph to the number of photoelectrons (npe) int npe = static_cast(naph_MAPMT*_QE_MAPMT); if ( _do_sampling ) npe = _rand->Poisson(npe); // Correct npe for the photocathode non-uniformity // !!! TODO !!! // Correct npe for the gain difference between MAPMTs npe = static_cast(npe*epsilon_MA); if ( !npe ) continue; // Convert npe to the number of ADC counts int nADC = static_cast(npe*_nADC_per_pe_MAPMT); if ( _do_sampling ) nADC = static_cast(_rand->Gaus(nADC, _electronics_response_spread_MAPMT)); // Convert nADC to a time-over-threshold int xTotDigi = static_cast (_tot_func_p1*log(_tot_func_p2*nADC+_tot_func_p3)); if (xTotDigi <= 0) continue; // Convert Geant4 hit time to deltaT in ADC counts int xDeltaTDigi = static_cast(static_cast(time)/_dbb_count); if ( _do_sampling ) xDeltaTDigi = static_cast(_rand->Gaus(xDeltaTDigi, _time_response_spread)); if (xDeltaTDigi < 0) xDeltaTDigi = 0; // Correct deltaT with fibre length xDeltaTDigi += static_cast (_attenMap.fibreDelay(xKey, x, y, "MA")/_dbb_count); // Set hit time int xTimeDigi = static_cast(static_cast(xTime)/_dbb_count); // Set bar hit EMRBarHitTmp barHitTmp; EMRBarHit barHit; barHit.SetChannel(bHit.GetChannel()); barHit.SetTot(xTotDigi); barHit.SetTime(xTimeDigi); // Discriminate primary hits (close to the trigger) from the rest if ( iPe < nPartEvents ) { barHit.SetDeltaT(xDeltaTDigi); barHitTmp._barhit = barHit; emr_events_tmp[iPe][iPlane]._barhits.push_back(barHitTmp); } else { barHit.SetDeltaT(0); emr_bar_hits_tmp.push_back(barHit); } } //------------ SAPMT signal simulation -------// if ( iPe == nPartEvents ) continue; int xOri = iPlane % 2; // Simulate SAPMT dark current // !!! TODO !!! // Sample naph with Poisson if ( _do_sampling ) naph_SAPMT = _rand->Poisson(naph_SAPMT); // Convert naph to the number of photoelectrons (npe) int npe = static_cast(naph_SAPMT*_QE_SAPMT); if ( _do_sampling ) npe = _rand->Poisson(npe); // Correct npe for the gain difference between SAPMTs npe = static_cast(npe*epsilon_SA); // Convert npe to the number of ADC counts int nADC = static_cast(npe*_nADC_per_pe_SAPMT); if ( _do_sampling ) nADC = static_cast(_rand->Gaus(nADC, _electronics_response_spread_SAPMT)); if (nADC < 0) nADC = 0; // Set signal baseline (8bit ADC) int baseline = _baseline[iPlane]; // Set electronics noise level - number of fluctuations within acquisition window int noise_level = _noise_level[iPlane]; // Set electronics noise position (upwards/downwards fluctuations) int noise_position = _noise_position[iPlane]; // Generate negative voltage pulse with random electronics noise TH1F *pulse_shape = new TH1F("pulse_shape", "pulse_shape", _acquisition_window, 0, _acquisition_window); double arrival_time = _signal_integration_window*2 + _arrival_time_shift; if ( static_cast(_rand->Uniform(0, 1.9999999)) ) arrival_time += _arrival_time_spread; arrival_time += _rand->Gaus(0, _arrival_time_gaus_width) +_rand->Uniform(-_arrival_time_uniform_width, _arrival_time_uniform_width); for (int i = 0; i < noise_level; i++) { double noise = _rand->Uniform(0, _acquisition_window); int bin = static_cast(noise+1); if ( pulse_shape->GetBinContent(bin) != 1 ) pulse_shape->Fill(noise); } if ( noise_position ) for (Int_t i = 1; i <= _acquisition_window; i++) pulse_shape->SetBinContent(i, -pulse_shape->GetBinContent(i)); for (int i = 0; i < nADC; i++) pulse_shape->Fill(_rand->Landau(arrival_time, _pulse_shape_landau_width)); for (int i = 1; i <= _acquisition_window; i++) pulse_shape->SetBinContent(i, baseline-pulse_shape->GetBinContent(i)); double pedestal_baseline = pulse_shape->Integral(1, _signal_integration_window) /_signal_integration_window; // Set pedestal area double xPedestalAreaDigi = pedestal_baseline*_signal_integration_window - pulse_shape->Integral(_signal_integration_window+1, _signal_integration_window*2); // Set pulse area double xAreaDigi = pedestal_baseline*_signal_integration_window - pulse_shape->Integral(arrival_time-10+1, arrival_time+_signal_integration_window-10); if (nADC == 0) xAreaDigi = 0; if (xAreaDigi < 0) xAreaDigi = 0; // Set hit arrival time int xArrivalTimeDigi = arrival_time; // Simulate time delay in cables // !!! TODO !!! // Sample pulse shape every 2ns (500 MHz clock) std::vector xSamplesDigi; for (int iBin = 0; iBin < _acquisition_window; iBin++) xSamplesDigi.push_back(pulse_shape->GetBinContent(iBin+1)); // Fill the fADC information emr_events_tmp[iPe][iPlane]._orientation = xOri; emr_events_tmp[iPe][iPlane]._charge = xAreaDigi; // ADC emr_events_tmp[iPe][iPlane]._pedestal_area = xPedestalAreaDigi; // ADC emr_events_tmp[iPe][iPlane]._time = xArrivalTimeDigi; // ADC emr_events_tmp[iPe][iPlane]._samples = xSamplesDigi; // ADC delete pulse_shape; } } } void MapCppEMRMCDigitization::process_candidates(EMRBarHitArray emr_bar_hits_tmp, EMREventPlaneHitVector& emr_candidates_tmp) const { // Sort the vector in ascending order of hit time (lambda sorting function) sort(emr_bar_hits_tmp.begin(), emr_bar_hits_tmp.end(), [] (const EMRBarHit& a, const EMRBarHit& b) { return a.GetTime() < b.GetTime(); }); // Find groups of hits defined by mat std::vector barHitGroupVector; EMRBarHitArray barHitGroup; int xHitTime = 0; for (size_t iHit = 1; iHit < emr_bar_hits_tmp.size()+1; iHit++) { int xDeltaT = _secondary_bunching_width + 1; if ( !xHitTime ) xHitTime = emr_bar_hits_tmp[iHit-1].GetTime(); if (iHit < emr_bar_hits_tmp.size()) xDeltaT = emr_bar_hits_tmp[iHit].GetTime() - xHitTime; if (xDeltaT > _secondary_bunching_width) { barHitGroup.push_back(emr_bar_hits_tmp[iHit-1]); if (barHitGroup.size() >= _secondary_n_low) { barHitGroupVector.push_back(barHitGroup); } barHitGroup.resize(0); xHitTime = 0; } else { barHitGroup.push_back(emr_bar_hits_tmp[iHit-1]); } } size_t nSeconPartEvents = barHitGroupVector.size(); // Resize the event array to accomodate one extra event per secondary track (n') emr_candidates_tmp = get_emr_events_tmp(nSeconPartEvents); // Associate each barHitGroup to a secondary trigger (daughter candidate) for (size_t iGroup = 0; iGroup < nSeconPartEvents; iGroup++) { EMRBarHitArray barHitArray = barHitGroupVector[iGroup]; for (size_t iHit = 0; iHit < barHitArray.size(); iHit++) { int xPlane = barHitArray[iHit].GetChannel()/_number_of_bars; EMRBarHitTmp barHitTmp; barHitTmp._barhit = barHitArray[iHit]; emr_candidates_tmp[iGroup][xPlane]._barhits.push_back(barHitTmp); } } } void MapCppEMRMCDigitization::fill(MAUS::Spill *spill, EMREventPlaneHitVector& emr_events_tmp, size_t nPartEvents) const { // Get the EMR recon events and the spill data ReconEventPArray *recEvts = spill->GetReconEvents(); EMRSpillData* emrSpill = spill->GetEMRSpillData(); for (size_t iPe = 0; iPe < emr_events_tmp.size(); iPe++) { EMREventTrack* evtTrack = new EMREventTrack; EMRPlaneHitArray plArray; for (size_t iPlane = 0; iPlane < _number_of_planes; iPlane++) { EMRPlaneHit *plHit = new EMRPlaneHit; EMRBarHitVector barHitVector = emr_events_tmp[iPe][iPlane]._barhits; plHit->SetPlane(iPlane); for (size_t iHit = 0; iHit < barHitVector.size(); iHit++) plHit->AddEMRBarHit(barHitVector[iHit]._barhit); plHit->SetOrientation(emr_events_tmp[iPe][iPlane]._orientation); plHit->SetDeltaT(emr_events_tmp[iPe][iPlane]._time); plHit->SetCharge(emr_events_tmp[iPe][iPlane]._charge); plHit->SetPedestalArea(emr_events_tmp[iPe][iPlane]._pedestal_area); plHit->SetSampleArray(emr_events_tmp[iPe][iPlane]._samples); if ( barHitVector.size() || plHit->GetCharge() > 0 ) { plArray.push_back(plHit); for ( size_t iHit = 0; iHit < barHitVector.size(); iHit++ ) emrSpill->AddEMRBarHit(barHitVector[iHit]._barhit); } else { delete plHit; } } evtTrack->SetEMRPlaneHitArray(plArray); if ( iPe < nPartEvents ) { evtTrack->SetType("mother"); evtTrack->SetTrackId(0); EMREvent* evt = new EMREvent; if ( plArray.size() ) { evt->AddEMREventTrack(evtTrack); } else { delete evtTrack; } size_t nRecEvents = spill->GetReconEventSize(); if ( nRecEvents > iPe ) { recEvts->at(iPe)->SetEMREvent(evt); } else { ReconEvent *recEvt = new ReconEvent; recEvt->SetPartEventNumber(iPe); recEvt->SetEMREvent(evt); recEvts->push_back(recEvt); } } else { evtTrack->SetType("candidate"); evtTrack->SetTrackId(iPe-nPartEvents); emrSpill->AddEMREventTrack(evtTrack); } } spill->SetReconEvents(recEvts); spill->SetEMRSpillData(emrSpill); } EMREventPlaneHitVector MapCppEMRMCDigitization::get_emr_events_tmp(size_t nPartEvts) const { EMREventPlaneHitVector emr_events_tmp; emr_events_tmp.resize(nPartEvts); for (size_t iPe = 0; iPe < nPartEvts; iPe++) { emr_events_tmp[iPe].resize(_number_of_planes); // number of planes for (size_t iPlane = 0; iPlane < _number_of_planes; iPlane++) { EMRPlaneHitTmp data; EMRBarHitVector barhits; std::vector xSamples; data._barhits = barhits; data._orientation = iPlane%2; data._time = -1; data._deltat = -1; data._charge = -1; data._pedestal_area = -1; data._samples = xSamples; emr_events_tmp[iPe][iPlane] = data; } } return emr_events_tmp; } } // namespace MAUS