LArSoft  v06_85_00
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)
 
std::vector< double > GetVisibleTimeOnlyCathode (double, int)
 

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)
 
G4double bi_exp (G4double t, G4double tau1, G4double tau2)
 
G4double sample_time (G4double tau1, G4double tau2)
 
void ProcessStep (const G4Step &step)
 

Private Attributes

G4EmSaturation * emSaturation
 
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
 

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 BuildThePhysicsTable(), emSaturation, ExcitationRatio, fd_break, fd_max, fFiniteRiseTime, ft0_break_point, ft0_max, ftf1_sampling_factor, fTrackSecondariesFirst, functions_vis, functions_vuv, phot::PhotonVisibilityService::IncludePropTime(), detinfo::LArProperties::ScintByParticleType(), scintillationByParticleType, phot::PhotonVisibilityService::SetDirectLightPropFunctions(), phot::PhotonVisibilityService::SetReflectedCOLightPropFunctions(), theFastIntegralTable, theSlowIntegralTable, and YieldFactor.

152  : G4VRestDiscreteProcess(processName, type)
153  {
154  G4cout<<"BEA1 timing distribution"<<G4endl;
155  SetProcessSubType(25);
156 
157  fTrackSecondariesFirst = false;
158  fFiniteRiseTime = false;
159 
160 
161  YieldFactor=1.0;
162  ExcitationRatio = 1.0;
163 
164  const detinfo::LArProperties* larp = lar::providerFrom<detinfo::LArPropertiesService>();
165  art::ServiceHandle<sim::LArG4Parameters> lgp;
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  gRandom->SetSeed(0);
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 }
phot::PhotonVisibilityService::SetDirectLightPropFunctions
void SetDirectLightPropFunctions(TF1 const *functions[8], double &d_break, double &d_max, double &tf1_sampling_factor) const
Definition: PhotonVisibilityService_service.cc:570
phot::PhotonVisibilityService::SetReflectedCOLightPropFunctions
void SetReflectedCOLightPropFunctions(TF1 const *functions[5], double &t0_max, double &t0_break_point) const
Definition: PhotonVisibilityService_service.cc:587
detinfo::LArProperties::ScintByParticleType
virtual bool ScintByParticleType() const =0
larg4::OpFastScintillation::OpFastScintillation ( const OpFastScintillation &  right)

Definition at line 188 of file OpFastScintillation.cxx.

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

189  : G4VRestDiscreteProcess(rhs.GetProcessName(), rhs.GetProcessType())
190  {
191  theSlowIntegralTable = rhs.GetSlowIntegralTable();
192  theFastIntegralTable = rhs.GetFastIntegralTable();
193 
194  fTrackSecondariesFirst = rhs.GetTrackSecondariesFirst();
195  fFiniteRiseTime = rhs.GetFiniteRiseTime();
196  YieldFactor = rhs.GetScintillationYieldFactor();
197  ExcitationRatio = rhs.GetScintillationExcitationRatio();
198  scintillationByParticleType = rhs.GetScintillationByParticleType();
199  emSaturation = rhs.GetSaturation();
200 
201  BuildThePhysicsTable();
202  }
larg4::OpFastScintillation::~OpFastScintillation ( )

Definition at line 209 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

210 {
211  if (theFastIntegralTable != NULL) {
212  theFastIntegralTable->clearAndDestroy();
213  delete theFastIntegralTable;
214  }
215  if (theSlowIntegralTable != NULL) {
216  theSlowIntegralTable->clearAndDestroy();
217  delete theSlowIntegralTable;
218  }
219 
220 }

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 230 of file OpFastScintillation.cxx.

References PostStepDoIt().

235 {
236  return OpFastScintillation::PostStepDoIt(aTrack, aStep);
237 }
larg4::OpFastScintillation::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: OpFastScintillation.cxx:243
G4double larg4::OpFastScintillation::bi_exp ( G4double  t,
G4double  tau1,
G4double  tau2 
)
inlineprivate

Definition at line 395 of file OpFastScintillation.hh.

Referenced by sample_time().

396 {
397  return std::exp(-1.0*t/tau2)*(1-std::exp(-1.0*t/tau1))/tau2/tau2*(tau1+tau2);
398 }
void larg4::OpFastScintillation::BuildThePhysicsTable ( )
protected

Definition at line 878 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

Referenced by GetScintillationByParticleType(), and OpFastScintillation().

