LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
larg4::OpFastScintillation Class Reference

#include "OpFastScintillation.hh"

Inheritance diagram for larg4::OpFastScintillation:

Public Member Functions

 OpFastScintillation (const G4String &processName="Scintillation", G4ProcessType type=fElectromagnetic)
 
 OpFastScintillation (const OpFastScintillation &right)
 
 ~OpFastScintillation ()
 
virtual G4bool IsApplicable (const G4ParticleDefinition &aParticleType)
 
G4double GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
 
G4double GetMeanLifeTime (const G4Track &aTrack, G4ForceCondition *)
 
virtual G4VParticleChange * PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &aTrack, const G4Step &aStep)
 
void SetTrackSecondariesFirst (const G4bool state)
 
void SetFiniteRiseTime (const G4bool state)
 
G4bool GetTrackSecondariesFirst () const
 
G4bool GetFiniteRiseTime () const
 
void SetScintillationYieldFactor (const G4double yieldfactor)
 
G4double GetScintillationYieldFactor () const
 
void SetScintillationExcitationRatio (const G4double excitationratio)
 
G4double GetScintillationExcitationRatio () const
 
G4PhysicsTable * GetFastIntegralTable () const
 
G4PhysicsTable * GetSlowIntegralTable () const
 
void AddSaturation (G4EmSaturation *sat)
 
void RemoveSaturation ()
 
G4EmSaturation * GetSaturation () const
 
void SetScintillationByParticleType (const G4bool)
 
G4bool GetScintillationByParticleType () const
 
void DumpPhysicsTable () const
 
std::vector< double > GetVUVTime (double, int) const
 
std::vector< double > GetVisibleTimeOnlyCathode (double, int) const
 

Protected Member Functions

void BuildThePhysicsTable ()
 
bool RecordPhotonsProduced (const G4Step &aStep, double N)
 

Protected Attributes

G4PhysicsTable * theSlowIntegralTable
 
G4PhysicsTable * theFastIntegralTable
 
G4bool fTrackSecondariesFirst
 
G4bool fFiniteRiseTime
 
G4double YieldFactor
 
G4double ExcitationRatio
 
G4bool scintillationByParticleType
 

Private Member Functions

G4double single_exp (G4double t, G4double tau2) const
 
G4double bi_exp (G4double t, G4double tau1, G4double tau2) const
 
G4double scint_time (const G4Step &aStep, G4double ScintillationTime, G4double ScintillationRiseTime) const
 
std::vector< double > propagation_time (G4ThreeVector x0, int OpChannel, int NPhotons, bool Reflected=false) const
 
G4double sample_time (G4double tau1, G4double tau2) const
 
double reemission_energy () const
 
void average_position (G4Step const &aStep, double *xzyPos) const
 
void ProcessStep (const G4Step &step)
 

Private Attributes

std::map< double, double > tpbemission
 
CLHEP::RandGeneral * rgen0
 
G4EmSaturation * emSaturation
 
TF1 * ParPropTimeTF1
 
float const * ReflT0s
 
TF1 const * functions_vuv [8]
 
TF1 const * functions_vis [5]
 
double fd_break
 
double fd_max
 
double ftf1_sampling_factor
 
double ft0_max
 
double ft0_break_point
 
bool bPropagate
 Whether propagation of photons is enabled. More...
 

Detailed Description

Definition at line 114 of file OpFastScintillation.hh.

Constructor & Destructor Documentation

larg4::OpFastScintillation::OpFastScintillation ( const G4String &  processName = "Scintillation",
G4ProcessType  type = fElectromagnetic 
)

Definition at line 151 of file OpFastScintillation.cxx.

References bPropagate, BuildThePhysicsTable(), emSaturation, ExcitationRatio, fd_break, fd_max, fFiniteRiseTime, ft0_break_point, ft0_max, ftf1_sampling_factor, fTrackSecondariesFirst, functions_vis, functions_vuv, phot::PhotonVisibilityService::IncludePropTime(), rgen0, detinfo::LArProperties::ScintByParticleType(), scintillationByParticleType, phot::PhotonVisibilityService::SetDirectLightPropFunctions(), phot::PhotonVisibilityService::SetReflectedCOLightPropFunctions(), theFastIntegralTable, theSlowIntegralTable, tpbemission, and YieldFactor.

152  : G4VRestDiscreteProcess(processName, type)
153  , bPropagate(!(art::ServiceHandle<sim::LArG4Parameters>()->NoPhotonPropagation()))
154  {
155  G4cout<<"BEA1 timing distribution"<<G4endl;
156  SetProcessSubType(25);
157 
158  fTrackSecondariesFirst = false;
159  fFiniteRiseTime = false;
160 
161 
162  YieldFactor=1.0;
163  ExcitationRatio = 1.0;
164 
165  const detinfo::LArProperties* larp = lar::providerFrom<detinfo::LArPropertiesService>();
166 
167  scintillationByParticleType = larp->ScintByParticleType();
168 
169  theFastIntegralTable = NULL;
170  theSlowIntegralTable = NULL;
171 
172  if (verboseLevel>0) {
173  G4cout << GetProcessName() << " is created " << G4endl;
174  }
175 
176  BuildThePhysicsTable();
177 
178  emSaturation = NULL;
179 
180  if (bPropagate) {
181  art::ServiceHandle<phot::PhotonVisibilityService> pvs;
182  if(pvs->IncludePropTime()) {
183  pvs->SetDirectLightPropFunctions(functions_vuv, fd_break, fd_max, ftf1_sampling_factor);
184  pvs->SetReflectedCOLightPropFunctions(functions_vis, ft0_max, ft0_break_point);
185  }
186  }
187 
188  tpbemission=lar::providerFrom<detinfo::LArPropertiesService>()->TpbEm();
189  const int nbins=tpbemission.size();
190  double * parent=new double[nbins];
191  int ii=0;
192  for( std::map<double, double>::iterator iter = tpbemission.begin(); iter != tpbemission.end(); ++iter)
193  {
194  parent[ii++]=(*iter).second;
195  }
196  rgen0 = new CLHEP::RandGeneral(parent,nbins);
197 
198  }
phot::PhotonVisibilityService::SetDirectLightPropFunctions
void SetDirectLightPropFunctions(TF1 const *functions[8], double &d_break, double &d_max, double &tf1_sampling_factor) const
Definition: PhotonVisibilityService_service.cc:616
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
larg4::OpFastScintillation::tpbemission
std::map< double, double > tpbemission
Definition: OpFastScintillation.hh:279
phot::PhotonVisibilityService::SetReflectedCOLightPropFunctions
void SetReflectedCOLightPropFunctions(TF1 const *functions[5], double &t0_max, double &t0_break_point) const
Definition: PhotonVisibilityService_service.cc:633
larg4::OpFastScintillation::bPropagate
bool bPropagate
Whether propagation of photons is enabled.
Definition: OpFastScintillation.hh:299
detinfo::LArProperties::ScintByParticleType
virtual bool ScintByParticleType() const =0
larg4::OpFastScintillation::OpFastScintillation ( const OpFastScintillation &  right)

Definition at line 202 of file OpFastScintillation.cxx.

References BuildThePhysicsTable(), emSaturation, ExcitationRatio, fFiniteRiseTime, fTrackSecondariesFirst, GetFastIntegralTable(), GetFiniteRiseTime(), GetSaturation(), GetScintillationByParticleType(), GetScintillationExcitationRatio(), GetScintillationYieldFactor(), GetSlowIntegralTable(), GetTrackSecondariesFirst(), scintillationByParticleType, theFastIntegralTable, theSlowIntegralTable, and YieldFactor.

203  : G4VRestDiscreteProcess(rhs.GetProcessName(), rhs.GetProcessType())
204  {
205  theSlowIntegralTable = rhs.GetSlowIntegralTable();
206  theFastIntegralTable = rhs.GetFastIntegralTable();
207 
208  fTrackSecondariesFirst = rhs.GetTrackSecondariesFirst();
209  fFiniteRiseTime = rhs.GetFiniteRiseTime();
210  YieldFactor = rhs.GetScintillationYieldFactor();
211  ExcitationRatio = rhs.GetScintillationExcitationRatio();
212  scintillationByParticleType = rhs.GetScintillationByParticleType();
213  emSaturation = rhs.GetSaturation();
214 
215  BuildThePhysicsTable();
216  }
larg4::OpFastScintillation::~OpFastScintillation ( )

