30 #include "Utilities/func/MathUtil.h" 75 double FiberTransmission(
double x);
76 double ScintTime(
double dEdx);
83 void CreatePhotoElectrons(
double meanPE,
double stepTime,
double dist,
double dEdx,
bool isScint,
RateInfo&
info,
double& lambda);
84 void ClusterPhotoElectrons(
int planeId,
int cellId,
int trackId);
180 uint16_t fPosCorrMapIndex = -1;
184 void SetPosCorrMapIndex(
const std::vector<TH2D *> *
const maps,
185 const unsigned int i);
188 void LoadPosCorrHists();
198 std::unique_ptr<std::vector<sim::PhotonSignal> >
fSignals;
212 : fFilterModuleLabel (pset.
get<
std::
string>(
"FilterModuleLabel"))
213 , fPSetNDOS (pset.
get<
fhicl::ParameterSet>(
"ndos"))
214 , fPSetND (pset.
get<
fhicl::ParameterSet>(
"nd"))
215 , fPSetFD (pset.
get<
fhicl::ParameterSet>(
"fd"))
216 , fPSetTB (pset.
get<
fhicl::ParameterSet>(
"tb"))
217 , fSignals (
new std::vector<
sim::PhotonSignal>)
235 produces< std::vector<sim::PhotonSignal> >();
237 produces<BrightnessLevel, art::InRun> ();
254 fScintTime = tfs->
make<TH1D>(
"hScintTime",
";time;photons", 1000, 0, 5000);
259 switch(geom->
DetId()){
273 assert(0 &&
"Unknown detector");
284 fAttenPars = pset.get< std::vector<double> >(
"AttenPars");
287 fStep = pset.get<
double >(
"Step");
296 fScintBetaProb = pset.get< std::vector<double> >(
"ScintBetaProb");
297 fScintBetaTau = pset.get< std::vector<double> >(
"ScintBetaTau");
300 fScintAlphaTau = pset.get< std::vector<double> >(
"ScintAlphaTau");
313 msglevel = pset.get<
int >(
"MessageLevel");
318 const double c = TMath::Ccgs()*1
e-9;
321 sp.find_file(pset.get<
std::string >(
"CollectionRateFile"), fCollectionRateFile);
324 if ( fCollectionRateFile.empty() ||
stat(fCollectionRateFile.c_str(), &sb)!=0 ){
326 throw cet::exception(
"NoCollectionRateFile") <<
"collection rate file " << fCollectionRateFile <<
" not found!\n" << __FILE__ <<
":" << __LINE__ <<
"\n";
329 sp.find_file(pset.get<
std::string >(
"SmearingHistoFile"), fSmearingHistoFile);
330 if ( fSmearingHistoFile.empty() ||
stat(fSmearingHistoFile.c_str(), &sb)!=0 ){
332 throw cet::exception(
"NoSmearingHistoFile") <<
"smearing histo file " << fSmearingHistoFile <<
" not found!\n" << __FILE__ <<
":" << __LINE__ <<
"\n";
335 if (fBrightnessMode !=
"None") {
336 fBrightnessFile = pset.get<
std::string >(
"BrightnessFile");
337 sp.find_file(pset.get<
std::string >(
"BrightnessFile"), fBrightnessFile);
338 if ( fBrightnessFile.empty() ||
stat(fBrightnessFile.c_str(), &sb)!=0 ){
340 throw cet::exception(
"NoBrightnessFile") <<
"fiber brightness histo file " << fBrightnessFile <<
" not found!\n" << __FILE__ <<
":" << __LINE__ <<
"\n";
345 sp.find_file(pset.get<
std::string >(
"TweakFile"), fTweakFile);
346 if ( fTweakFile.empty() ||
stat(fTweakFile.c_str(), &sb)!=0 ) {
348 throw cet::exception(
"NoTweakFile") <<
"tweak file " << fTweakFile <<
" not found!\n" << __FILE__ <<
":" << __LINE__ <<
"\n";
352 sp.find_file(pset.get<
std::string >(
"PosCorrFile"), fPosCorrFile);
353 if(fPosCorrFile.empty()){
355 std::cout <<
"Fiber position efficiency correction disabled\n";
358 if(
stat(fPosCorrFile.c_str(), &sb) != 0)
360 << fPosCorrFile <<
" not found!\n" << __FILE__ <<
":" << __LINE__
367 <<
"can't load histos\n" 368 << __FILE__ <<
":" << __LINE__ <<
"\n";
373 if (fBrightnessMode ==
"ByCell"){
374 std::unique_ptr<BrightnessLevel> pfb(
new BrightnessLevel(fBrightnessFile, fBrightnessMapName, fBrightnessValueName,fBrightnessValueMapName));
375 run.
put(std::move(pfb));
377 else if (fBrightnessMode ==
"ByBin"){
378 std::unique_ptr<BrightnessLevel> pfb(
new BrightnessLevel(fBrightnessFile, fBrightnessMapName, fBrightnessValueName));
379 run.
put(std::move(pfb));
402 +((uint64_t)evt.
subRun() << 32)
403 +((uint64_t)evt.
event() );
405 for (
unsigned int ilist = 0; ilist < flslistHandle->size(); ++ilist) {
406 const std::vector<sim::FLSHit>
hits = flslistHandle->at(ilist).fHits;
407 for (
unsigned int ihit = 0; ihit < hits.size(); ++ihit) {
418 std::unique_ptr<std::vector<sim::PhotonSignal>> pscol(
new std::vector<sim::PhotonSignal>);
422 evt.
put(std::move( pscol ));
430 const double Xa,
const double Ya,
431 const double Xb,
const double Yb)
433 double sum = 0, dem = 0;
437 const double stepsize = 0.01;
439 const double dy = Yb - Ya;
440 const double dx = Xb - Xa;
441 const double seglen =
hypot(dx, dy);
442 const int nstep =
std::max(1,
int(seglen/stepsize));
443 const double thisstepsize = seglen/nstep;
444 const double thisstepstart = thisstepsize/2;
447 const double s = thisstepstart +
step*thisstepsize;
448 const double sfrac = seglen == 0? 0: s/seglen;
455 const double x = -(Xa + sfrac*
dx);
457 const double y = Ya + sfrac*
dy;
459 double eff = themap->GetBinContent(themap->FindBin(x, y));
467 themap->GetBinXYZ(themap->FindBin(x, y),
binx,
biny,
binz);
468 eff = themap->GetBinContent(binx+1, biny);
470 if(eff <= 0) eff = themap->GetBinContent(binx-1,
biny);
471 if(eff <= 0) eff = themap->GetBinContent(binx,
biny+1);
472 if(eff <= 0) eff = themap->GetBinContent(binx,
biny-1);
474 if(eff <= 0)
std::cerr <<
"GetPosCorr should not happen: " << eff
475 <<
" efficiency with " << themap->GetName()
476 <<
" point " << -x <<
" " << y
477 <<
" from path " << Xa <<
" " << Ya <<
" " << Xb <<
" " << Yb
478 <<
" length " << seglen <<
"\n";
482 sum += thisstepsize == 0? eff: thisstepsize*
eff;
488 return seglen == 0? sum: dem == 0? 0: sum/dem;
492 std::vector<TH2D *> & hold,
const bool verbose)
494 for(
int n = 1; ;
n++){
495 TH2D *
h =
dynamic_cast<TH2D *
>(f->Get(Form(
"%s_%d", name,
n)));
497 h->SetDirectory(NULL);
501 if(verbose)
std::cout <<
"ImprovedTransport: using " 502 << hold.size() <<
" " << humanname <<
"\n";
506 <<
"ImprovedTransport: No " << humanname <<
"\n";
512 if(f == NULL || f->IsZombie())
514 <<
"ImprovedTransport: Unable to open " 515 <<
fPosCorrFile <<
"\n" << __FILE__ <<
":" << __LINE__ <<
"\n";
520 "big horizontal cell maps",
523 "big horizontal cell maps with clumped fibers",
526 "small horizontal cell maps",
529 "small horizontal cell maps with clumped fibers",
532 "small horizontal cell maps with bubbles",
535 "small horizontal cell maps with bubbles and clumped fibers",
538 "big vertical cell maps",
541 "big vertical cell maps with fibers in corners",
544 "small vertical cell maps",
547 "small vertical cell maps with fibers in corners",
557 const std::vector<TH2D *> *
const maps,
const unsigned int i)
596 const bool horizontal = planeId%2 == 0;
598 const bool largecell = cellId%16 != 0 && cellId%16 != 15;
600 const bool topcell = horizontal && cellId%32 == 31;
608 const bool hasbubble = topcell &&
627 const bool clump = horizontal &&
fPosCorrHasher.Rndm() < clumpfrac;
632 int dum1 = 0, dum2 = 0;
633 const double frombottom =
642 const double chardist = 100.0 * (1 + 0.012*(!largecell));
649 const bool saggedcorner = sagged?
fPosCorrHasher.Rndm() < 0.5:
false;
651 std::vector<TH2D *> * maps =
674 TH2D * themap = (*maps)[mapi];
681 const double scale = 1/0.134;
685 printf(
"GetPosCorr gave %f for path %f %f %f %f len %f when %s %s %s %s\n",
686 eff, Xa, Ya, Xb, Yb,
hypot(Yb-Ya, Xb-Xa),
687 horizontal?
"horizontal":
"vertical",
688 sagged?
"sagged":
"notsagged",
689 hasbubble?
"bubble":
"nobubble",
690 largecell?
"large":
"small");
701 const double MeVPerGeV = 1000.;
723 double seglength =
util::pythag(Xb - Xa, Yb - Ya, Zb - Za);
726 double DZ =
fGeo->
Plane((
unsigned int)(planeId))->
Cell ((
unsigned int)(cellId ))->
HalfL();
761 double dEdx = (seglength > 0) ? MeVPerGeV*hit.
GetEdep()/seglength : 0.0;
768 for (
int istep = 0; istep <
nsteps; ++istep) {
770 double stepZ = Za + ((Zb-Za)/nsteps)*(istep+0.5);
771 double stepT = Ta + ((Tb-Ta)/nsteps)*(istep+0.5);
775 double tweakCorr(1.0);
780 else tweakCorr =
fTweak_Y.Interpolate(DZ - stepZ);
785 else tweakCorr =
fTweak_YMC.Interpolate(DZ - stepZ);
789 else tweakCorr =
fTweak_Y.Interpolate(DZ - stepZ);
795 std::list<RateInfo>::iterator rate_it;
797 double collectionRate = rate_it->rate;
798 double currZ = stepZ + rate_it->z;
799 double fraction(1.0);
802 if (
std::fabs(currZ) - 0.5*rate_it->zWidth > DZ )
continue;
803 else if (
std::fabs(currZ) + 0.5*rate_it->zWidth < DZ ) fraction = 1.0;
804 else fraction = (DZ - (
std::fabs(currZ) - 0.5*rate_it->zWidth ))/rate_it->zWidth;
805 collectionRate *= fraction;
809 double muGreenScint =
fQuantumEff*collectionRate*brightnessCorr*tweakCorr*
810 positionCorr*factor*scintillationPhotons;
811 double muGreenCkv =
fQuantumEff*collectionRate*brightnessCorr*tweakCorr*
812 positionCorr*factor*cerenkovPhotons;
814 double dist0 = DZ - currZ + pigtail;
815 double dist1 = 3*DZ + currZ + pigtail;
835 for(
unsigned int i = start;
i <
fSignals->size(); ++
i){
836 (*fSignals)[
i].SetPoissonLambda(lambda/N);
845 photon.SetTrackId(trackId);
846 photon.SetPoissonLambda(lambda);
861 double rand0(0), rand1(0), rand2(0);
862 for (
int iPE=0; iPE < nPE; ++iPE) {
863 double time = stepTime;
874 if (isScint) scintTime =
ScintTime(dEdx);
879 double fiberspread = 0;
881 double weight(-1.0), uweight(0.0);
882 while (
weight < uweight) {
885 fiberspread = dist*rand0;
920 std::list< double >::iterator
it;
922 double currTime = *
it;
929 meanTime = meanTime*nPE + currTime;
959 <<
"ImprovedTransport: Unable to open " 961 << __FILE__ <<
":" << __LINE__ <<
"\n";
971 TH2D *htemp = (TH2D*)
f1.Get(
"dT_dZ_CollectionRate")->Clone();
982 TH1D *htemp2 = (TH1D*)
f2.Get(
"Dt_per_z")->Clone();
993 std::vector< double > plane_brightness;
995 plane_brightness.push_back(1.0);
1002 TH1D* hBrightness1D(0);
1003 TH2D* hBrightness2D(0);
1005 hBrightness2D = (TH2D*)
f3.Get(
"BrightnessByCell");
1006 if (!hBrightness2D) {
1012 hBrightness2D = (TH2D*)
f3.Get(
"BrightnessByModule");
1013 if (!hBrightness2D) {
1019 hBrightness1D = (TH1D*)
f3.Get(
"BrightnessByPlane");
1020 if (!hBrightness1D) {
1031 std::cerr <<
"ImprovedTransport: " <<
fBrightnessMode <<
" is not a recognized brightness correction mode. Valid choices are None, ByCell, ByModule, or ByPlane " <<
std::endl;
1055 tmpTweak = (TH1D*)
f4.Get(
"Tweak_X")->Clone();
1064 tmpTweak = (TH1D*)
f4.Get(
"Tweak_Y")->Clone();
1076 std::stringstream hName;
1082 info_pos.
rate = 0.0;
1085 hName.str(
"timePDF_Pos_");
1087 info_pos.
timePDF.SetName(hName.str().c_str());
1088 info_pos.
timePDF.SetBins(nbinsY, ylo, yhi);
1089 info_pos.
timePDF.SetDirectory(0);
1092 info_neg.
rate = 0.0;
1095 hName.str(
"timePDF_Neg_");
1097 info_neg.
timePDF.SetName(hName.str().c_str());
1098 info_neg.
timePDF.SetBins(nbinsY, ylo, yhi);
1099 info_neg.
timePDF.SetDirectory(0);
1103 info_pos.
rate += collectionRate;
1104 info_pos.
timePDF.SetBinContent(iT, collectionRate);
1106 info_neg.
rate += collectionRate;
1107 info_neg.
timePDF.SetBinContent(iT, collectionRate);
1130 mf::LogWarning(
"AttenPars") <<
"No AttenPars entries. Fiber transmission set to 1!";
1134 mf::LogWarning(
"AttenPars") <<
"Odd number of AttenPars entries. Fiber transmission set to 1!";
1159 mf::LogWarning(
"ScintTime") <<
"No ScintBetaProb entries. No scintillation delays applied!";
1162 else if (!nProbAlpha) {
1163 mf::LogWarning(
"ScintTime") <<
"No ScintAlphaProb entries. No scintillation delays applied!";
1166 else if (nProbBeta != nTauBeta) {
1167 mf::LogWarning(
"ScintTime") <<
"Number of entries in ScintBetaProb and ScintBetaTau not the same. No scintillation delays applied!";
1170 else if (nProbAlpha != nTauAlpha) {
1171 mf::LogWarning(
"ScintTime") <<
"Number of entries in ScintBetaAlpha and ScintBetaAlpha not the same. No scintillation delays applied!";
1175 mf::LogWarning(
"ScintTime") <<
"Beta dEdx must be strictly smaller than alpha dEdx. No scintillation delays applied!";
1182 bool useBeta = (u < probBeta);
1185 for (
int iC = 0; iC < nProbBeta; ++iC) {
1193 for (
int iC = 0; iC < nProbAlpha; ++iC) {
double fCerenkovEff
Absorption and re-emission efficiency for Cerenkov photons.
geo::CellUniqueId GetCellUniqueId() const
Unique Cell ID.
bool fApplyFiberSmearing
Turn on fiber smearing.
T max(const caf::Proxy< T > &a, T b)
std::string fBrightnessValueMapName
const XML_Char XML_Encoding * info
float GetEntryX() const
Entry point of the particle (position, time and energy)
ImprovedTransport(fhicl::ParameterSet const &pset)
bool fApplyTweak
Turn on the data driven attenuation correction for each view.
fhicl::ParameterSet fPSetTB
void ClusterPhotoElectrons(int planeId, int cellId, int trackId)
SubRunNumber_t subRun() const
TH2D fCollectionRateHisto
Collection rate for given dT vs. dZ.
double fStep
step size for walking along the track (cm)
std::string fCollectionRateFile
Relative path to the collection rate histo file.
std::vector< TH2D * > fSmallHorizBubbleClumpMap
std::vector< TH2D * > fSmallHorizMap
int GetPlaneID() const
Plane ID.
fvar< T > fabs(const fvar< T > &x)
std::string fBrightnessValueName
CLHEP::RandPoissonQ * fPoisson
std::enable_if_t< std::is_arithmetic< T >::value, T > hypot(T x, T y)
double fYViewFactor
scale factor for the y-view
TH1D fTweak_XMC
A data/mc correction for the X-View in the muon catcher.
TH1D fTweak_Y
A data/mc correction for the Y-View.
std::vector< TH2D * > fSmallVertMap
double fXViewFactor
scale factor for the x-view
void produce(art::Event &evt)
double GetPosCorr(const sim::FLSHit &hit)
int GetCellID() const
Cell ID.
double fTimeClustering
groups photons closer than this many ns apart
const CellGeo * Cell(int icell) const
double FiberTransmission(double x)
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
double fLightLevelScale
shift overall light yield by multiplicative factor
void LoadPosCorrHists()
Reads in the efficiency maps from fPosCorrFile.
Vertical planes which measure X.
unsigned int Ncells() const
Number of cells in this plane.
a struct to hold collection rate info
bool fWriteEmptySignals
Write an empty PhotonSignal if there's 0 PE.
double fFiberSpeedOfLight
Speed of light in the fiber.
base_engine_t & createEngine(seed_t seed)
float GetNCerenkov() const
Get total N Cerenkov photons.
::xsd::cxx::tree::exception< char > exception
art::ServiceHandle< geo::Geometry > fGeo
Handle to the geometry service.
A single unit of energy deposition in the liquid scintillator.
const PlaneGeo * Plane(unsigned int i) const
This slightly less simple photon transport uses a template for photon collection in time and in dista...
art::ProductID put(std::unique_ptr< PROD > &&)
BrightnessLevel * fBrightnessLevel
DEFINE_ART_MODULE(TestTMapFile)
float GetExitX() const
Exit point of the particle (position, time and energy)
::xsd::cxx::tree::time< char, simple_type > time
double fFiberIndexOfRefraction
n for fiber core
a class for transporting photons in a roughly realistic way
TH1D fTweak_YMC
A data/mc correction for the Y-View in the muon catcher.
double fScintBetaDEDX
Approximate dE/dx of betas used to determine PDF.
double fCherenkovTau
Time constant for directly excited PC to PPO scintillation.
void SetPlaneCell(int plane, int cell)
std::unique_ptr< std::vector< sim::PhotonSignal > > fSignals
Produced vector of PhotonSignals.
const CellUniqueId & Id() const
std::list< RateInfo > fRateInfo
List of RateInfo object.
std::vector< double > fScintAlphaTau
Time constants of each exponential component for alphas.
ProductID put(std::unique_ptr< PROD > &&product)
fhicl::ParameterSet fPSetFD
base_engine_t & getEngine() const
View_t View() const
Which coordinate does this plane measure.
std::vector< double > fAttenPars
from Leon's fit
double getBrightnessLevel(unsigned int plane, unsigned int cell) const
Far Detector at Ash River, MN.
std::vector< double > fScintBetaProb
Cumulative probability for each exponential component for betas.
double fDriftGradient
gradient for light level detector aging systematic [fraction/year]
virtual ~ImprovedTransport()
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
std::string fBrightnessMode
correct the fiber brightness by cell, by module, or not at all
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
fvar< T > exp(const fvar< T > &x)
printf("%d Experimental points found\n", nlines)
void SetPosCorrMapIndex(const std::vector< TH2D * > *const maps, const unsigned int i)
Near Detector in the NuMI cavern.
CLHEP::RandExponential * fExp
int Duration()
in units of seconds
std::string fTweakFile
Relative path to the data driven attenuation correction for each view.
bool fSimulateDrift
whether to vary light level as a function of time
EventNumber_t event() const
std::vector< TH2D * > fBigHorizClumpMap
std::vector< TH2D * > fBigHorizMap
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
fhicl::ParameterSet fPSetND
std::list< double > fPETimes
List of photo-electron times (each representing one observed PE)
float GetEdepBirks() const
Get total Energy with Birks suppression deposited into the cell for the whole FLSHit.
uint16_t fPosCorrMapIndex
double fEmissionTau
exponential time constant for green photon emission
double fScintAlphaDEDX
Approximate dE/dx of alphas used to determine PDF.
double ScintTime(double dEdx)
double GetPigtail(const CellUniqueId &id) const
Return length of fiber in cm from end of cell to apd.
bool fUseScintTime
Turn on sampling from PDF of scintillator decay time, see arXiv:1704.02291 for more details...
std::vector< std::vector< double > > fBrightnessMap
The brightness correction by cell.
#define expect(expr, value)
bool fApplyBirks
Use EdepBirks from FLSHits instead of calculating here.
void CreatePhotoElectrons(double meanPE, double stepTime, double dist, double dEdx, bool isScint, RateInfo &info, double &lambda)
unsigned int GetRandomNumberSeed()
std::vector< TH2D * > fBigVertCornerMap
unsigned long long int CellUniqueId
std::vector< TH2D * > fSmallVertCornerMap
T * make(ARGS...args) const
TH1D fFiberSmearingHisto
Smearing histo for photons bouncing around in fiber.
int GetTrackID() const
Track ID.
std::vector< double > fScintAlphaProb
Cumulative probability for each exponential component for alphas.
static double line_efficiency(TH2D *themap, const double Xa, const double Ya, const double Xb, const double Yb)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
double pythag(double x, double y)
2D Euclidean distance
std::vector< double > fScintBetaTau
Time constants of each exponential component for betas.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::string fSmearingHistoFile
Relative path to the fiber smearing histo file.
fhicl::ParameterSet fPSetNDOS
assert(nhit_max >=nhit_nbins)
double fPhotonsPerMeV
blue photon production scale
void GetRandom(TH1 *histo, double &x, double &y, double &z)
std::vector< TH2D * > fSmallHorizClumpMap
std::string fBrightnessFile
Relative path to the brightness correction file.
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
unsigned int NPlanes() const
CLHEP::RandGaussQ * fGauss
double fQuantumEff
APD quantum efficiency.
void load_hset(TFile *f, const char *name, const char *humanname, std::vector< TH2D * > &hold, const bool verbose)
const CellGeo * IdToCell(const CellUniqueId &id, int *iplane, int *icell) const
void SetTimeMean(double t)
std::string fFilterModuleLabel
module that filtered the fls hits
int fDriftReference
drift reference time (ie. where the slope is fixed to 1)
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string fBrightnessMapName
fvar< T > ceil(const fvar< T > &x)
void beginRun(art::Run &run)
TH1D fTweak_X
A data/mc correction for the X-View.
void StepAlongHit(sim::FLSHit &hit)
std::vector< TH2D * > fSmallHorizBubbleMap
std::vector< TH2D * > fBigVertMap
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.