879 {
880  if (theFastIntegralTable && theSlowIntegralTable) return;
881 
882  const G4MaterialTable* theMaterialTable =
883  G4Material::GetMaterialTable();
884  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
885 
886  // create new physics table
887 
888  if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
889  if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
890 
891  // loop for materials
892 
893  for (G4int i=0 ; i < numOfMaterials; i++)
894  {
895  G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
896  new G4PhysicsOrderedFreeVector();
897  G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
898  new G4PhysicsOrderedFreeVector();
899 
900  // Retrieve vector of scintillation wavelength intensity for
901  // the material from the material's optical properties table.
902 
903  G4Material* aMaterial = (*theMaterialTable)[i];
904 
905  G4MaterialPropertiesTable* aMaterialPropertiesTable =
906  aMaterial->GetMaterialPropertiesTable();
907 
908  if (aMaterialPropertiesTable) {
909 
910  G4MaterialPropertyVector* theFastLightVector =
911  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
912 
913  if (theFastLightVector) {
914 
915  // Retrieve the first intensity point in vector
916  // of (photon energy, intensity) pairs
917 
918  G4double currentIN = (*theFastLightVector)[0];
919 
920  if (currentIN >= 0.0) {
921 
922  // Create first (photon energy, Scintillation
923  // Integral pair
924 
925  G4double currentPM = theFastLightVector->Energy(0);
926 
927  G4double currentCII = 0.0;
928 
929  aPhysicsOrderedFreeVector->
930  InsertValues(currentPM , currentCII);
931 
932  // Set previous values to current ones prior to loop
933 
934  G4double prevPM = currentPM;
935  G4double prevCII = currentCII;
936  G4double prevIN = currentIN;
937 
938  // loop over all (photon energy, intensity)
939  // pairs stored for this material
940 
941  for (size_t i = 1;
942  i < theFastLightVector->GetVectorLength();
943  i++)
944  {
945  currentPM = theFastLightVector->Energy(i);
946  currentIN = (*theFastLightVector)[i];
947 
948  currentCII = 0.5 * (prevIN + currentIN);
949 
950  currentCII = prevCII +
951  (currentPM - prevPM) * currentCII;
952 
953  aPhysicsOrderedFreeVector->
954  InsertValues(currentPM, currentCII);
955 
956  prevPM = currentPM;
957  prevCII = currentCII;
958  prevIN = currentIN;
959  }
960 
961  }
962  }
963 
964  G4MaterialPropertyVector* theSlowLightVector =
965  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
966 
967  if (theSlowLightVector) {
968 
969  // Retrieve the first intensity point in vector
970  // of (photon energy, intensity) pairs
971 
972  G4double currentIN = (*theSlowLightVector)[0];
973 
974  if (currentIN >= 0.0) {
975 
976  // Create first (photon energy, Scintillation
977  // Integral pair
978 
979  G4double currentPM = theSlowLightVector->Energy(0);
980 
981  G4double currentCII = 0.0;
982 
983  bPhysicsOrderedFreeVector->
984  InsertValues(currentPM , currentCII);
985 
986  // Set previous values to current ones prior to loop
987 
988  G4double prevPM = currentPM;
989  G4double prevCII = currentCII;
990  G4double prevIN = currentIN;
991 
992  // loop over all (photon energy, intensity)
993  // pairs stored for this material
994 
995  for (size_t i = 1;
996  i < theSlowLightVector->GetVectorLength();
997  i++)
998  {
999  currentPM = theSlowLightVector->Energy(i);
1000  currentIN = (*theSlowLightVector)[i];
1001 
1002  currentCII = 0.5 * (prevIN + currentIN);
1003 
1004  currentCII = prevCII +
1005  (currentPM - prevPM) * currentCII;
1006 
1007  bPhysicsOrderedFreeVector->
1008  InsertValues(currentPM, currentCII);
1009 
1010  prevPM = currentPM;
1011  prevCII = currentCII;
1012  prevIN = currentIN;
1013  }
1014 
1015  }
1016  }
1017  }
1018 
1019  // The scintillation integral(s) for a given material
1020  // will be inserted in the table(s) according to the
1021  // position of the material in the material table.
1022 
1023  theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
1024  theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
1025 
1026  }
1027 }
void larg4::OpFastScintillation::DumpPhysicsTable ( ) const
inline

Definition at line 363 of file OpFastScintillation.hh.

References theFastIntegralTable, and theSlowIntegralTable.

Referenced by GetScintillationByParticleType().

364 {
365  if (theFastIntegralTable) {
366  G4int PhysicsTableSize = theFastIntegralTable->entries();
367  G4PhysicsOrderedFreeVector *v;
368 
369  for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
370  {
371  v = (G4PhysicsOrderedFreeVector*)(*theFastIntegralTable)[i];
372  v->DumpValues();
373  }
374  }
375 
376  if (theSlowIntegralTable) {
377  G4int PhysicsTableSize = theSlowIntegralTable->entries();
378  G4PhysicsOrderedFreeVector *v;
379 
380  for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
381  {
382  v = (G4PhysicsOrderedFreeVector*)(*theSlowIntegralTable)[i];
383  v->DumpValues();
384  }
385  }
386 }
G4PhysicsTable * larg4::OpFastScintillation::GetFastIntegralTable ( ) const
inline

Definition at line 357 of file OpFastScintillation.hh.

References theFastIntegralTable.

Referenced by OpFastScintillation().

358 {
359  return theFastIntegralTable;
360 }
G4bool larg4::OpFastScintillation::GetFiniteRiseTime ( ) const
inline

Definition at line 321 of file OpFastScintillation.hh.

References fFiniteRiseTime.

Referenced by OpFastScintillation().

322 {
323  return fFiniteRiseTime;
324 }
G4double larg4::OpFastScintillation::GetMeanFreePath ( const G4Track &  aTrack,
G4double  ,
G4ForceCondition *  condition 
)

Definition at line 1046 of file OpFastScintillation.cxx.

1049 {
1050  *condition = StronglyForced;
1051 
1052  return DBL_MAX;
1053 
1054 }
G4double larg4::OpFastScintillation::GetMeanLifeTime ( const G4Track &  aTrack,
G4ForceCondition *  condition 
)

Definition at line 1060 of file OpFastScintillation.cxx.

1062 {
1063  *condition = Forced;
1064 
1065  return DBL_MAX;
1066 
1067 }
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 345 of file OpFastScintillation.hh.

References ExcitationRatio.

Referenced by OpFastScintillation().

346 {
347  return ExcitationRatio;
348 }
G4double larg4::OpFastScintillation::GetScintillationYieldFactor ( ) const
inline

Definition at line 333 of file OpFastScintillation.hh.

References YieldFactor.

Referenced by OpFastScintillation().

334 {
335  return YieldFactor;
336 }
G4PhysicsTable * larg4::OpFastScintillation::GetSlowIntegralTable ( ) const
inline

Definition at line 351 of file OpFastScintillation.hh.

References theSlowIntegralTable.

Referenced by OpFastScintillation().

352 {
353  return theSlowIntegralTable;
354 }
G4bool larg4::OpFastScintillation::GetTrackSecondariesFirst ( ) const
inline

Definition at line 315 of file OpFastScintillation.hh.

References fTrackSecondariesFirst.

Referenced by OpFastScintillation().

316 {
317  return fTrackSecondariesFirst;
318 }
std::vector< double > larg4::OpFastScintillation::GetVisibleTimeOnlyCathode ( double  t0,
int  number_photons 
)

Definition at line 1165 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 RecordPhotonsProduced().

1165  {
1166  //-----Distances in cm and times in ns-----//
1167  //gRandom->SetSeed(0);
1168 
1169  std::vector<double> arrival_time_distrb;
1170  arrival_time_distrb.clear();
1171  // Parametrization data:
1172 
1173  if(t0 < 8 || t0 > ft0_max) {
1174  G4cout<<"WARNING: Parametrization of Cathode-Only reflected Light not fully reliable"<<G4endl;
1175  G4cout<<"Too close/far to the PMT -> set 0 Visible photons(?)!!!!!!"<<G4endl;
1176  return arrival_time_distrb;
1177  }
1178  //signals (remember this is transportation) no longer than 1us
1179  const double signal_t_range = 1000.;
1180  double pars_landau[3] = {functions_vis[1]->Eval(t0), functions_vis[2]->Eval(t0),
1181  pow(10.,functions_vis[0]->Eval(t0))};
1182  double pars_expo[2] = {functions_vis[3]->Eval(t0), functions_vis[4]->Eval(t0)};
1183  if(t0 > ft0_break_point) {
1184  pars_landau[0] = -0.798934 + 1.06216*t0;
1185  pars_landau[1] = functions_vis[2]->Eval(ft0_break_point);
1186  pars_landau[2] = pow(10.,functions_vis[0]->Eval(ft0_break_point));
1187  pars_expo[0] = functions_vis[3]->Eval(ft0_break_point);
1188  pars_expo[1] = functions_vis[4]->Eval(ft0_break_point);
1189  }
1190 
1191  // Defining the two functions (Landau + Exponential) describing the timing vs t0
1192  TF1 *flandau = new TF1("flandau","[2]*TMath::Landau(x,[0],[1])",0,signal_t_range/2);
1193  flandau->SetParameters(pars_landau);
1194  TF1 *fexpo = new TF1("fexpo","expo",0, signal_t_range/2);
1195  fexpo->SetParameters(pars_expo);
1196  //this is to find the intersection point between the two functions:
1197  TF1 *fint = new TF1("fint",finter_d,flandau->GetMaximumX(),2*t0,5);
1198  double parsInt[5] = {pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1199  fint->SetParameters(parsInt);
1200  double t_int = fint->GetMinimumX();
1201  double minVal = fint->Eval(t_int);
1202  //the functions must intersect!!!
1203  if(minVal>0.015)
1204  G4cout<<"WARNING: Parametrization of Direct Light discontinuous (landau + expo)!!!!!!"<<G4endl;
1205 
1206  TF1 *fVisTiming = new TF1("fTiming",LandauPlusExpoFinal,0,signal_t_range,6);
1207  double parsfinal[6] = {t_int, pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1208  fVisTiming->SetParameters(parsfinal);
1209  // Set the number of points used to sample the function
1210 
1211  int f_sampling = 1000;
1212  if(t0 < 20)
1213  f_sampling *= ftf1_sampling_factor;
1214  fVisTiming->SetNpx(f_sampling);
1215 
1216  for(int i=0; i<number_photons; i++)
1217  arrival_time_distrb.push_back(fVisTiming->GetRandom());
1218 
1219  //deleting ...
1220 
1221  delete flandau;
1222  delete fexpo;
1223  delete fint;
1224  delete fVisTiming;
1225  G4cout<<"Timing distribution BEAMAUS!"<<G4endl;
1226  return arrival_time_distrb;
1227 }
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:1229
larg4::LandauPlusExpoFinal
double LandauPlusExpoFinal(double *x, double *par)
Definition: OpFastScintillation.cxx:1236
std::vector< double > larg4::OpFastScintillation::GetVUVTime ( double  distance,
int  number_photons 
)

Definition at line 1095 of file OpFastScintillation.cxx.

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

Referenced by GetScintillationByParticleType(), and RecordPhotonsProduced().

1095  {
1096 
1097  //-----Distances in cm and times in ns-----//
1098  //gRandom->SetSeed(0);
1099  std::vector<double> arrival_time_distrb;
1100  arrival_time_distrb.clear();
1101 
1102  // Parametrization data:
1103 
1104  if(distance < 10 || distance > fd_max) {
1105  G4cout<<"WARNING: Parametrization of Direct Light not fully reliable"<<G4endl;
1106  G4cout<<"Too close/far to the PMT -> set 0 VUV photons(?)!!!!!!"<<G4endl;
1107  return arrival_time_distrb;
1108  }
1109  //signals (remember this is transportation) no longer than 1us
1110  const double signal_t_range = 1000.;
1111  const double vuv_vgroup = 10.13;//cm/ns
1112  double t_direct = distance/vuv_vgroup;
1113  // Defining the two functions (Landau + Exponential) describing the timing vs distance
1114  double pars_landau[3] = {functions_vuv[1]->Eval(distance), functions_vuv[2]->Eval(distance),
1115  pow(10.,functions_vuv[0]->Eval(distance))};
1116  if(distance > fd_break) {
1117  pars_landau[0]=functions_vuv[6]->Eval(distance);
1118  pars_landau[1]=functions_vuv[2]->Eval(fd_break);
1119  pars_landau[2]=pow(10.,functions_vuv[5]->Eval(distance));
1120  }
1121  TF1 *flandau = new TF1("flandau","[2]*TMath::Landau(x,[0],[1])",0,signal_t_range/2);
1122  flandau->SetParameters(pars_landau);
1123  double pars_expo[2] = {functions_vuv[3]->Eval(distance), functions_vuv[4]->Eval(distance)};
1124  if(distance > (fd_break - 50.)) {
1125  pars_expo[0] = functions_vuv[7]->Eval(distance);
1126  pars_expo[1] = functions_vuv[4]->Eval(fd_break - 50.);
1127  }
1128  TF1 *fexpo = new TF1("fexpo","expo",0, signal_t_range/2);
1129  fexpo->SetParameters(pars_expo);
1130  //this is to find the intersection point between the two functions:
1131  TF1 *fint = new TF1("fint",finter_d,flandau->GetMaximumX(),3*t_direct,5);
1132  double parsInt[5] = {pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1133  fint->SetParameters(parsInt);
1134  double t_int = fint->GetMinimumX();
1135  double minVal = fint->Eval(t_int);
1136  //the functions must intersect!!!
1137  if(minVal>0.015)
1138  G4cout<<"WARNING: Parametrization of Direct Light discontinuous (landau + expo)!!!!!!"<<G4endl;
1139 
1140  TF1 *fVUVTiming = new TF1("fTiming",LandauPlusExpoFinal,0,signal_t_range,6);
1141  double parsfinal[6] = {t_int, pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
1142  fVUVTiming->SetParameters(parsfinal);
1143  // Set the number of points used to sample the function
1144 
1145  int f_sampling = 1000;
1146  if(distance < 50)
1147  f_sampling *= ftf1_sampling_factor;
1148  fVUVTiming->SetNpx(f_sampling);
1149 
1150  for(int i=0; i<number_photons; i++)
1151  arrival_time_distrb.push_back(fVUVTiming->GetRandom());
1152 
1153  //deleting ...
1154  delete flandau;
1155  delete fexpo;
1156  delete fint;
1157  delete fVUVTiming;
1158  G4cout<<"BEAMAUS timing distribution hecha"<<G4endl;
1159  return arrival_time_distrb;
1160 }
larg4::finter_d
double finter_d(double *x, double *par)
Definition: OpFastScintillation.cxx:1229
larg4::LandauPlusExpoFinal
double LandauPlusExpoFinal(double *x, double *par)
Definition: OpFastScintillation.cxx:1236
G4bool larg4::OpFastScintillation::IsApplicable ( const G4ParticleDefinition &  aParticleType)
inlinevirtual

Definition at line 294 of file OpFastScintillation.hh.

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

295 {
296  if (aParticleType.GetParticleName() == "opticalphoton") return false;
297  if (aParticleType.IsShortLived()) return false;
298 
299  return true;
300 }
G4VParticleChange * larg4::OpFastScintillation::PostStepDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)
virtual

Definition at line 243 of file OpFastScintillation.cxx.

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

Referenced by AtRestDoIt().

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

Definition at line 353 of file OpFastScintillation.cxx.

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

Referenced by RecordPhotonsProduced().

354  {
355  if(step.GetTotalEnergyDeposit() <= 0) return;
356 
357  OpDetPhotonTable::Instance()->AddEnergyDeposit
358  (-1,
359  -1,
360  (double)(step.GetTotalEnergyDeposit()/CLHEP::MeV), //energy in MeV
361  (float)(step.GetPreStepPoint()->GetPosition().x()/CLHEP::cm),
362  (float)(step.GetPreStepPoint()->GetPosition().y()/CLHEP::cm),
363  (float)(step.GetPreStepPoint()->GetPosition().z()/CLHEP::cm),
364  (float)(step.GetPostStepPoint()->GetPosition().x()/CLHEP::cm),
365  (float)(step.GetPostStepPoint()->GetPosition().y()/CLHEP::cm),
366  (float)(step.GetPostStepPoint()->GetPosition().z()/CLHEP::cm),
367  (double)(step.GetPreStepPoint()->GetGlobalTime()),
368  (double)(step.GetPostStepPoint()->GetGlobalTime()),
369  //step.GetTrack()->GetTrackID(),
370  ParticleListAction::GetCurrentTrackID(),
371  step.GetTrack()->GetParticleDefinition()->GetPDGEncoding(),
372  step.GetPreStepPoint()->GetPhysicalVolume()->GetName()
373  );
374  }
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 vol="EMPTY")
Definition: OpDetPhotonTable.cxx:154
larg4::OpDetPhotonTable::Instance
static OpDetPhotonTable * Instance(bool LitePhotons=false)
Definition: OpDetPhotonTable.cxx:36
bool larg4::OpFastScintillation::RecordPhotonsProduced ( const G4Step &  aStep,
double  N 
)
protected

Definition at line 378 of file OpFastScintillation.cxx.

References larg4::OpDetPhotonTable::AddOpDetBacktrackerRecord(), larg4::OpDetPhotonTable::AddPhoton(), sim::OpDetBacktrackerRecord::AddScintillationPhotons(), energy, sim::OnePhoton::Energy, ExcitationRatio, fFiniteRiseTime, sim::LArG4Parameters::FillSimEnergyDeposits(), phot::PhotonVisibilityService::GetAllVisibilities(), geo::OpDetGeo::GetCenter(), larg4::ParticleListAction::GetCurrentTrackID(), phot::PhotonVisibilityService::GetReflT0s(), phot::PhotonVisibilityService::GetTimingTF1(), GetVisibleTimeOnlyCathode(), GetVUVTime(), phot::PhotonVisibilityService::IncludeParPropTime(), phot::PhotonVisibilityService::IncludePropTime(), sim::OnePhoton::InitialPosition, larg4::IonizationAndScintillation::Instance(), larg4::OpDetPhotonTable::Instance(), min, sim::OnePhoton::MotherTrackID, phot::PhotonVisibilityService::NOpChannels(), sim::LArG4Parameters::NoPhotonPropagation(), geo::GeometryCore::OpDetGeoFromOpDet(), ProcessStep(), sample_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().

379 {
380 
381  // Get the pointer to the fast scintillation table
382  OpDetPhotonTable * fst = OpDetPhotonTable::Instance();
383  OpDetPhotonTable* litefst = OpDetPhotonTable::Instance();
384 
385  // Get the pointer to the visibility service
386  art::ServiceHandle<phot::PhotonVisibilityService> pvs;
387  art::ServiceHandle<sim::LArG4Parameters> lgp;
388 
389  const G4Track * aTrack = aStep.GetTrack();
390 
391  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
392  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
393 
394  const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
395  const G4Material* aMaterial = aTrack->GetMaterial();
396 
397  G4int materialIndex = aMaterial->GetIndex();
398  G4int tracknumber = aStep.GetTrack()->GetTrackID();
399 
400  G4ThreeVector x0 = pPreStepPoint->GetPosition();
401  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
402  //G4double t0 = pPreStepPoint->GetGlobalTime() - fGlobalTimeOffset;
403  G4double t0 = pPreStepPoint->GetGlobalTime();
404 
405 
406  G4MaterialPropertiesTable* aMaterialPropertiesTable =
407  aMaterial->GetMaterialPropertiesTable();
408 
409  // Get the visibility vector for this point
410  size_t NOpChannels = 0;
411  NOpChannels = pvs->NOpChannels();
412 
413 
414  G4MaterialPropertyVector* Fast_Intensity =
415  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
416  G4MaterialPropertyVector* Slow_Intensity =
417  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
418 
419  if (!Fast_Intensity && !Slow_Intensity )
420  return 1;
421 
422  if(lgp->FillSimEnergyDeposits())
423  ProcessStep(aStep);
424 
425  if(lgp->NoPhotonPropagation())
426  return 0;
427 
428  G4int nscnt = 1;
429  if (Fast_Intensity && Slow_Intensity) nscnt = 2;
430 
431 
432  double Num = 0;
433  double YieldRatio=0;
434 
435 
436  if (scintillationByParticleType) {
437  // The scintillation response is a function of the energy
438  // deposited by particle types.
439 
440  // Get the definition of the current particle
441  G4ParticleDefinition *pDef = aParticle->GetDefinition();
442 
443  // Obtain the G4MaterialPropertyVectory containing the
444  // scintillation light yield as a function of the deposited
445  // energy for the current particle type
446 
447  // Protons
448  if(pDef==G4Proton::ProtonDefinition())
449  {
450  YieldRatio = aMaterialPropertiesTable->
451  GetConstProperty("PROTONYIELDRATIO");
452 
453  }
454 
455  // Muons
456  else if(pDef==G4MuonPlus::MuonPlusDefinition()||pDef==G4MuonMinus::MuonMinusDefinition())
457  {
458  YieldRatio = aMaterialPropertiesTable->
459  GetConstProperty("MUONYIELDRATIO");
460  }
461 
462  // Pions
463  else if(pDef==G4PionPlus::PionPlusDefinition()||pDef==G4PionMinus::PionMinusDefinition())
464  {
465  YieldRatio = aMaterialPropertiesTable->
466  GetConstProperty("PIONYIELDRATIO");
467  }
468 
469  // Kaons
470  else if(pDef==G4KaonPlus::KaonPlusDefinition()||pDef==G4KaonMinus::KaonMinusDefinition())
471  {
472  YieldRatio = aMaterialPropertiesTable->
473  GetConstProperty("KAONYIELDRATIO");
474  }
475 
476  // Alphas
477  else if(pDef==G4Alpha::AlphaDefinition())
478  {
479  YieldRatio = aMaterialPropertiesTable->
480  GetConstProperty("ALPHAYIELDRATIO");
481  }
482 
483  // Electrons (must also account for shell-binding energy
484  // attributed to gamma from standard PhotoElectricEffect)
485  else if(pDef==G4Electron::ElectronDefinition() ||
486  pDef==G4Gamma::GammaDefinition())
487  {
488  YieldRatio = aMaterialPropertiesTable->
489  GetConstProperty("ELECTRONYIELDRATIO");
490  }
491 
492  // Default for particles not enumerated/listed above
493  else
494  {
495  YieldRatio = aMaterialPropertiesTable->
496  GetConstProperty("ELECTRONYIELDRATIO");
497  }
498 
499  // If the user has not specified yields for (p,d,t,a,carbon)
500  // then these unspecified particles will default to the
501  // electron's scintillation yield
502  if(YieldRatio==0){
503  {
504 
505  YieldRatio = aMaterialPropertiesTable->
506  GetConstProperty("ELECTRONYIELDRATIO");
507 
508  }
509  }
510  }
511 
512  double const xyz[3] = { x0[0]/CLHEP::cm, x0[1]/CLHEP::cm, x0[2]/CLHEP::cm };
513  float const* Visibilities = pvs->GetAllVisibilities(xyz);
514 
515  float const* ReflVisibilities = nullptr;
516  float const* ReflT0s = nullptr;
517 
518  TF1* ParPropTimeTF1 = nullptr;
519 
520  if(pvs->StoreReflected()) {
521  ReflVisibilities = pvs->GetAllVisibilities(xyz,true);
522  if(pvs->StoreReflT0())
523  ReflT0s = pvs->GetReflT0s(xyz);
524  }
525  if(pvs->IncludeParPropTime())
526  {
527  ParPropTimeTF1 = pvs->GetTimingTF1(xyz);
528  }
529  /*
530  // For Kazu to debug # photons generated using csv file, by default should be commented out
531  // but don't remove as it's useful. Separated portion of code relevant to this block
532  // is labeled as "CASE-DEBUG DO NOT REMOVE THIS COMMENT"
533  // (2017 May 2nd ... "kazuhiro at nevis dot columbia dot edu")
534  //
535  static bool first_time=false;
536  if(!first_time) {
537  std::cout<<"LOGMEid,pdg,mass,ke,te,de,x,y,z,t,dr,mean_npe,gen_npe,det_npe"<<std::endl;
538  first_time=true;
539  }
540 
541  std::cout<<"LOGME"
542  <<aStep.GetTrack()->GetTrackID()<<","
543  <<aParticle->GetDefinition()->GetPDGEncoding()<<","
544  <<aParticle->GetDefinition()->GetPDGMass()/CLHEP::MeV<<","
545  <<pPreStepPoint->GetKineticEnergy()<<","
546  <<pPreStepPoint->GetTotalEnergy()<<","
547  <<aStep.GetTotalEnergyDeposit()<<","
548  <<xyz[0]<<","<<xyz[1]<<","<<xyz[2]<<","<<t0<<","
549  <<aStep.GetDeltaPosition().mag()<<","
550  <<MeanNumberOfPhotons<<","<<std::flush;
551 
552  double gen_photon_ctr=0;
553  double det_photon_ctr=0;
554  */
555  for (G4int scnt = 1; scnt <= nscnt; scnt++) {
556 
557  G4double ScintillationTime = 0.*CLHEP::ns;
558  G4double ScintillationRiseTime = 0.*CLHEP::ns;
559  G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
560 
561  if (scnt == 1) {
562  if (nscnt == 1) {
563  if(Fast_Intensity){
564  ScintillationTime = aMaterialPropertiesTable->
565  GetConstProperty("FASTTIMECONSTANT");
566  if (fFiniteRiseTime) {
567  ScintillationRiseTime = aMaterialPropertiesTable->
568  GetConstProperty("FASTSCINTILLATIONRISETIME");
569  }
570  ScintillationIntegral =
571  (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
572  }
573  if(Slow_Intensity){
574  ScintillationTime = aMaterialPropertiesTable->
575  GetConstProperty("SLOWTIMECONSTANT");
576  if (fFiniteRiseTime) {
577  ScintillationRiseTime = aMaterialPropertiesTable->
578  GetConstProperty("SLOWSCINTILLATIONRISETIME");
579  }
580  ScintillationIntegral =
581  (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
582  }
583  }//endif nscnt=1
584  else {
585  if(YieldRatio==0)
586  YieldRatio = aMaterialPropertiesTable->
587  GetConstProperty("YIELDRATIO");
588 
589 
590  if ( ExcitationRatio == 1.0 ) {
591  Num = std::min(YieldRatio,1.0)*MeanNumberOfPhotons;
592  }
593  else {
594  Num = std::min(ExcitationRatio,1.0)*MeanNumberOfPhotons;
595  }
596  ScintillationTime = aMaterialPropertiesTable->
597  GetConstProperty("FASTTIMECONSTANT");
598  if (fFiniteRiseTime) {
599  ScintillationRiseTime = aMaterialPropertiesTable->
600  GetConstProperty("FASTSCINTILLATIONRISETIME");
601  }
602  ScintillationIntegral =
603  (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
604  }//endif nscnt!=1
605  }//end scnt=1
606 
607  else {
608  Num = MeanNumberOfPhotons - Num;
609  ScintillationTime = aMaterialPropertiesTable->
610  GetConstProperty("SLOWTIMECONSTANT");
611  if (fFiniteRiseTime) {
612  ScintillationRiseTime = aMaterialPropertiesTable->
613  GetConstProperty("SLOWSCINTILLATIONRISETIME");
614  }
615  ScintillationIntegral =
616  (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
617  }
618 
619  if (!ScintillationIntegral) continue;
620 
621  //gen_photon_ctr += Num; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
622 
623  // Max Scintillation Integral
624 
625  // G4double CIImax = ScintillationIntegral->GetMaxValue();
626 
627 
628  //std::cout << "++++++++++++" << Num << "++++++++++" << std::endl;
629 
630 
631  // here we go: now if visibilities are invalid, we are in trouble
632  //if (!Visibilities && (NOpChannels > 0)) {
633  // throw cet::exception("OpFastScintillator")
634  // << "Photon library does not cover point ( " << xyz[0] << ", "
635  // << xyz[1] << ", " << xyz[2] << " ) cm.\n";
636  //}
637 
638  if(!Visibilities){
639  }else{
640  std::map<int, int> DetectedNum;
641 
642  std::map<int, int> ReflDetectedNum;
643 
644  for(size_t OpDet=0; OpDet!=NOpChannels; OpDet++)
645  {
646  G4int DetThisPMT = G4int(G4Poisson(Visibilities[OpDet] * Num));
647 
648  if(DetThisPMT>0)
649  {
650  DetectedNum[OpDet]=DetThisPMT;
651  // mf::LogInfo("OpFastScintillation") << "FastScint: " <<
652  // // it->second<<" " << Num << " " << DetThisPMT;
653 
654  //det_photon_ctr += DetThisPMT; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
655  }
656  if(pvs->StoreReflected()) {
657  G4int ReflDetThisPMT = G4int(G4Poisson(ReflVisibilities[OpDet] * Num));
658  if(ReflDetThisPMT>0)
659  {
660  ReflDetectedNum[OpDet]=ReflDetThisPMT;
661  }
662  }
663 
664  }
665  // Now we run through each PMT figuring out num of detected photons
666 
667  if(lgp->UseLitePhotons())
668  {
669  std::map<int, std::map<int, int>> StepPhotonTable;
670  // And then add these to the total collection for the event
671  for(std::map<int,int>::const_iterator itdetphot = DetectedNum.begin();
672  itdetphot!=DetectedNum.end(); ++itdetphot)
673  {
674  std::map<int, int> StepPhotons;
675  for (G4int i = 0; i < itdetphot->second; ++i)
676  {
677  G4double deltaTime = aStep.GetStepLength() /
678  ((pPreStepPoint->GetVelocity()+ pPostStepPoint->GetVelocity())/2.);
679 
680 
681  if (ScintillationRiseTime==0.0) {
682  deltaTime = deltaTime -
683  ScintillationTime * std::log( G4UniformRand() );
684  } else {
685  deltaTime = deltaTime +
686  sample_time(ScintillationRiseTime, ScintillationTime);
687  }
688 
689  if(pvs->IncludeParPropTime()) {
690  deltaTime += ParPropTimeTF1[itdetphot->first].GetRandom();
691  }
692 
693  G4double aSecondaryTime = t0 + deltaTime;
694  float Time = aSecondaryTime;
695  int ticks = static_cast<int>(Time);
696  StepPhotons[ticks]++;
697  }
698  StepPhotonTable[itdetphot->first] = StepPhotons;
699  //Iterate over Step Photon Table to add photons to OpDetBacktrackerRecords.
700 
701  sim::OpDetBacktrackerRecord tmpOpDetBTRecord(itdetphot->first);
702  //int thisG4TrackID = (aStep.GetTrack())->GetTrackID();
703  int thisG4TrackID = ParticleListAction::GetCurrentTrackID();
704  CLHEP::Hep3Vector prePoint = (aStep.GetPreStepPoint())->GetPosition();
705  CLHEP::Hep3Vector postPoint = (aStep.GetPostStepPoint())->GetPosition();
706  //Note the use of xO (letter O) instead of x0. This is to differentiate the positions here with the earlier declared double* x0
707  double xO = ( ( (prePoint.getX() + postPoint.getX() ) / 2.0) / CLHEP::cm );
708  double yO = ( ( (prePoint.getY() + postPoint.getY() ) / 2.0) / CLHEP::cm );
709  double zO = ( ( (prePoint.getZ() + postPoint.getZ() ) / 2.0) / CLHEP::cm );
710  double const xyzPos[3] = {xO,yO,zO};
711 // double energy = ( aStep.GetTotalEnergyDeposit() / CLHEP::MeV );
712  double energy = larg4::IonizationAndScintillation::Instance()->VisibleEnergyDeposit()/CLHEP::MeV;
713  //Loop over StepPhotons to get number of photons detected at each time for this channel and G4Step.
714  for(std::map<int,int>::iterator stepPhotonsIt = StepPhotons.begin(); stepPhotonsIt != StepPhotons.end(); ++stepPhotonsIt)
715  {
716  int photonTime = stepPhotonsIt->first;
717  int numPhotons = stepPhotonsIt->second;
718  tmpOpDetBTRecord.AddScintillationPhotons(thisG4TrackID, photonTime, numPhotons, xyzPos, energy);
719  }
720  //Add OpDetBackTrackerRecord. (opdetphotonTABLE->instance().addOpDetBacktrackerRecord(sim::OpDetBacktrackerRecord BTRrecord)
721  litefst->AddOpDetBacktrackerRecord(tmpOpDetBTRecord);
722  }
723  litefst->AddPhoton(&StepPhotonTable);
724  }
725  else
726  {
727  // We need to know the positions of the different optical detectors
728  art::ServiceHandle<geo::Geometry> geo;
729  // And then add these to the total collection for the event
730  for(std::map<int,int>::const_iterator itdetphot = DetectedNum.begin();
731  itdetphot!=DetectedNum.end(); ++itdetphot)
732  {
733  G4double propagation_time = 0;
734  std::vector<double> arrival_time_dist_vuv;
735  if(pvs->IncludePropTime()) {
736  // Get VUV photons arrival time distribution from the parametrization
737  double OpDetCenter[3];
738  geo->OpDetGeoFromOpDet(itdetphot->first).GetCenter(OpDetCenter);
739  G4ThreeVector OpDetPoint(OpDetCenter[0]*CLHEP::cm,OpDetCenter[1]*CLHEP::cm,OpDetCenter[2]*CLHEP::cm);
740  double distance_in_cm = (x0 - OpDetPoint).mag()/CLHEP::cm; // this must be in CENTIMETERS!
741  arrival_time_dist_vuv = GetVUVTime(distance_in_cm, itdetphot->second);//in ns
742  if(!(arrival_time_dist_vuv.size()>0)) {
743  //G4cout<<"Number of VUV detected photons = "<<itdetphot->second<<" ... too small, => No photons simulated!"<<G4endl;
744  continue;
745  }
746  }
747 
748  for (G4int i = 0; i < itdetphot->second; ++i)
749  {
750  G4double deltaTime = aStep.GetStepLength() /
751  ((pPreStepPoint->GetVelocity()+
752  pPostStepPoint->GetVelocity())/2.);
753 
754  if (ScintillationRiseTime==0.0)
755  {
756  deltaTime = deltaTime -
757  ScintillationTime * std::log( G4UniformRand() );
758  }
759  else
760  {
761  deltaTime = deltaTime +
762  sample_time(ScintillationRiseTime, ScintillationTime);
763  }
764 
765  G4double aSecondaryTime = t0 + deltaTime;
766  if(pvs->IncludePropTime())
767  propagation_time = arrival_time_dist_vuv.at(i)*CLHEP::ns;
768 
769  // The sim photon in this case stores its production point and time
770  TVector3 PhotonPosition(x0[0],x0[1],x0[2]);
771 
772  // We don't know anything about the momentum dir, so set it to be Z
773  float Energy = 9.7*CLHEP::eV;
774  float Time = aSecondaryTime;
775  if(pvs->IncludePropTime())
776  Time += propagation_time;
777 
778  // Make a photon object for the collection
779  sim::OnePhoton PhotToAdd;
780  PhotToAdd.InitialPosition = PhotonPosition;
781  PhotToAdd.Energy = Energy;
782  PhotToAdd.Time = Time;
783  PhotToAdd.SetInSD = false;
784  PhotToAdd.MotherTrackID = tracknumber;
785 
786  fst->AddPhoton(itdetphot->first, std::move(PhotToAdd));
787  }
788  }
789  // repeat the same for the reflected light, in case it has been sored
790  if(pvs->StoreReflected())
791  {
792  for(std::map<int,int>::const_iterator itdetphot = ReflDetectedNum.begin();
793  itdetphot!=ReflDetectedNum.end(); ++itdetphot)
794  {
795  G4double propagation_time = 0;
796  std::vector<double> arrival_time_dist_vis;
797  if(pvs->IncludePropTime()) {
798  // Get Visible photons arrival time distribution from the parametrization
799  double OpDetCenter[3];
800  geo->OpDetGeoFromOpDet(itdetphot->first).GetCenter(OpDetCenter);
801  G4ThreeVector OpDetPoint(OpDetCenter[0]*CLHEP::cm,OpDetCenter[1]*CLHEP::cm,OpDetCenter[2]*CLHEP::cm);
802  // currently unused:
803  // double distance_in_cm = (x0 - OpDetPoint).mag()/CLHEP::cm; // this must be in CENTIMETERS!
804  double t0_vis_lib = ReflT0s[itdetphot->first];
805  arrival_time_dist_vis = GetVisibleTimeOnlyCathode(t0_vis_lib, itdetphot->second);
806  if(!(arrival_time_dist_vis.size()>0)) {
807  //G4cout<<"Number of Visible detected photons = "<<itdetphot->second<<" ... too small, => No photons simulated!"<<G4endl;
808  continue;
809  }
810  }
811  for (G4int i = 0; i < itdetphot->second; ++i)
812  {
813  G4double deltaTime = aStep.GetStepLength() /
814  ((pPreStepPoint->GetVelocity()+
815  pPostStepPoint->GetVelocity())/2.);
816 
817  if (ScintillationRiseTime==0.0)
818  {
819  deltaTime = deltaTime -
820  ScintillationTime * std::log( G4UniformRand() );
821  }
822  else
823  {
824  deltaTime = deltaTime +
825  sample_time(ScintillationRiseTime, ScintillationTime);
826  }
827  G4double aSecondaryTime = t0 + deltaTime;
828  if(pvs->IncludePropTime())
829  propagation_time = arrival_time_dist_vis.at(i)*CLHEP::ns;
830 
831  TVector3 PhotonPosition(x0[0],x0[1],x0[2]);
832 
833  // We don't know anything about the momentum dir, so set it to be Z
834  std::map<double,double> tpbemission=lar::providerFrom<detinfo::LArPropertiesService>()->TpbEm();
835  const int nbins=tpbemission.size();
836 
837  double * parent=new double[nbins];
838 
839  int ii=0;
840  for( std::map<double, double>::iterator iter = tpbemission.begin(); iter != tpbemission.end(); ++iter)
841  { parent[ii++]=(*iter).second;
842  }
843 
844  CLHEP::RandGeneral rgen0(parent,nbins);
845 
846  double energy_reemited = rgen0.fire()*((*(--tpbemission.end())).first-(*tpbemission.begin()).first)+(*tpbemission.begin()).first;
847 
848  float Energy = energy_reemited*CLHEP::eV;
849  float Time = aSecondaryTime;
850  if(pvs->IncludePropTime())
851  Time += propagation_time;
852 
853  // Make a photon object for the collection
854  sim::OnePhoton PhotToAdd;
855  PhotToAdd.InitialPosition = PhotonPosition;
856  PhotToAdd.Energy = Energy;
857  PhotToAdd.Time = Time;
858  PhotToAdd.MotherTrackID = tracknumber;
859  PhotToAdd.SetInSD = false;
860  fst->AddPhoton(itdetphot->first, std::move(PhotToAdd));
861  delete [] parent;
862  }
863  }
864  }
865  }
866  }
867  }
868 
869  //std::cout<<gen_photon_ctr<<","<<det_photon_ctr<<std::endl; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
870  return 0;
871 }
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
geo::GeometryCore::OpDetGeoFromOpDet
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Number of OpDets in the whole detector.
Definition: GeometryCore.cxx:1799
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
larg4::OpFastScintillation::sample_time
G4double sample_time(G4double tau1, G4double tau2)
Definition: OpFastScintillation.cxx:1069
larg4::OpFastScintillation::ProcessStep
void ProcessStep(const G4Step &step)
Definition: OpFastScintillation.cxx:353
geo::OpDetGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: OpDetGeo.cxx:47
sim::LArG4Parameters::NoPhotonPropagation
bool NoPhotonPropagation() const
Definition: LArG4Parameters.h:68
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:309
phot::PhotonVisibilityService::GetTimingTF1
TF1 * GetTimingTF1(double const *xyz) const
Definition: PhotonVisibilityService_service.cc:492
energy
double energy
Definition: plottest35.C:25
larg4::IonizationAndScintillation::Instance
static IonizationAndScintillation * Instance()
Definition: IonizationAndScintillation.cxx:38
larg4::OpFastScintillation::GetVisibleTimeOnlyCathode
std::vector< double > GetVisibleTimeOnlyCathode(double, int)
Definition: OpFastScintillation.cxx:1165
sim::LArG4Parameters::FillSimEnergyDeposits
bool FillSimEnergyDeposits() const
Definition: LArG4Parameters.h:66
larg4::OpDetPhotonTable::Instance
static OpDetPhotonTable * Instance(bool LitePhotons=false)
Definition: OpDetPhotonTable.cxx:36
sim::OnePhoton::InitialPosition
TVector3 InitialPosition
Definition: SimPhotons.h:49
min
Int_t min
Definition: plot.C:26
geo
Namespace collecting geometry-related classes utilities.
Definition: TrackContainmentAlg.hh:20
larg4::OpFastScintillation::GetVUVTime
std::vector< double > GetVUVTime(double, int)
Definition: OpFastScintillation.cxx:1095
sim::LArG4Parameters::UseLitePhotons
bool UseLitePhotons() const
Definition: LArG4Parameters.h:64
phot::PhotonVisibilityService::GetReflT0s
float const * GetReflT0s(double const *xyz) const
Definition: PhotonVisibilityService_service.cc:442
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 
)
private

Definition at line 1069 of file OpFastScintillation.cxx.

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

Referenced by RecordPhotonsProduced().

1070 {
1071 // tau1: rise time and tau2: decay time
1072 
1073  while(1) {
1074  // two random numbers
1075  G4double ran1 = G4UniformRand();
1076  G4double ran2 = G4UniformRand();
1077  //
1078  // exponential distribution as envelope function: very efficient
1079  //
1080  G4double d = (tau1+tau2)/tau2;
1081  // make sure the envelope function is
1082  // always larger than the bi-exponential
1083  G4double t = -1.0*tau2*std::log(1-ran1);
1084  G4double g = d*single_exp(t,tau2);
1085  if (ran2 <= bi_exp(t,tau1,tau2)/g) return t;
1086  }
1087  return -1.0;
1088 }
larg4::OpFastScintillation::single_exp
G4double single_exp(G4double t, G4double tau2)
Definition: OpFastScintillation.hh:389
d
Float_t d
Definition: plot.C:237
larg4::OpFastScintillation::bi_exp
G4double bi_exp(G4double t, G4double tau1, G4double tau2)
Definition: OpFastScintillation.hh:395
void larg4::OpFastScintillation::SetFiniteRiseTime ( const G4bool  state)
inline

Definition at line 309 of file OpFastScintillation.hh.

References fFiniteRiseTime.

310 {
311  fFiniteRiseTime = state;
312 }
void larg4::OpFastScintillation::SetScintillationByParticleType ( const G4bool  scintType)

Definition at line 1032 of file OpFastScintillation.cxx.

References emSaturation, RemoveSaturation(), and scintillationByParticleType.

Referenced by GetSaturation().

1033 {
1034  if (emSaturation) {
1035  G4Exception("OpFastScintillation::SetScintillationByParticleType", "Scint02",
1036  JustWarning, "Redefinition: Birks Saturation is replaced by ScintillationByParticleType!");
1037  RemoveSaturation();
1038  }
1039  scintillationByParticleType = scintType;
1040 }
void larg4::OpFastScintillation::SetScintillationExcitationRatio ( const G4double  excitationratio)
inline

Definition at line 339 of file OpFastScintillation.hh.

References ExcitationRatio.

340 {
341  ExcitationRatio = excitationratio;
342 }
void larg4::OpFastScintillation::SetScintillationYieldFactor ( const G4double  yieldfactor)
inline

Definition at line 327 of file OpFastScintillation.hh.

References YieldFactor.

328 {
329  YieldFactor = yieldfactor;
330 }
void larg4::OpFastScintillation::SetTrackSecondariesFirst ( const G4bool  state)
inline

Definition at line 303 of file OpFastScintillation.hh.

References fTrackSecondariesFirst.

304 {
305  fTrackSecondariesFirst = state;
306 }
G4double larg4::OpFastScintillation::single_exp ( G4double  t,
G4double  tau2 
)
inlineprivate

Definition at line 389 of file OpFastScintillation.hh.

Referenced by sample_time().

390 {
391  return std::exp(-1.0*t/tau2)/tau2;
392 }

Member Data Documentation

G4EmSaturation* larg4::OpFastScintillation::emSaturation
private

Definition at line 272 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 276 of file OpFastScintillation.hh.

Referenced by GetVUVTime(), and OpFastScintillation().

double larg4::OpFastScintillation::fd_max
private

Definition at line 277 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 279 of file OpFastScintillation.hh.

Referenced by GetVisibleTimeOnlyCathode(), and OpFastScintillation().

double larg4::OpFastScintillation::ft0_max
private

Definition at line 279 of file OpFastScintillation.hh.

Referenced by GetVisibleTimeOnlyCathode(), and OpFastScintillation().

double larg4::OpFastScintillation::ftf1_sampling_factor
private

Definition at line 278 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 275 of file OpFastScintillation.hh.

Referenced by GetVisibleTimeOnlyCathode(), and OpFastScintillation().

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

Definition at line 274 of file OpFastScintillation.hh.

Referenced by GetVUVTime(), and OpFastScintillation().

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().

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: