32 #include "Utilities/func/MathUtil.h" 103 produces< std::vector<sim::PhotonSignal> >();
116 << __FILE__ <<
":" << __LINE__ <<
"\n";
140 <<
"can't load smearing histo\n" 141 << __FILE__ <<
":" << __LINE__ <<
"\n";
162 for(
unsigned int i = 0;
i < flslistHandle->size(); ++
i){
168 std::unique_ptr<std::vector<sim::PhotonSignal> > photlist(
new std::vector<sim::PhotonSignal>);
170 const double MeVPerGeV = 1000.;
171 const double c = TMath::Ccgs()*1
e-9;
180 while ( iter_hitlist != hitlist.
end() ) {
183 std::vector<sim::FLSHit>
hits = (*iter_hitlist)->fHits;
185 std::vector<sim::FLSHit>::iterator iter_hit = hits.begin();
186 while ( iter_hit != hits.end() ) {
191 double Xa = (*iter_hit).GetEntryX();
192 double Ya = (*iter_hit).GetEntryY();
193 double Za = (*iter_hit).GetEntryZ();
194 double Xb = (*iter_hit).GetExitX();
195 double Yb = (*iter_hit).GetExitY();
196 double Zb = (*iter_hit).GetExitZ();
199 double Ta = (*iter_hit).GetEntryT();
200 double Tb = (*iter_hit).GetExitT();
202 int planeId = (*iter_hit).GetPlaneID();
203 int cellId = (*iter_hit).GetCellID();
204 int trackId = (*iter_hit).GetTrackID();
207 double seglength =
util::pythag(Xb - Xa, Yb - Ya, Zb - Za);
210 double v = ((Tb-Ta>0) ? seglength/(Tb-Ta) :
c);
215 double DZ = fGeo->
Plane((
unsigned int)(planeId))->
Cell ((
unsigned int)(cellId ))->
HalfL();
218 std::list<double> currtimes;
221 for (
double currS=0; currS<=seglength; currS+=
fStep) {
225 if ( dS > seglength-currS ) dS = seglength-currS;
228 double u = (seglength) ? (currS + dS / 2) / seglength : 0;
231 double currZ = Za + (Zb - Za) * u;
237 double dist0 = DZ - currZ + pigtail;
238 double dist1 = 3*DZ + currZ + pigtail;
243 double muGreenT = green_scale * (seglength ? dS / seglength : 1);
244 if (
fApplyBirks) muGreenT *= iter_hit->GetEdepBirks();
245 else muGreenT *= iter_hit->GetEdep();
250 int nGreen0 = poisson.
fire(muGreen0);
251 int nGreen1 = poisson.
fire(muGreen1);
253 muTot += muGreen0+muGreen1;
267 for (Int_t iP=0; iP<nGreen0+nGreen1; iP++) {
274 double dist = ( iP<nGreen0 ? dist0 : dist1 );
275 double fiberspread = 0;
280 while (weight<uweight) {
282 double rand0(0.0), rand1(0.0), rand2(0.0);
284 fiberspread = dist*rand0;
288 double Delta_dist = fiberspread * cn ;
290 uweight = flat.
fire(0,1);
301 currtimes.push_back(time);
308 const unsigned int nTimes = currtimes.size();
312 std::list<double>::iterator
it = currtimes.begin();
324 while ( it != currtimes.end() ) {
328 sum += (double)(*it);
336 photlist->push_back(photon);
348 photlist->push_back(photon);
359 mf::LogDebug(
"PhotonTransport") <<
"Read " << nHit <<
" FLS hits";
360 mf::LogDebug(
"PhotonTransport") <<
"Created " << photlist->size() <<
" photon signals";
362 evt.
put(std::move(photlist));
382 std::cerr <<
"In PhotonTransport: No AttenPars entries. Fiber transmission set to 1!" <<
std::endl;
386 std::cerr <<
"In PhotonTransport: Odd number of AttenPars entries. Fiber transmission set to 1!" <<
std::endl;
391 for (Int_t iP=0;iP<nPar;iP+=2) {
407 <<
"PhotonTransport: Unable to open " 409 <<
"\nCannot apply fiber smearing as requested!\n" 410 << __FILE__ <<
":" << __LINE__ <<
"\n";
415 TH1D *htemp = (TH1D*)
f.Get(
"Dt_per_z")->Clone();
417 std::cerr <<
"PhotonTransport: Histogram Dt_per_z not found in " void SetPoissonLambda(double l)
a class for transporting photons in a roughly realistic way
bool fApplyFiberSmearing
apply fiber smearing
const CellGeo * Cell(int icell) const
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
double fTimeClustering
groups photons closer than this many ns apart
base_engine_t & createEngine(seed_t seed)
::xsd::cxx::tree::exception< char > exception
const PlaneGeo * Plane(unsigned int i) const
This slightly less simple photon transport uses a template for photon collection in time and in dista...
DEFINE_ART_MODULE(TestTMapFile)
void SetPlaneCell(int plane, int cell)
const CellUniqueId & Id() const
ProductID put(std::unique_ptr< PROD > &&product)
double fEmissionTau
exponential time constant for green photon emission
base_engine_t & getEngine() const
void push_back(Ptr< U > const &p)
T get(std::string const &key) const
double fFiberIndexOfRefraction
n for fiber core
fvar< T > exp(const fvar< T > &x)
double FiberTransmission(double x)
static constexpr Double_t gauss
virtual ~PhotonTransport()
std::vector< double > fAttenPars
from Leon's fit
data_t::const_iterator const_iterator
double GetPigtail(const CellUniqueId &id) const
Return length of fiber in cm from end of cell to apd.
unsigned int GetRandomNumberSeed()
std::string fSmearingHistoFile
relative path to the smearing histo file
void produce(art::Event &evt)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
double pythag(double x, double y)
2D Euclidean distance
double fPhotonsPerMeV
The length (in cm) of the fibre from the end of the cell to the APD.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
double fTimeSpread
spread in time due to everything
void GetRandom(TH1 *histo, double &x, double &y, double &z)
PhotonTransport(fhicl::ParameterSet const &pset)
std::string fFilterModuleLabel
module that filtered the fls hits
bool fApplyBirks
Use EdepBirks from FLSHits instead of calculating here.
bool fNoAttenNoPoissonMode
run special no-attenuation, no-Poisson-fluct mode?
double fQuantumEff
APD quantum efficiency.
double fStep
step size for walking along the track (cm)
void SetTimeMean(double t)
Encapsulate the geometry of one entire detector (near, far, ndos)
double fCollectionEff
fiber collection efficiency