Definition at line 223 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

224  {
225  if (theFastIntegralTable != NULL) {
226  theFastIntegralTable->clearAndDestroy();
227  delete theFastIntegralTable;
228  }
229  if (theSlowIntegralTable != NULL) {
230  theSlowIntegralTable->clearAndDestroy();
231  delete theSlowIntegralTable;
232  }
233 
234  }

Member Function Documentation

void larg4::OpFastScintillation::AddSaturation ( G4EmSaturation *  sat)
inline

Definition at line 209 of file OpFastScintillation.hh.

References emSaturation.

209 { emSaturation = sat; }
G4VParticleChange * larg4::OpFastScintillation::AtRestDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)
virtual

Definition at line 244 of file OpFastScintillation.cxx.

References PostStepDoIt().

249  {
250  return OpFastScintillation::PostStepDoIt(aTrack, aStep);
251  }
larg4::OpFastScintillation::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: OpFastScintillation.cxx:257
void larg4::OpFastScintillation::average_position ( G4Step const &  aStep,
double *  xzyPos 
) const
private

Definition at line 1049 of file OpFastScintillation.cxx.

Referenced by RecordPhotonsProduced().

1050  {
1051  CLHEP::Hep3Vector prePoint = (aStep.GetPreStepPoint())->GetPosition();
1052  CLHEP::Hep3Vector postPoint = (aStep.GetPostStepPoint())->GetPosition();
1053  xyzPos[0] = ( ( (prePoint.getX() + postPoint.getX() ) / 2.0) / CLHEP::cm );
1054  xyzPos[1] = ( ( (prePoint.getY() + postPoint.getY() ) / 2.0) / CLHEP::cm );
1055  xyzPos[2] = ( ( (prePoint.getZ() + postPoint.getZ() ) / 2.0) / CLHEP::cm );
1056  }
G4double larg4::OpFastScintillation::bi_exp ( G4double  t,
G4double  tau1,
G4double  tau2 
) const
inlineprivate

Definition at line 412 of file OpFastScintillation.hh.

Referenced by sample_time().

413 {
414  return std::exp(-1.0*t/tau2)*(1-std::exp(-1.0*t/tau1))/tau2/tau2*(tau1+tau2);
415 }
void larg4::OpFastScintillation::BuildThePhysicsTable ( )
protected

Definition at line 759 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

Referenced by GetScintillationByParticleType(), and OpFastScintillation().

760  {
761  if (theFastIntegralTable && theSlowIntegralTable) return;
762 
763  const G4MaterialTable* theMaterialTable =
764  G4Material::GetMaterialTable();
765  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
766 
767  // create new physics table
768 
769  if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
770  if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
771 
772  // loop for materials
773 
774  for (G4int i=0 ; i < numOfMaterials; i++)
775  {
776  G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
777  new G4PhysicsOrderedFreeVector();
778  G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
779  new G4PhysicsOrderedFreeVector();
780 
781  // Retrieve vector of scintillation wavelength intensity for
782  // the material from the material's optical properties table.
783 
784  G4Material* aMaterial = (*theMaterialTable)[i];
785 
786  G4MaterialPropertiesTable* aMaterialPropertiesTable =
787  aMaterial->GetMaterialPropertiesTable();
788 
789  if (aMaterialPropertiesTable) {
790 
791  G4MaterialPropertyVector* theFastLightVector =
792  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
793 
794  if (theFastLightVector) {
795 
796  // Retrieve the first intensity point in vector
797  // of (photon energy, intensity) pairs
798 
799  G4double currentIN = (*theFastLightVector)[0];
800 
801  if (currentIN >= 0.0) {
802 
803  // Create first (photon energy, Scintillation
804  // Integral pair
805 
806  G4double currentPM = theFastLightVector->Energy(0);
807 
808  G4double currentCII = 0.0;
809 
810  aPhysicsOrderedFreeVector->
811  InsertValues(currentPM , currentCII);
812 
813  // Set previous values to current ones prior to loop
814 
815  G4double prevPM = currentPM;
816  G4double prevCII = currentCII;
817  G4double prevIN = currentIN;
818 
819  // loop over all (photon energy, intensity)
820  // pairs stored for this material
821 
822  for (size_t i = 1;
823  i < theFastLightVector->GetVectorLength();
824  i++)
825  {
826  currentPM = theFastLightVector->Energy(i);
827  currentIN = (*theFastLightVector)[i];
828 
829  currentCII = 0.5 * (prevIN + currentIN);
830 
831  currentCII = prevCII +
832  (currentPM - prevPM) * currentCII;
833 
834  aPhysicsOrderedFreeVector->
835  InsertValues(currentPM, currentCII);
836 
837  prevPM = currentPM;
838  prevCII = currentCII;
839  prevIN = currentIN;
840  }
841 
842  }
843  }
844 
845  G4MaterialPropertyVector* theSlowLightVector =
846  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
847 
848  if (theSlowLightVector) {
849 
850  // Retrieve the first intensity point in vector
851  // of (photon energy, intensity) pairs
852 
853  G4double currentIN = (*theSlowLightVector)[0];
854 
855  if (currentIN >= 0.0) {
856 
857  // Create first (photon energy, Scintillation
858  // Integral pair
859 
860  G4double currentPM = theSlowLightVector->Energy(0);
861 
862  G4double currentCII = 0.0;
863 
864  bPhysicsOrderedFreeVector->
865  InsertValues(currentPM , currentCII);
866 
867  // Set previous values to current ones prior to loop
868 
869  G4double prevPM = currentPM;
870  G4double prevCII = currentCII;
871  G4double prevIN = currentIN;
872 
873  // loop over all (photon energy, intensity)
874  // pairs stored for this material
875 
876  for (size_t i = 1;
877  i < theSlowLightVector->GetVectorLength();
878  i++)
879  {
880  currentPM = theSlowLightVector->Energy(i);
881  currentIN = (*theSlowLightVector)[i];
882 
883  currentCII = 0.5 * (prevIN + currentIN);
884 
885  currentCII = prevCII +
886  (currentPM - prevPM) * currentCII;
887 
888  bPhysicsOrderedFreeVector->
889  InsertValues(currentPM, currentCII);
890 
891  prevPM = currentPM;
892  prevCII = currentCII;
893  prevIN = currentIN;
894  }
895 
896  }
897  }
898  }
899 
900  // The scintillation integral(s) for a given material
901  // will be inserted in the table(s) according to the
902  // position of the material in the material table.
903 
904  theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
905  theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
906 
907  }
908  }
void larg4::OpFastScintillation::DumpPhysicsTable ( ) const
inline

Definition at line 380 of file OpFastScintillation.hh.

References theFastIntegralTable, and theSlowIntegralTable.

Referenced by GetScintillationByParticleType().

381 {
382  if (theFastIntegralTable) {
383  G4int PhysicsTableSize = theFastIntegralTable->entries();
384  G4PhysicsOrderedFreeVector *v;
385 
386  for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
387  {
388  v = (G4PhysicsOrderedFreeVector*)(*theFastIntegralTable)[i];
389  v->DumpValues();
390  }
391  }
392 
393  if (theSlowIntegralTable) {
394  G4int PhysicsTableSize = theSlowIntegralTable->entries();
395  G4PhysicsOrderedFreeVector *v;
396 
397  for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
398  {
399  v = (G4PhysicsOrderedFreeVector*)(*theSlowIntegralTable)[i];
400  v->DumpValues();
401  }
402  }
403 }
G4PhysicsTable * larg4::OpFastScintillation::GetFastIntegralTable ( ) const
inline

Definition at line 374 of file OpFastScintillation.hh.

References theFastIntegralTable.

Referenced by OpFastScintillation().

375 {
376  return theFastIntegralTable;
377 }
G4bool larg4::OpFastScintillation::GetFiniteRiseTime ( ) const
inline

Definition at line 338 of file OpFastScintillation.hh.

References fFiniteRiseTime.

Referenced by OpFastScintillation().

339 {
340  return fFiniteRiseTime;
341 }
G4double larg4::OpFastScintillation::GetMeanFreePath ( const G4Track &  aTrack,
G4double  ,
G4ForceCondition *  condition 
)

Definition at line 927 of file OpFastScintillation.cxx.

930  {
931  *condition = StronglyForced;
932 
933  return DBL_MAX;
934 
935  }
G4double larg4::OpFastScintillation::GetMeanLifeTime ( const G4Track &  aTrack,
G4ForceCondition *  condition 
)

Definition at line 941 of file OpFastScintillation.cxx.

943  {
944  *condition = Forced;
945 
946  return DBL_MAX;
947 
948  }
G4EmSaturation* larg4::OpFastScintillation::GetSaturation ( ) const
inline

Definition at line 215 of file OpFastScintillation.hh.

References emSaturation, and SetScintillationByParticleType().

Referenced by OpFastScintillation().

215 { return emSaturation; }
G4bool larg4::OpFastScintillation::GetScintillationByParticleType ( ) const
inline

Definition at line 222 of file OpFastScintillation.hh.

References BuildThePhysicsTable(), DumpPhysicsTable(), GetVisibleTimeOnlyCathode(), GetVUVTime(), RecordPhotonsProduced(), and scintillationByParticleType.

Referenced by OpFastScintillation().

223  { return scintillationByParticleType; }
G4double larg4::OpFastScintillation::GetScintillationExcitationRatio ( ) const
inline

Definition at line 362 of file OpFastScintillation.hh.

References ExcitationRatio.

Referenced by OpFastScintillation().

363 {
364  return ExcitationRatio;
365 }
G4double larg4::OpFastScintillation::GetScintillationYieldFactor ( ) const
inline

Definition at line 350 of file OpFastScintillation.hh.

References YieldFactor.

Referenced by OpFastScintillation().

351 {
352  return YieldFactor;
353 }
G4PhysicsTable * larg4::OpFastScintillation::GetSlowIntegralTable ( ) const
inline

Definition at line 368 of file OpFastScintillation.hh.

References theSlowIntegralTable.

Referenced by OpFastScintillation().

369 {
370  return theSlowIntegralTable;
371 }
G4bool larg4::OpFastScintillation::GetTrackSecondariesFirst ( ) const
inline

Definition at line 332 of file OpFastScintillation.hh.

References fTrackSecondariesFirst.

Referenced by OpFastScintillation().

333 {
334  return fTrackSecondariesFirst;
335 }
std::vector< double > larg4::OpFastScintillation::GetVisibleTimeOnlyCathode ( double  t0,
int  number_photons 
) const

Definition at line 1136 of file OpFastScintillation.cxx.

References larg4::finter_d(), ft0_break_point, ft0_max, ftf1_sampling_factor, functions_vis, and larg4::LandauPlusExpoFinal().

Referenced by GetScintillationByParticleType(), and propagation_time().

1136  {
1137  //-----Distances in cm and times in ns-----//
1138  //gRandom->SetSeed(0);
1139 
1140  std::vector<double> arrival_time_distrb;
1141  arrival_time_distrb.clear();
1142  // Parametrization data:
1143 
1144  if(t0 < 8 || t0 > ft0_max) {
1145  G4cout<<"WARNING: Parametrization of Cathode-Only reflected Light not fully reliable"<<G4endl;
1146  G4cout<<"Too close/far to the PMT -> set 0 Visible photons(?)!!!!!!"<<G4endl;
1147  return arrival_time_distrb;
1148  }
1149  //signals (remember this is transportation) no longer than 1us
1150  const double signal_t_range = 1000.;
1151  double pars_landau[3] = {functions_vis[1]->Eval(t0), functions_vis[2]->Eval(t0),
1152  pow(10.,functions_vis[0]->Eval(t0))};
1153  double pars_expo[2] = {functions_vis[3]->Eval(t0), functions_vis[4]->Eval(t0)};
1154  if(t0 > ft0_break_point) {
1155  pars_landau[0] = -0.798934 + 1.06216*t0;
1156  pars_landau[1] = functions_vis[2]->Eval(ft0_break_point);
1157  pars_landau[2] = pow(10.,functions_vis[0]->Eval(ft0_break_point));
1158  pars_expo[0] = functions_vis[3]->Eval(ft0_break_point);
1159  pars_expo[1] = functions_vis[4]->Eval(ft0_break_point);
1160  }
1161 
1162  // Defining the two functions (Landau + Exponential) describing the timing vs t0
1163  TF1 *flandau = new TF1("flandau","[2]*TMath::Landau(x,[0],[1])",0,signal_t_range/2);
1164  flandau->SetParameters(pars_landau);
1165  TF1 *fexpo = new TF1("fexpo","expo",0, signal_t_range/2);
1166  fexpo->SetParameters(pars_expo);
1167  //this is to find the intersection point between the two functions:
1168  TF1 *fint = new TF1("fint",finter_d,flandau->GetMaximumX(),2*t0,5);
1169  double parsInt[5] = {pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1170  fint->SetParameters(parsInt);
1171  double t_int = fint->GetMinimumX();
1172  double minVal = fint->Eval(t_int);
1173  //the functions must intersect!!!
1174  if(minVal>0.015)
1175  G4cout<<"WARNING: Parametrization of Direct Light discontinuous (landau + expo)!!!!!!"<<G4endl;
1176 
1177  TF1 *fVisTiming = new TF1("fTiming",LandauPlusExpoFinal,0,signal_t_range,6);
1178  double parsfinal[6] = {t_int, pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1179  fVisTiming->SetParameters(parsfinal);
1180  // Set the number of points used to sample the function
1181 
1182  int f_sampling = 1000;
1183  if(t0 < 20)
1184  f_sampling *= ftf1_sampling_factor;
1185  fVisTiming->SetNpx(f_sampling);
1186 
1187  for(int i=0; i<number_photons; i++)
1188  arrival_time_distrb.push_back(fVisTiming->GetRandom());
1189 
1190  //deleting ...
1191 
1192  delete flandau;
1193  delete fexpo;
1194  delete fint;
1195  delete fVisTiming;
1196  G4cout<<"Timing distribution BEAMAUS!"<<G4endl;
1197  return arrival_time_distrb;
1198  }
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
larg4::finter_d
double finter_d(double *x, double *par)
Definition: OpFastScintillation.cxx:1200
larg4::LandauPlusExpoFinal
double LandauPlusExpoFinal(double *x, double *par)
Definition: OpFastScintillation.cxx:1207
std::vector< double > larg4::OpFastScintillation::GetVUVTime ( double  distance,
int  number_photons 
) const

Definition at line 1066 of file OpFastScintillation.cxx.

References fd_break, fd_max, larg4::finter_d(), ftf1_sampling_factor, functions_vuv, and larg4::LandauPlusExpoFinal().

Referenced by GetScintillationByParticleType(), and propagation_time().

1066  {
1067 
1068  //-----Distances in cm and times in ns-----//
1069  //gRandom->SetSeed(0);
1070  std::vector<double> arrival_time_distrb;
1071  arrival_time_distrb.clear();
1072 
1073  // Parametrization data:
1074 
1075  if(distance < 10 || distance > fd_max) {
1076  G4cout<<"WARNING: Parametrization of Direct Light not fully reliable"<<G4endl;
1077  G4cout<<"Too close/far to the PMT -> set 0 VUV photons(?)!!!!!!"<<G4endl;
1078  return arrival_time_distrb;
1079  }
1080  //signals (remember this is transportation) no longer than 1us
1081  const double signal_t_range = 1000.;
1082  const double vuv_vgroup = 10.13;//cm/ns
1083  double t_direct = distance/vuv_vgroup;
1084  // Defining the two functions (Landau + Exponential) describing the timing vs distance
1085  double pars_landau[3] = {functions_vuv[1]->Eval(distance), functions_vuv[2]->Eval(distance),
1086  pow(10.,functions_vuv[0]->Eval(distance))};
1087  if(distance > fd_break) {
1088  pars_landau[0]=functions_vuv[6]->Eval(distance);
1089  pars_landau[1]=functions_vuv[2]->Eval(fd_break);
1090  pars_landau[2]=pow(10.,functions_vuv[5]->Eval(distance));
1091  }
1092  TF1 *flandau = new TF1("flandau","[2]*TMath::Landau(x,[0],[1])",0,signal_t_range/2);
1093  flandau->SetParameters(pars_landau);
1094  double pars_expo[2] = {functions_vuv[3]->Eval(distance), functions_vuv[4]->Eval(distance)};
1095  if(distance > (fd_break - 50.)) {
1096  pars_expo[0] = functions_vuv[7]->Eval(distance);
1097  pars_expo[1] = functions_vuv[4]->Eval(fd_break - 50.);
1098  }
1099  TF1 *fexpo = new TF1("fexpo","expo",0, signal_t_range/2);
1100  fexpo->SetParameters(pars_expo);
1101  //this is to find the intersection point between the two functions:
1102  TF1 *fint = new TF1("fint",finter_d,flandau->GetMaximumX(),3*t_direct,5);
1103  double parsInt[5] = {pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1104  fint->SetParameters(parsInt);
1105  double t_int = fint->GetMinimumX();
1106  double minVal = fint->Eval(t_int);
1107  //the functions must intersect!!!
1108  if(minVal>0.015)
1109  G4cout<<"WARNING: Parametrization of Direct Light discontinuous (landau + expo)!!!!!!"<<G4endl;
1110 
1111  TF1 *fVUVTiming = new TF1("fTiming",LandauPlusExpoFinal,0,signal_t_range,6);
1112  double parsfinal[6] = {t_int, pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1113  fVUVTiming->SetParameters(parsfinal);
1114  // Set the number of points used to sample the function
1115 
1116  int f_sampling = 1000;
1117  if(distance < 50)
1118  f_sampling *= ftf1_sampling_factor;
1119  fVUVTiming->SetNpx(f_sampling);
1120 
1121  for(int i=0; i<number_photons; i++)
1122  arrival_time_distrb.push_back(fVUVTiming->GetRandom());
1123 
1124  //deleting ...
1125  delete flandau;
1126  delete fexpo;
1127  delete fint;
1128  delete fVUVTiming;
1129  G4cout<<"BEAMAUS timing distribution hecha"<<G4endl;
1130  return arrival_time_distrb;
1131  }
larg4::finter_d
double finter_d(double *x, double *par)
Definition: OpFastScintillation.cxx:1200
larg4::LandauPlusExpoFinal
double LandauPlusExpoFinal(double *x, double *par)
Definition: OpFastScintillation.cxx:1207
G4bool larg4::OpFastScintillation::IsApplicable ( const G4ParticleDefinition &  aParticleType)
inlinevirtual

Definition at line 311 of file OpFastScintillation.hh.

Referenced by larg4::FastOpticalPhysics::ConstructProcess().

312 {
313  if (aParticleType.GetParticleName() == "opticalphoton") return false;
314  if (aParticleType.IsShortLived()) return false;
315 
316  return true;
317 }
G4VParticleChange * larg4::OpFastScintillation::PostStepDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)
virtual

Definition at line 257 of file OpFastScintillation.cxx.

References larg4::IonizationAndScintillation::Instance(), larg4::IonizationAndScintillation::NumberScintillationPhotons(), RecordPhotonsProduced(), and larg4::IonizationAndScintillation::Reset().

Referenced by AtRestDoIt().

263  {
264  aParticleChange.Initialize(aTrack);
265 
266  // Check that we are in a material with a properties table, if not
267  // just return
268  const G4Material* aMaterial = aTrack.GetMaterial();
269  G4MaterialPropertiesTable* aMaterialPropertiesTable =
270  aMaterial->GetMaterialPropertiesTable();
271  if (!aMaterialPropertiesTable)
272  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
273 
274  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
275 
276  G4ThreeVector x0 = pPreStepPoint->GetPosition();
277  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
278 
279 
281  // This is the old G4 way - but we do things differently - Ben J, Oct Nov 2012.
283  //
284  // if (MeanNumberOfPhotons > 10.)
285  // {
286  // G4double sigma = ResolutionScale * std::sqrt(MeanNumberOfPhotons);
287  // NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
288  // }
289  // else
290  // {
291  // NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
292  // }
293  //
294  //
295  //
296  // if (NumPhotons <= 0)
297  // {
298  // // return unchanged particle and no secondaries
299  //
300  // aParticleChange.SetNumberOfSecondaries(0);
301  //
302  // return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
303  // }
304  //
305  //
306  // aParticleChange.SetNumberOfSecondaries(NumPhotons);
307  //
308  //
309  // if (fTrackSecondariesFirst) {
310  // if (aTrack.GetTrackStatus() == fAlive )
311  // aParticleChange.ProposeTrackStatus(fSuspend);
312  // }
313  //
314  //
315  //
316  //
318  //
319 
320 
322  // The fast sim way - Ben J, Nov 2012
324  //
325  //
326 
327  // We don't want to produce any trackable G4 secondaries
328  aParticleChange.SetNumberOfSecondaries(0);
329 
330 
331  // Retrieve the Scintillation Integral for this material
332  // new G4PhysicsOrderedFreeVector allocated to hold CII's
333 
334 
335  // Some explanation for later improvements to scint yield code:
336  //
337  // What does G4 do here?
338  // It produces light in 2 steps, fast (scnt=1) then slow (scnt=2)
339  //
340  // The ratio of slow photons to fast photons is related by the yieldratio
341  // parameter. G4's poisson fluctuating scheme is a bit different to ours
342  // - we should check that they are equivalent.
343  //
344  // G4 poisson fluctuates the number of initial photons then divides them
345  // with a constant factor between fast + slow, whereas we poisson
346  // fluctuate separateyly the fast and slow detection numbers.
347  //
348 
349  // get the number of photons produced from the IonizationAndScintillation
350  // singleton
351  larg4::IonizationAndScintillation::Instance()->Reset(&aStep);
352  double MeanNumberOfPhotons = larg4::IonizationAndScintillation::Instance()->NumberScintillationPhotons();
353 // double stepEnergy = larg4::IonizationAndScintillation::Instance()->VisibleEnergyDeposit()/CLHEP::MeV;
354  RecordPhotonsProduced(aStep, MeanNumberOfPhotons);//, stepEnergy);
355 
356  if (verboseLevel>0) {
357  G4cout << "\n Exiting from OpFastScintillation::DoIt -- NumberOfSecondaries = "
358  << aParticleChange.GetNumberOfSecondaries() << G4endl;
359  }
360 
361 
362  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
363  }
larg4::OpFastScintillation::RecordPhotonsProduced
bool RecordPhotonsProduced(const G4Step &aStep, double N)
Definition: OpFastScintillation.cxx:392
larg4::IonizationAndScintillation::Instance
static IonizationAndScintillation * Instance()
Definition: IonizationAndScintillation.cxx:38
void larg4::OpFastScintillation::ProcessStep ( const G4Step &  step)
private

Definition at line 367 of file OpFastScintillation.cxx.

References larg4::OpDetPhotonTable::AddEnergyDeposit(), larg4::ParticleListAction::GetCurrentTrackID(), and larg4::OpDetPhotonTable::Instance().

Referenced by RecordPhotonsProduced().

368  {
369  if(step.GetTotalEnergyDeposit() <= 0) return;
370 
371  OpDetPhotonTable::Instance()->AddEnergyDeposit
372  (-1,
373  -1,
374  (double)(step.GetTotalEnergyDeposit()/CLHEP::MeV), //energy in MeV
375  (float)(step.GetPreStepPoint()->GetPosition().x()/CLHEP::cm),
376  (float)(step.GetPreStepPoint()->GetPosition().y()/CLHEP::cm),
377  (float)(step.GetPreStepPoint()->GetPosition().z()/CLHEP::cm),
378  (float)(step.GetPostStepPoint()->GetPosition().x()/CLHEP::cm),
379  (float)(step.GetPostStepPoint()->GetPosition().y()/CLHEP::cm),
380  (float)(step.GetPostStepPoint()->GetPosition().z()/CLHEP::cm),
381  (double)(step.GetPreStepPoint()->GetGlobalTime()),
382  (double)(step.GetPostStepPoint()->GetGlobalTime()),
383  //step.GetTrack()->GetTrackID(),
384  ParticleListAction::GetCurrentTrackID(),
385  step.GetTrack()->GetParticleDefinition()->GetPDGEncoding(),
386  step.GetPreStepPoint()->GetPhysicalVolume()->GetName()
387  );
388  }
larg4::OpDetPhotonTable::AddEnergyDeposit
void AddEnergyDeposit(int n_elec, int n_photon, double energy, float start_x, float start_y, float start_z, float end_x, float end_y, float end_z, double start_time, double end_time, int trackid, int pdgcode, std::string const &vol="EMPTY")
Definition: OpDetPhotonTable.cxx:205
larg4::OpDetPhotonTable::Instance
static OpDetPhotonTable * Instance(bool LitePhotons=false)
Definition: OpDetPhotonTable.cxx:36
std::vector< double > larg4::OpFastScintillation::propagation_time ( G4ThreeVector  x0,
int  OpChannel,
int  NPhotons,
bool  Reflected = false 
) const
private

Definition at line 971 of file OpFastScintillation.cxx.

References geo::OpDetGeo::GetCenter(), GetVisibleTimeOnlyCathode(), GetVUVTime(), phot::PhotonVisibilityService::IncludeParPropTime(), phot::PhotonVisibilityService::IncludePropTime(), geo::GeometryCore::OpDetGeoFromOpDet(), ParPropTimeTF1, and ReflT0s.

Referenced by RecordPhotonsProduced().

972  {
973  static art::ServiceHandle<geo::Geometry> geo;
974  static art::ServiceHandle<phot::PhotonVisibilityService> pvs;
975 
976  // Initialize vector of the right length with all 0's
977  std::vector<double> arrival_time_dist(NPhotons, 0);
978 
979 
980  if (pvs->IncludeParPropTime() && pvs->IncludePropTime()) {
981  throw cet::exception("OpFastScintillation") << "Cannot have both propagation time models simultaneously.";
982  }
983 
984  else if (pvs->IncludeParPropTime() && !(ParPropTimeTF1 && (ParPropTimeTF1[OpChannel].GetNdim()==1)) )
985  {
986  //Warning: TF1::GetNdim()==1 will tell us if the TF1 is really defined or it is the default one.
987  //This will fix a segfault when using timing and interpolation.
988  G4cout << "WARNING: Requested parameterized timing, but no function found. Not applying propagation time." << G4endl;
989  }
990 
991  else if (pvs->IncludeParPropTime()) {
992  if (Reflected)
993  throw cet::exception("OpFastScintillation") << "No parameterized propagation time for reflected light";
994 
995  for (int i = 0; i < NPhotons; i++) {
996  arrival_time_dist[i] = ParPropTimeTF1[OpChannel].GetRandom();
997  }
998  }
999 
1000  else if (pvs->IncludePropTime()) {
1001  // Get VUV photons arrival time distribution from the parametrization
1002  double OpDetCenter[3];
1003  geo->OpDetGeoFromOpDet(OpChannel).GetCenter(OpDetCenter);
1004  G4ThreeVector OpDetPoint(OpDetCenter[0]*CLHEP::cm,OpDetCenter[1]*CLHEP::cm,OpDetCenter[2]*CLHEP::cm);
1005  if (!Reflected) {
1006  double distance_in_cm = (x0 - OpDetPoint).mag()/CLHEP::cm; // this must be in CENTIMETERS!
1007  arrival_time_dist = GetVUVTime(distance_in_cm, NPhotons);//in ns
1008  }
1009  else {
1010  double t0_vis_lib = ReflT0s[OpChannel];
1011  arrival_time_dist = GetVisibleTimeOnlyCathode(t0_vis_lib, NPhotons);
1012  }
1013  }
1014 
1015  return arrival_time_dist;
1016  }
geo::GeometryCore::OpDetGeoFromOpDet
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Number of OpDets in the whole detector.
Definition: GeometryCore.cxx:1922
geo::OpDetGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: OpDetGeo.cxx:47
larg4::OpFastScintillation::GetVisibleTimeOnlyCathode
std::vector< double > GetVisibleTimeOnlyCathode(double, int) const
Definition: OpFastScintillation.cxx:1136
larg4::OpFastScintillation::GetVUVTime
std::vector< double > GetVUVTime(double, int) const
Definition: OpFastScintillation.cxx:1066
geo
Namespace collecting geometry-related classes utilities.
Definition: TrackContainmentAlg.hh:20
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool larg4::OpFastScintillation::RecordPhotonsProduced ( const G4Step &  aStep,
double  N 
)
protected

Definition at line 392 of file OpFastScintillation.cxx.

References larg4::OpDetPhotonTable::AddLitePhoton(), larg4::OpDetPhotonTable::AddOpDetBacktrackerRecord(), larg4::OpDetPhotonTable::AddPhoton(), sim::OpDetBacktrackerRecord::AddScintillationPhotons(), average_position(), bPropagate, sim::OnePhoton::Energy, ExcitationRatio, fFiniteRiseTime, sim::LArG4Parameters::FillSimEnergyDeposits(), phot::PhotonVisibilityService::GetAllVisibilities(), larg4::ParticleListAction::GetCurrentTrackID(), phot::PhotonVisibilityService::GetReflT0s(), phot::PhotonVisibilityService::GetTimingTF1(), phot::PhotonVisibilityService::IncludeParPropTime(), sim::OnePhoton::InitialPosition, larg4::IonizationAndScintillation::Instance(), larg4::OpDetPhotonTable::Instance(), min, sim::OnePhoton::MotherTrackID, phot::PhotonVisibilityService::NOpChannels(), ParPropTimeTF1, ProcessStep(), propagation_time(), reemission_energy(), ReflT0s, scint_time(), scintillationByParticleType, sim::OnePhoton::SetInSD, phot::PhotonVisibilityService::StoreReflected(), phot::PhotonVisibilityService::StoreReflT0(), theFastIntegralTable, theSlowIntegralTable, sim::OnePhoton::Time, sim::LArG4Parameters::UseLitePhotons(), and larg4::IonizationAndScintillation::VisibleEnergyDeposit().

Referenced by GetScintillationByParticleType(), and PostStepDoIt().

393  {
394  // make sure that whatever happens afterwards, the energy deposition is stored
395  art::ServiceHandle<sim::LArG4Parameters> lgp;
396  if(lgp->FillSimEnergyDeposits())
397  ProcessStep(aStep);
398 
399 
400  // Get the pointer to the fast scintillation table
401  OpDetPhotonTable * fst = OpDetPhotonTable::Instance();
402 
403  const G4Track * aTrack = aStep.GetTrack();
404 
405  G4StepPoint const* pPreStepPoint = aStep.GetPreStepPoint();
406  // unused G4StepPoint const* pPostStepPoint = aStep.GetPostStepPoint();
407 
408  const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
409  const G4Material* aMaterial = aTrack->GetMaterial();
410 
411  G4int materialIndex = aMaterial->GetIndex();
412  G4int tracknumber = aStep.GetTrack()->GetTrackID();
413 
414  G4ThreeVector x0 = pPreStepPoint->GetPosition();
415  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
416  //G4double t0 = pPreStepPoint->GetGlobalTime() - fGlobalTimeOffset;
417  G4double t0 = pPreStepPoint->GetGlobalTime();
418 
419 
420  G4MaterialPropertiesTable* aMaterialPropertiesTable =
421  aMaterial->GetMaterialPropertiesTable();
422 
423  G4MaterialPropertyVector* Fast_Intensity =
424  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
425  G4MaterialPropertyVector* Slow_Intensity =
426  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
427 
428  if (!Fast_Intensity && !Slow_Intensity )
429  return 1;
430 
431  if(!bPropagate)
432  return 0;
433 
434  // Get the visibility vector for this point
435  art::ServiceHandle<phot::PhotonVisibilityService> pvs;
436  size_t const NOpChannels = pvs->NOpChannels();
437 
438 
439  G4int nscnt = 1;
440  if (Fast_Intensity && Slow_Intensity) nscnt = 2;
441 
442 
443  double Num = 0;
444  double YieldRatio=0;
445 
446 
447  if (scintillationByParticleType) {
448  // The scintillation response is a function of the energy
449  // deposited by particle types.
450 
451  // Get the definition of the current particle
452  G4ParticleDefinition *pDef = aParticle->GetDefinition();
453 
454  // Obtain the G4MaterialPropertyVectory containing the
455  // scintillation light yield as a function of the deposited
456  // energy for the current particle type
457 
458  // Protons
459  if(pDef==G4Proton::ProtonDefinition())
460  {
461  YieldRatio = aMaterialPropertiesTable->
462  GetConstProperty("PROTONYIELDRATIO");
463 
464  }
465 
466  // Muons
467  else if(pDef==G4MuonPlus::MuonPlusDefinition()||pDef==G4MuonMinus::MuonMinusDefinition())
468  {
469  YieldRatio = aMaterialPropertiesTable->
470  GetConstProperty("MUONYIELDRATIO");
471  }
472 
473  // Pions
474  else if(pDef==G4PionPlus::PionPlusDefinition()||pDef==G4PionMinus::PionMinusDefinition())
475  {
476  YieldRatio = aMaterialPropertiesTable->
477  GetConstProperty("PIONYIELDRATIO");
478  }
479 
480  // Kaons
481  else if(pDef==G4KaonPlus::KaonPlusDefinition()||pDef==G4KaonMinus::KaonMinusDefinition())
482  {
483  YieldRatio = aMaterialPropertiesTable->
484  GetConstProperty("KAONYIELDRATIO");
485  }
486 
487  // Alphas
488  else if(pDef==G4Alpha::AlphaDefinition())
489  {
490  YieldRatio = aMaterialPropertiesTable->
491  GetConstProperty("ALPHAYIELDRATIO");
492  }
493 
494  // Electrons (must also account for shell-binding energy
495  // attributed to gamma from standard PhotoElectricEffect)
496  else if(pDef==G4Electron::ElectronDefinition() ||
497  pDef==G4Gamma::GammaDefinition())
498  {
499  YieldRatio = aMaterialPropertiesTable->
500  GetConstProperty("ELECTRONYIELDRATIO");
501  }
502 
503  // Default for particles not enumerated/listed above
504  else
505  {
506  YieldRatio = aMaterialPropertiesTable->
507  GetConstProperty("ELECTRONYIELDRATIO");
508  }
509 
510  // If the user has not specified yields for (p,d,t,a,carbon)
511  // then these unspecified particles will default to the
512  // electron's scintillation yield
513  if(YieldRatio==0){
514  {
515 
516  YieldRatio = aMaterialPropertiesTable->
517  GetConstProperty("ELECTRONYIELDRATIO");
518 
519  }
520  }
521  }
522 
523  double const xyz[3] = { x0[0]/CLHEP::cm, x0[1]/CLHEP::cm, x0[2]/CLHEP::cm };
524  float const* Visibilities = pvs->GetAllVisibilities(xyz);
525 
526  float const* ReflVisibilities = nullptr;
527 
528 
529  // Store timing information in the object for use in propagation_time method
530  if(pvs->StoreReflected()) {
531  ReflVisibilities = pvs->GetAllVisibilities(xyz,true);
532  if(pvs->StoreReflT0())
533  ReflT0s = pvs->GetReflT0s(xyz);
534  }
535  if(pvs->IncludeParPropTime())
536  {
537  ParPropTimeTF1 = pvs->GetTimingTF1(xyz);
538  }
539 
540  /*
541  // For Kazu to debug # photons generated using csv file, by default should be commented out
542  // but don't remove as it's useful. Separated portion of code relevant to this block
543  // is labeled as "CASE-DEBUG DO NOT REMOVE THIS COMMENT"
544  // (2017 May 2nd ... "kazuhiro at nevis dot columbia dot edu")
545  //
546  static bool first_time=false;
547  if(!first_time) {
548  std::cout<<"LOGMEid,pdg,mass,ke,te,de,x,y,z,t,dr,mean_npe,gen_npe,det_npe"<<std::endl;
549  first_time=true;
550  }
551 
552  std::cout<<"LOGME"
553  <<aStep.GetTrack()->GetTrackID()<<","
554  <<aParticle->GetDefinition()->GetPDGEncoding()<<","
555  <<aParticle->GetDefinition()->GetPDGMass()/CLHEP::MeV<<","
556  <<pPreStepPoint->GetKineticEnergy()<<","
557  <<pPreStepPoint->GetTotalEnergy()<<","
558  <<aStep.GetTotalEnergyDeposit()<<","
559  <<xyz[0]<<","<<xyz[1]<<","<<xyz[2]<<","<<t0<<","
560  <<aStep.GetDeltaPosition().mag()<<","
561  <<MeanNumberOfPhotons<<","<<std::flush;
562 
563  double gen_photon_ctr=0;
564  double det_photon_ctr=0;
565  */
566  for (G4int scnt = 1; scnt <= nscnt; scnt++) {
567 
568  G4double ScintillationTime = 0.*CLHEP::ns;
569  G4double ScintillationRiseTime = 0.*CLHEP::ns;
570  G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
571 
572  if (scnt == 1) {
573  if (nscnt == 1) {
574  if(Fast_Intensity){
575  ScintillationTime = aMaterialPropertiesTable->
576  GetConstProperty("FASTTIMECONSTANT");
577  if (fFiniteRiseTime) {
578  ScintillationRiseTime = aMaterialPropertiesTable->
579  GetConstProperty("FASTSCINTILLATIONRISETIME");
580  }
581  ScintillationIntegral =
582  (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
583  }
584  if(Slow_Intensity){
585  ScintillationTime = aMaterialPropertiesTable->
586  GetConstProperty("SLOWTIMECONSTANT");
587  if (fFiniteRiseTime) {
588  ScintillationRiseTime = aMaterialPropertiesTable->
589  GetConstProperty("SLOWSCINTILLATIONRISETIME");
590  }
591  ScintillationIntegral =
592  (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
593  }
594  }//endif nscnt=1
595  else {
596  if(YieldRatio==0)
597  YieldRatio = aMaterialPropertiesTable->
598  GetConstProperty("YIELDRATIO");
599 
600 
601  if ( ExcitationRatio == 1.0 ) {
602  Num = std::min(YieldRatio,1.0)*MeanNumberOfPhotons;
603  }
604  else {
605  Num = std::min(ExcitationRatio,1.0)*MeanNumberOfPhotons;
606  }
607  ScintillationTime = aMaterialPropertiesTable->
608  GetConstProperty("FASTTIMECONSTANT");
609  if (fFiniteRiseTime) {
610  ScintillationRiseTime = aMaterialPropertiesTable->
611  GetConstProperty("FASTSCINTILLATIONRISETIME");
612  }
613  ScintillationIntegral =
614  (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
615  }//endif nscnt!=1
616  }//end scnt=1
617 
618  else {
619  Num = MeanNumberOfPhotons - Num;
620  ScintillationTime = aMaterialPropertiesTable->
621  GetConstProperty("SLOWTIMECONSTANT");
622  if (fFiniteRiseTime) {
623  ScintillationRiseTime = aMaterialPropertiesTable->
624  GetConstProperty("SLOWSCINTILLATIONRISETIME");
625  }
626  ScintillationIntegral =
627  (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
628  }
629 
630  if (!ScintillationIntegral) continue;
631 
632  //gen_photon_ctr += Num; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
633 
634  // Max Scintillation Integral
635 
636  // G4double CIImax = ScintillationIntegral->GetMaxValue();
637 
638 
639  //std::cout << "++++++++++++" << Num << "++++++++++" << std::endl;
640 
641 
642  // here we go: now if visibilities are invalid, we are in trouble
643  //if (!Visibilities && (NOpChannels > 0)) {
644  // throw cet::exception("OpFastScintillator")
645  // << "Photon library does not cover point ( " << xyz[0] << ", "
646  // << xyz[1] << ", " << xyz[2] << " ) cm.\n";
647  //}
648 
649  if(!Visibilities){
650  }else{
651  std::map<int, int> DetectedNum;
652 
653  std::map<int, int> ReflDetectedNum;
654 
655  for(size_t OpDet=0; OpDet!=NOpChannels; OpDet++)
656  {
657  G4int DetThisPMT = G4int(G4Poisson(Visibilities[OpDet] * Num));
658 
659  if(DetThisPMT>0)
660  {
661  DetectedNum[OpDet]=DetThisPMT;
662  // mf::LogInfo("OpFastScintillation") << "FastScint: " <<
663  // // it->second<<" " << Num << " " << DetThisPMT;
664 
665  //det_photon_ctr += DetThisPMT; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
666  }
667  if(pvs->StoreReflected()) {
668  G4int ReflDetThisPMT = G4int(G4Poisson(ReflVisibilities[OpDet] * Num));
669  if(ReflDetThisPMT>0)
670  {
671  ReflDetectedNum[OpDet]=ReflDetThisPMT;
672  }
673  }
674 
675  }
676 
677 
678  // Now we run through each PMT figuring out num of detected photons
679  for (int Reflected = 0; Reflected <= 1; Reflected++) {
680  // Only do the reflected loop if we have reflected visibilities
681  if (Reflected && !pvs->StoreReflected())
682  continue;
683 
684  std::map<int,int>::const_iterator itstart;
685  std::map<int,int>::const_iterator itend;
686  if (Reflected) {
687  itstart = ReflDetectedNum.begin();
688  itend = ReflDetectedNum.end();
689  }
690  else{
691  itstart = DetectedNum.begin();
692  itend = DetectedNum.end();
693  }
694 
695  for(std::map<int,int>::const_iterator itdetphot=itstart; itdetphot!=itend; ++itdetphot)
696  {
697  int OpChannel = itdetphot->first;
698  int NPhotons = itdetphot->second;
699 
700  // Set up the OpDetBTR information
701  sim::OpDetBacktrackerRecord tmpOpDetBTRecord(OpChannel);
702  int thisG4TrackID = ParticleListAction::GetCurrentTrackID();
703  double xyzPos[3];
704  average_position(aStep, xyzPos);
705  double Edeposited = larg4::IonizationAndScintillation::Instance()->VisibleEnergyDeposit()/CLHEP::MeV;
706 
707  // Get the transport time distribution
708  std::vector<double> arrival_time_dist = propagation_time(x0, OpChannel, NPhotons, Reflected);
709 
710  // Loop through the photons
711  for (G4int i = 0; i < NPhotons; ++i)
712  {
713  G4double Time = t0
714  + scint_time(aStep, ScintillationTime, ScintillationRiseTime)
715  + arrival_time_dist[i]*CLHEP::ns;
716 
717  // Always store the BTR
718  tmpOpDetBTRecord.AddScintillationPhotons(thisG4TrackID, Time, 1, xyzPos, Edeposited);
719 
720  // Store as lite photon or as OnePhoton
721  if(lgp->UseLitePhotons())
722  {
723  fst->AddLitePhoton(OpChannel, static_cast<int>(Time), 1, Reflected);
724  }
725  else {
726  // The sim photon in this case stores its production point and time
727  TVector3 PhotonPosition( x0[0], x0[1], x0[2] );
728 
729  float PhotonEnergy = 0;
730  if (Reflected) PhotonEnergy = reemission_energy()*CLHEP::eV;
731  else PhotonEnergy = 9.7*CLHEP::eV;
732 
733  // Make a photon object for the collection
734  sim::OnePhoton PhotToAdd;
735  PhotToAdd.InitialPosition = PhotonPosition;
736  PhotToAdd.Energy = PhotonEnergy;
737  PhotToAdd.Time = Time;
738  PhotToAdd.SetInSD = false;
739  PhotToAdd.MotherTrackID = tracknumber;
740 
741  fst->AddPhoton(OpChannel, std::move(PhotToAdd), Reflected);
742  }
743  }
744  fst->AddOpDetBacktrackerRecord(tmpOpDetBTRecord, Reflected);
745  }
746  }
747  }
748  }
749 
750  //std::cout<<gen_photon_ctr<<","<<det_photon_ctr<<std::endl; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
751  return 0;
752  }
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
larg4::OpFastScintillation::average_position
void average_position(G4Step const &aStep, double *xzyPos) const
Definition: OpFastScintillation.cxx:1049
larg4::OpFastScintillation::propagation_time
std::vector< double > propagation_time(G4ThreeVector x0, int OpChannel, int NPhotons, bool Reflected=false) const
Definition: OpFastScintillation.cxx:971
larg4::OpFastScintillation::ProcessStep
void ProcessStep(const G4Step &step)
Definition: OpFastScintillation.cxx:367
larg4::OpFastScintillation::scint_time
G4double scint_time(const G4Step &aStep, G4double ScintillationTime, G4double ScintillationRiseTime) const
Definition: OpFastScintillation.cxx:950
sim::OpDetBacktrackerRecord
Energy deposited on a readout Optical Detector by simulated tracks.
Definition: OpDetBacktrackerRecord.h:127
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:25
phot::PhotonVisibilityService::GetAllVisibilities
float const * GetAllVisibilities(double const *xyz, bool wantReflected=false) const
Definition: PhotonVisibilityService_service.cc:355
phot::PhotonVisibilityService::GetTimingTF1
TF1 * GetTimingTF1(double const *xyz) const
Definition: PhotonVisibilityService_service.cc:538
larg4::IonizationAndScintillation::Instance
static IonizationAndScintillation * Instance()
Definition: IonizationAndScintillation.cxx:38
sim::LArG4Parameters::FillSimEnergyDeposits
bool FillSimEnergyDeposits() const
Definition: LArG4Parameters.h:67
larg4::OpDetPhotonTable::Instance
static OpDetPhotonTable * Instance(bool LitePhotons=false)
Definition: OpDetPhotonTable.cxx:36
larg4::OpFastScintillation::bPropagate
bool bPropagate
Whether propagation of photons is enabled.
Definition: OpFastScintillation.hh:299
sim::OnePhoton::InitialPosition
TVector3 InitialPosition
Definition: SimPhotons.h:49
min
Int_t min
Definition: plot.C:26
sim::LArG4Parameters::UseLitePhotons
bool UseLitePhotons() const
Definition: LArG4Parameters.h:65
phot::PhotonVisibilityService::GetReflT0s
float const * GetReflT0s(double const *xyz) const
Definition: PhotonVisibilityService_service.cc:488
double larg4::OpFastScintillation::reemission_energy ( ) const
private

Definition at line 1043 of file OpFastScintillation.cxx.

References rgen0, and tpbemission.

Referenced by RecordPhotonsProduced().

1044  {
1045  return rgen0->fire()*((*(--tpbemission.end())).first-(*tpbemission.begin()).first)+(*tpbemission.begin()).first;
1046  }
larg4::OpFastScintillation::tpbemission
std::map< double, double > tpbemission
Definition: OpFastScintillation.hh:279
void larg4::OpFastScintillation::RemoveSaturation ( )
inline

Definition at line 212 of file OpFastScintillation.hh.

References emSaturation.

Referenced by SetScintillationByParticleType().

212 { emSaturation = NULL; }
G4double larg4::OpFastScintillation::sample_time ( G4double  tau1,
G4double  tau2 
) const
private

Definition at line 1021 of file OpFastScintillation.cxx.

References bi_exp(), d, and single_exp().

Referenced by scint_time().

1022  {
1023 // tau1: rise time and tau2: decay time
1024 
1025  while(1) {
1026  // two random numbers
1027  G4double ran1 = G4UniformRand();
1028  G4double ran2 = G4UniformRand();
1029  //
1030  // exponential distribution as envelope function: very efficient
1031  //
1032  G4double d = (tau1+tau2)/tau2;
1033  // make sure the envelope function is
1034  // always larger than the bi-exponential
1035  G4double t = -1.0*tau2*std::log(1-ran1);
1036  G4double g = d*single_exp(t,tau2);
1037  if (ran2 <= bi_exp(t,tau1,tau2)/g) return t;
1038  }
1039  return -1.0;
1040  }
larg4::OpFastScintillation::bi_exp
G4double bi_exp(G4double t, G4double tau1, G4double tau2) const
Definition: OpFastScintillation.hh:412
larg4::OpFastScintillation::single_exp
G4double single_exp(G4double t, G4double tau2) const
Definition: OpFastScintillation.hh:406
d
Float_t d
Definition: plot.C:237
G4double larg4::OpFastScintillation::scint_time ( const G4Step &  aStep,
G4double  ScintillationTime,
G4double  ScintillationRiseTime 
) const
private

Definition at line 950 of file OpFastScintillation.cxx.

References sample_time().

Referenced by RecordPhotonsProduced().

953  {
954  G4StepPoint const* pPreStepPoint = aStep.GetPreStepPoint();
955  G4StepPoint const* pPostStepPoint = aStep.GetPostStepPoint();
956  G4double avgVelocity = (pPreStepPoint->GetVelocity() + pPostStepPoint->GetVelocity())/2.;
957 
958  G4double deltaTime = aStep.GetStepLength() / avgVelocity;
959 
960  if (ScintillationRiseTime==0.0) {
961  deltaTime = deltaTime -
962  ScintillationTime * std::log( G4UniformRand() );
963  } else {
964  deltaTime = deltaTime +
965  sample_time(ScintillationRiseTime, ScintillationTime);
966  }
967  return deltaTime;
968  }
larg4::OpFastScintillation::sample_time
G4double sample_time(G4double tau1, G4double tau2) const
Definition: OpFastScintillation.cxx:1021
void larg4::OpFastScintillation::SetFiniteRiseTime ( const G4bool  state)
inline

Definition at line 326 of file OpFastScintillation.hh.

References fFiniteRiseTime.

327 {
328  fFiniteRiseTime = state;
329 }
void larg4::OpFastScintillation::SetScintillationByParticleType ( const G4bool  scintType)

Definition at line 913 of file OpFastScintillation.cxx.

References emSaturation, RemoveSaturation(), and scintillationByParticleType.

Referenced by GetSaturation().

914  {
915  if (emSaturation) {
916  G4Exception("OpFastScintillation::SetScintillationByParticleType", "Scint02",
917  JustWarning, "Redefinition: Birks Saturation is replaced by ScintillationByParticleType!");
918  RemoveSaturation();
919  }
920  scintillationByParticleType = scintType;
921  }
void larg4::OpFastScintillation::SetScintillationExcitationRatio ( const G4double  excitationratio)
inline

Definition at line 356 of file OpFastScintillation.hh.

References ExcitationRatio.

357 {
358  ExcitationRatio = excitationratio;
359 }
void larg4::OpFastScintillation::SetScintillationYieldFactor ( const G4double  yieldfactor)
inline

Definition at line 344 of file OpFastScintillation.hh.

References YieldFactor.

345 {
346  YieldFactor = yieldfactor;
347 }
void larg4::OpFastScintillation::SetTrackSecondariesFirst ( const G4bool  state)
inline

Definition at line 320 of file OpFastScintillation.hh.

References fTrackSecondariesFirst.

321 {
322  fTrackSecondariesFirst = state;
323 }
G4double larg4::OpFastScintillation::single_exp ( G4double  t,
G4double  tau2 
) const
inlineprivate

Definition at line 406 of file OpFastScintillation.hh.

Referenced by sample_time().

407 {
408  return std::exp(-1.0*t/tau2)/tau2;
409 }

Member Data Documentation

bool larg4::OpFastScintillation::bPropagate
private

Whether propagation of photons is enabled.

Definition at line 299 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

G4EmSaturation* larg4::OpFastScintillation::emSaturation
private

Definition at line 284 of file OpFastScintillation.hh.

Referenced by AddSaturation(), GetSaturation(), OpFastScintillation(), RemoveSaturation(), and SetScintillationByParticleType().

G4double larg4::OpFastScintillation::ExcitationRatio
protected

Definition at line 260 of file OpFastScintillation.hh.

Referenced by GetScintillationExcitationRatio(), OpFastScintillation(), RecordPhotonsProduced(), and SetScintillationExcitationRatio().

double larg4::OpFastScintillation::fd_break
private

Definition at line 291 of file OpFastScintillation.hh.

Referenced by GetVUVTime(), and OpFastScintillation().

double larg4::OpFastScintillation::fd_max
private

Definition at line 292 of file OpFastScintillation.hh.

Referenced by GetVUVTime(), and OpFastScintillation().

G4bool larg4::OpFastScintillation::fFiniteRiseTime
protected

Definition at line 256 of file OpFastScintillation.hh.

Referenced by GetFiniteRiseTime(), OpFastScintillation(), RecordPhotonsProduced(), and SetFiniteRiseTime().

double larg4::OpFastScintillation::ft0_break_point
private

Definition at line 294 of file OpFastScintillation.hh.

Referenced by GetVisibleTimeOnlyCathode(), and OpFastScintillation().

double larg4::OpFastScintillation::ft0_max
private

Definition at line 294 of file OpFastScintillation.hh.

Referenced by GetVisibleTimeOnlyCathode(), and OpFastScintillation().

double larg4::OpFastScintillation::ftf1_sampling_factor
private

Definition at line 293 of file OpFastScintillation.hh.

Referenced by GetVisibleTimeOnlyCathode(), GetVUVTime(), and OpFastScintillation().

G4bool larg4::OpFastScintillation::fTrackSecondariesFirst
protected

Definition at line 255 of file OpFastScintillation.hh.

Referenced by GetTrackSecondariesFirst(), OpFastScintillation(), and SetTrackSecondariesFirst().

TF1 const* larg4::OpFastScintillation::functions_vis[5]
private

Definition at line 290 of file OpFastScintillation.hh.

Referenced by GetVisibleTimeOnlyCathode(), and OpFastScintillation().

TF1 const* larg4::OpFastScintillation::functions_vuv[8]
private

Definition at line 289 of file OpFastScintillation.hh.

Referenced by GetVUVTime(), and OpFastScintillation().

TF1* larg4::OpFastScintillation::ParPropTimeTF1
private

Definition at line 286 of file OpFastScintillation.hh.

Referenced by propagation_time(), and RecordPhotonsProduced().

float const* larg4::OpFastScintillation::ReflT0s
private

Definition at line 287 of file OpFastScintillation.hh.

Referenced by propagation_time(), and RecordPhotonsProduced().

CLHEP::RandGeneral* larg4::OpFastScintillation::rgen0
private

Definition at line 280 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and reemission_energy().

G4bool larg4::OpFastScintillation::scintillationByParticleType
protected

Definition at line 262 of file OpFastScintillation.hh.

Referenced by GetScintillationByParticleType(), OpFastScintillation(), RecordPhotonsProduced(), and SetScintillationByParticleType().

G4PhysicsTable* larg4::OpFastScintillation::theFastIntegralTable
protected

Definition at line 253 of file OpFastScintillation.hh.

Referenced by BuildThePhysicsTable(), DumpPhysicsTable(), GetFastIntegralTable(), OpFastScintillation(), RecordPhotonsProduced(), and ~OpFastScintillation().

G4PhysicsTable* larg4::OpFastScintillation::theSlowIntegralTable
protected

Definition at line 252 of file OpFastScintillation.hh.

Referenced by BuildThePhysicsTable(), DumpPhysicsTable(), GetSlowIntegralTable(), OpFastScintillation(), RecordPhotonsProduced(), and ~OpFastScintillation().

std::map<double,double> larg4::OpFastScintillation::tpbemission
private

Definition at line 279 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and reemission_energy().

G4double larg4::OpFastScintillation::YieldFactor
protected

Definition at line 258 of file OpFastScintillation.hh.

Referenced by GetScintillationYieldFactor(), OpFastScintillation(), and SetScintillationYieldFactor().

The documentation for this class was generated from the following files: