00001 #ifndef LARPROPERTIES_CC
00002 #define LARPROPERTIES_CC
00003
00004 #include "LArProperties.hh"
00005
00006 namespace larutil {
00007
00008 LArProperties* LArProperties::_me = 0;
00009
00010 LArProperties::LArProperties(bool default_load) : LArUtilBase()
00011 {
00012 _name = "LArProperties";
00013 if(default_load){
00014 _file_name = Form("%s/LArUtil/dat/%s",
00015 getenv("LARLIGHT_CORE_DIR"),
00016 kUTIL_DATA_FILENAME[LArUtilConfig::Detector()].c_str());
00017 _tree_name = kTREENAME_LARPROPERTIES;
00018 LoadData();
00019 }
00020 }
00021
00022 void LArProperties::ClearData()
00023 {
00024 fEfield.clear();
00025 fTemperature = larlight::DATA::INVALID_DOUBLE;
00026 fElectronlifetime = larlight::DATA::INVALID_DOUBLE;
00027 fRadiationLength = larlight::DATA::INVALID_DOUBLE;
00028 fArgon39DecayRate = larlight::DATA::INVALID_DOUBLE;
00029 fZ = larlight::DATA::INVALID_DOUBLE;
00030 fA = larlight::DATA::INVALID_DOUBLE;
00031 fI = larlight::DATA::INVALID_DOUBLE;
00032 fSa = larlight::DATA::INVALID_DOUBLE;
00033 fSx0 = larlight::DATA::INVALID_DOUBLE;
00034 fSx1 = larlight::DATA::INVALID_DOUBLE;
00035 fScbar = larlight::DATA::INVALID_DOUBLE;
00036 fFastScintSpectrum.clear();
00037 fFastScintEnergies.clear();
00038 fSlowScintSpectrum.clear();
00039 fSlowScintEnergies.clear();
00040 fRIndexSpectrum.clear();
00041 fRIndexEnergies.clear();
00042 fAbsLengthSpectrum.clear();
00043 fAbsLengthEnergies.clear();
00044 fRayleighSpectrum.clear();
00045 fRayleighEnergies.clear();
00046 fScintByParticleType=false;
00047 fProtonScintYield=larlight::DATA::INVALID_DOUBLE;
00048 fProtonScintYieldRatio=larlight::DATA::INVALID_DOUBLE;
00049 fPionScintYield=larlight::DATA::INVALID_DOUBLE;
00050 fPionScintYieldRatio=larlight::DATA::INVALID_DOUBLE;
00051 fMuonScintYield=larlight::DATA::INVALID_DOUBLE;
00052 fMuonScintYieldRatio=larlight::DATA::INVALID_DOUBLE;
00053 fKaonScintYield=larlight::DATA::INVALID_DOUBLE;
00054 fKaonScintYieldRatio=larlight::DATA::INVALID_DOUBLE;
00055 fElectronScintYield=larlight::DATA::INVALID_DOUBLE;
00056 fElectronScintYieldRatio=larlight::DATA::INVALID_DOUBLE;
00057 fAlphaScintYield=larlight::DATA::INVALID_DOUBLE;
00058 fAlphaScintYieldRatio=larlight::DATA::INVALID_DOUBLE;
00059 fScintYield = larlight::DATA::INVALID_DOUBLE;
00060 fScintResolutionScale = larlight::DATA::INVALID_DOUBLE;
00061 fScintFastTimeConst = larlight::DATA::INVALID_DOUBLE;
00062 fScintSlowTimeConst = larlight::DATA::INVALID_DOUBLE;
00063 fScintYieldRatio = larlight::DATA::INVALID_DOUBLE;
00064 fScintBirksConstant = larlight::DATA::INVALID_DOUBLE;
00065 fEnableCerenkovLight = false;
00066 fReflectiveSurfaceNames.clear();
00067 fReflectiveSurfaceEnergies.clear();
00068 fReflectiveSurfaceReflectances.clear();
00069 fReflectiveSurfaceDiffuseFractions.clear();
00070 }
00071
00072 bool LArProperties::ReadTree()
00073 {
00074 ClearData();
00075
00076 TChain* ch = new TChain(_tree_name.c_str());
00077 ch->AddFile(_file_name.c_str());
00078
00079 std::string error_msg("");
00080 if(!(ch->GetBranch("fEfield"))) error_msg += " fEfield\n";
00081 if(!(ch->GetBranch("fTemperature"))) error_msg += " fTemperature\n";
00082 if(!(ch->GetBranch("fElectronlifetime"))) error_msg += " fElectronlifetime\n";
00083 if(!(ch->GetBranch("fRadiationLength"))) error_msg += " fRadiationLength\n";
00084 if(!(ch->GetBranch("fArgon39DecayRate"))) error_msg += " fArgon39DecayRate\n";
00085
00086 if(!(ch->GetBranch("fZ"))) error_msg += " fZ\n";
00087 if(!(ch->GetBranch("fA"))) error_msg += " fA\n";
00088 if(!(ch->GetBranch("fI"))) error_msg += " fI\n";
00089 if(!(ch->GetBranch("fSa"))) error_msg += " fSa\n";
00090 if(!(ch->GetBranch("fSx0"))) error_msg += " fSx0\n";
00091 if(!(ch->GetBranch("fSx1"))) error_msg += " fSx1\n";
00092 if(!(ch->GetBranch("fScbar"))) error_msg += " fScbar\n";
00093
00094 if(!(ch->GetBranch("fFastScintSpectrum"))) error_msg += " fFastScintSpectrum\n";
00095 if(!(ch->GetBranch("fFastScintEnergies"))) error_msg += " fFastScintEnergies\n";
00096 if(!(ch->GetBranch("fSlowScintSpectrum"))) error_msg += " fSlowScintSpectrum\n";
00097 if(!(ch->GetBranch("fSlowScintEnergies"))) error_msg += " fSlowScintEnergies\n";
00098 if(!(ch->GetBranch("fRIndexSpectrum"))) error_msg += " fRIndexSpectrum\n";
00099 if(!(ch->GetBranch("fRIndexEnergies"))) error_msg += " fRIndexEnergies\n";
00100 if(!(ch->GetBranch("fAbsLengthSpectrum"))) error_msg += " fAbsLengthSpectrum\n";
00101 if(!(ch->GetBranch("fAbsLengthEnergies"))) error_msg += " fAbsLengthEnergies\n";
00102 if(!(ch->GetBranch("fRayleighSpectrum"))) error_msg += " fRayleighSpectrum\n";
00103 if(!(ch->GetBranch("fRayleighEnergies"))) error_msg += " fRayleighEnergies\n";
00104
00105 if(!(ch->GetBranch("fScintByParticleType"))) error_msg += " fScintByParticleType\n";
00106
00107 if(!(ch->GetBranch("fProtonScintYield"))) error_msg += " fProtonScintYield\n";
00108 if(!(ch->GetBranch("fProtonScintYieldRatio"))) error_msg += " fProtonScintYieldRatio\n";
00109 if(!(ch->GetBranch("fMuonScintYield"))) error_msg += " fMuonScintYield\n";
00110 if(!(ch->GetBranch("fMuonScintYieldRatio"))) error_msg += " fMuonScintYieldRatio\n";
00111 if(!(ch->GetBranch("fPionScintYield"))) error_msg += " fPionScintYield\n";
00112 if(!(ch->GetBranch("fPionScintYieldRatio"))) error_msg += " fPionScintYieldRatio\n";
00113 if(!(ch->GetBranch("fKaonScintYield"))) error_msg += " fKaonScintYield\n";
00114 if(!(ch->GetBranch("fKaonScintYieldRatio"))) error_msg += " fKaonScintYieldRatio\n";
00115 if(!(ch->GetBranch("fElectronScintYield"))) error_msg += " fElectronScintYield\n";
00116 if(!(ch->GetBranch("fElectronScintYieldRatio"))) error_msg += " fElectronScintYieldRatio\n";
00117 if(!(ch->GetBranch("fAlphaScintYield"))) error_msg += " fAlphaScintYield\n";
00118 if(!(ch->GetBranch("fAlphaScintYieldRatio"))) error_msg += " fAlphaScintYieldRatio\n";
00119
00120 if(!(ch->GetBranch("fScintYield"))) error_msg += " fScintYield\n";
00121 if(!(ch->GetBranch("fScintResolutionScale"))) error_msg += " fScintResolutionScale\n";
00122 if(!(ch->GetBranch("fScintFastTimeConst"))) error_msg += " fScintFastTimeConst\n";
00123 if(!(ch->GetBranch("fScintSlowTimeConst"))) error_msg += " fScintSlowTimeConst\n";
00124 if(!(ch->GetBranch("fScintYieldRatio"))) error_msg += " fScintYieldRatio\n";
00125 if(!(ch->GetBranch("fScintBirksConstant"))) error_msg += " fScintBirksConstant\n";
00126
00127 if(!(ch->GetBranch("fEnableCerenkovLight"))) error_msg += " fEnableCerenkovLight\n";
00128
00129 if(!(ch->GetBranch("fReflectiveSurfaceNames")))
00130 error_msg += " fReflectiveSurfaceNames\n";
00131 if(!(ch->GetBranch("fReflectiveSurfaceEnergies")))
00132 error_msg += " fReflectiveSurfaceEnergies\n";
00133 if(!(ch->GetBranch("fReflectiveSurfaceReflectances")))
00134 error_msg += " fReflectiveSurfaceReflectances\n";
00135 if(!(ch->GetBranch("fReflectiveSurfaceDiffuseFractions")))
00136 error_msg += " fReflectiveSurfaceDiffuseFractions\n";
00137
00138 if(!error_msg.empty()) {
00139
00140 throw LArUtilException(Form("Missing following TBranches...\n%s",error_msg.c_str()));
00141
00142 return false;
00143 }
00144
00145 std::vector<Double_t> *pEfield=nullptr;
00146 ch->SetBranchAddress("fEfield",&pEfield);
00147
00148 ch->SetBranchAddress("fTemperature",&fTemperature);
00149 ch->SetBranchAddress("fElectronlifetime",&fElectronlifetime);
00150 ch->SetBranchAddress("fRadiationLength",&fRadiationLength);
00151 ch->SetBranchAddress("fArgon39DecayRate",&fArgon39DecayRate);
00152
00153 ch->SetBranchAddress("fZ",&fZ);
00154 ch->SetBranchAddress("fA",&fA);
00155 ch->SetBranchAddress("fI",&fI);
00156 ch->SetBranchAddress("fSa",&fSa);
00157 ch->SetBranchAddress("fSk",&fSk);
00158 ch->SetBranchAddress("fSx0",&fSx0);
00159 ch->SetBranchAddress("fSx1",&fSx1);
00160 ch->SetBranchAddress("fScbar",&fScbar);
00161
00162
00163 std::vector<Double_t> *pFastScintSpectrum=nullptr;
00164 std::vector<Double_t> *pFastScintEnergies=nullptr;
00165 std::vector<Double_t> *pSlowScintSpectrum=nullptr;
00166 std::vector<Double_t> *pSlowScintEnergies=nullptr;
00167 std::vector<Double_t> *pRIndexSpectrum=nullptr;
00168 std::vector<Double_t> *pRIndexEnergies=nullptr;
00169 std::vector<Double_t> *pAbsLengthSpectrum=nullptr;
00170 std::vector<Double_t> *pAbsLengthEnergies=nullptr;
00171 std::vector<Double_t> *pRayleighSpectrum=nullptr;
00172 std::vector<Double_t> *pRayleighEnergies=nullptr;
00173
00174 ch->SetBranchAddress("fFastScintSpectrum",&pFastScintSpectrum);
00175 ch->SetBranchAddress("fFastScintEnergies",&pFastScintEnergies);
00176 ch->SetBranchAddress("fSlowScintSpectrum",&pSlowScintSpectrum);
00177 ch->SetBranchAddress("fSlowScintEnergies",&pSlowScintEnergies);
00178 ch->SetBranchAddress("fRIndexSpectrum",&pRIndexSpectrum);
00179 ch->SetBranchAddress("fRIndexEnergies",&pRIndexEnergies);
00180 ch->SetBranchAddress("fAbsLengthSpectrum",&pAbsLengthSpectrum);
00181 ch->SetBranchAddress("fAbsLengthEnergies",&pAbsLengthEnergies);
00182 ch->SetBranchAddress("fRayleighSpectrum",&pRayleighSpectrum);
00183 ch->SetBranchAddress("fRayleighEnergies",&pRayleighEnergies);
00184
00185 ch->SetBranchAddress("fScintByParticleType",&fScintByParticleType);
00186
00187 ch->SetBranchAddress("fProtonScintYield", &fProtonScintYield);
00188 ch->SetBranchAddress("fProtonScintYieldRatio", &fProtonScintYieldRatio);
00189 ch->SetBranchAddress("fMuonScintYield", &fMuonScintYield);
00190 ch->SetBranchAddress("fMuonScintYieldRatio", &fMuonScintYieldRatio);
00191 ch->SetBranchAddress("fPionScintYield", &fPionScintYield);
00192 ch->SetBranchAddress("fPionScintYieldRatio", &fPionScintYieldRatio);
00193 ch->SetBranchAddress("fKaonScintYield", &fKaonScintYield);
00194 ch->SetBranchAddress("fKaonScintYieldRatio", &fKaonScintYieldRatio);
00195 ch->SetBranchAddress("fElectronScintYield", &fElectronScintYield);
00196 ch->SetBranchAddress("fElectronScintYieldRatio", &fElectronScintYieldRatio);
00197 ch->SetBranchAddress("fAlphaScintYield", &fAlphaScintYield);
00198 ch->SetBranchAddress("fAlphaScintYieldRatio", &fAlphaScintYieldRatio);
00199
00200 ch->SetBranchAddress("fScintYield", &fScintYield);
00201 ch->SetBranchAddress("fScintResolutionScale", &fScintResolutionScale);
00202 ch->SetBranchAddress("fScintFastTimeConst", &fScintFastTimeConst);
00203 ch->SetBranchAddress("fScintSlowTimeConst", &fScintSlowTimeConst);
00204 ch->SetBranchAddress("fScintYieldRatio", &fScintYieldRatio);
00205 ch->SetBranchAddress("fScintBirksConstant", &fScintBirksConstant);
00206
00207 ch->SetBranchAddress("fEnableCerenkovLight", &fEnableCerenkovLight);
00208
00209 std::vector<std::string> *pReflectiveSurfaceNames=nullptr;
00210 std::vector<Double_t> *pReflectiveSurfaceEnergies=nullptr;
00211 std::vector<std::vector<Double_t> > *pReflectiveSurfaceReflectances=nullptr;
00212 std::vector<std::vector<Double_t> > *pReflectiveSurfaceDiffuseFractions=nullptr;
00213
00214 ch->SetBranchAddress("fReflectiveSurfaceNames", &pReflectiveSurfaceNames);
00215 ch->SetBranchAddress("fReflectiveSurfaceEnergies", &pReflectiveSurfaceEnergies);
00216 ch->SetBranchAddress("fReflectiveSurfaceReflectances", &pReflectiveSurfaceReflectances);
00217 ch->SetBranchAddress("fReflectiveSurfaceDiffuseFractions", &pReflectiveSurfaceDiffuseFractions);
00218
00219 ch->GetEntry(0);
00220
00221
00222
00223 for(size_t i=0; i<pEfield->size(); ++i)
00224 fEfield.push_back(pEfield->at(i));
00225
00226 size_t n_entries = pFastScintSpectrum->size();
00227 fFastScintSpectrum.reserve(n_entries);
00228 fFastScintEnergies.reserve(n_entries);
00229 for(size_t i=0; i<n_entries; ++i) {
00230 fFastScintSpectrum.push_back(pFastScintSpectrum->at(i));
00231 fFastScintEnergies.push_back(pFastScintEnergies->at(i));
00232 }
00233 n_entries = pSlowScintSpectrum->size();
00234 fSlowScintSpectrum.reserve(n_entries);
00235 fSlowScintEnergies.reserve(n_entries);
00236 for(size_t i=0; i<n_entries; ++i) {
00237 fSlowScintSpectrum.push_back(pSlowScintSpectrum->at(i));
00238 fSlowScintEnergies.push_back(pSlowScintEnergies->at(i));
00239 }
00240 n_entries = pRIndexSpectrum->size();
00241 fRIndexSpectrum.reserve(n_entries);
00242 fRIndexEnergies.reserve(n_entries);
00243 for(size_t i=0; i<n_entries; ++i) {
00244 fRIndexSpectrum.push_back(pRIndexSpectrum->at(i));
00245 fRIndexEnergies.push_back(pRIndexEnergies->at(i));
00246 }
00247 n_entries = pAbsLengthSpectrum->size();
00248 fAbsLengthSpectrum.reserve(n_entries);
00249 fAbsLengthEnergies.reserve(n_entries);
00250 for(size_t i=0; i<n_entries; ++i) {
00251 fAbsLengthSpectrum.push_back(pAbsLengthSpectrum->at(i));
00252 fAbsLengthEnergies.push_back(pAbsLengthEnergies->at(i));
00253 }
00254 n_entries = pRayleighSpectrum->size();
00255 fRayleighSpectrum.reserve(n_entries);
00256 fRayleighEnergies.reserve(n_entries);
00257 for(size_t i=0; i<n_entries; ++i) {
00258 fRayleighSpectrum.push_back(pRayleighSpectrum->at(i));
00259 fRayleighEnergies.push_back(pRayleighEnergies->at(i));
00260 }
00261
00262
00263 size_t n_surface = pReflectiveSurfaceNames->size();
00264 fReflectiveSurfaceNames.reserve(n_surface);
00265 fReflectiveSurfaceEnergies.reserve(n_surface);
00266 fReflectiveSurfaceReflectances.reserve(n_surface);
00267 fReflectiveSurfaceDiffuseFractions.reserve(n_surface);
00268 for(size_t i=0; i<n_surface; ++i) {
00269
00270 fReflectiveSurfaceNames.push_back(pReflectiveSurfaceNames->at(i));
00271 fReflectiveSurfaceEnergies.push_back(pReflectiveSurfaceEnergies->at(i));
00272 fReflectiveSurfaceReflectances.push_back(pReflectiveSurfaceReflectances->at(i));
00273 fReflectiveSurfaceDiffuseFractions.push_back(pReflectiveSurfaceDiffuseFractions->at(i));
00274
00275 }
00276
00277 delete ch;
00278 return true;
00279 }
00280
00281 Double_t LArProperties::Density(Double_t temperature) const
00282 {
00283
00284 if(temperature == 0.)
00285 temperature = Temperature();
00286
00287 Double_t density = -0.00615*temperature + 1.928;
00288
00289 return density;
00290 }
00291
00292 Double_t LArProperties::Efield(UInt_t planegap) const
00293 {
00294 if(planegap >= fEfield.size())
00295 throw LArUtilException("requesting Electric field in a plane gap that is not defined");
00296
00297 return fEfield.at(planegap);
00298 }
00299
00300
00301 Double_t LArProperties::DriftVelocity(Double_t efield, Double_t temperature) const {
00302
00303
00304
00305
00306
00307
00308
00309
00310 if(efield == 0.)
00311 efield = Efield();
00312
00313 if(efield > 4.0) {
00314
00315 std::ostringstream msg;
00316 msg <<"DriftVelocity Warning! : E-field value of "
00317 << efield
00318 << " kV/cm is outside of range covered by drift"
00319 << " velocity parameterization. Returned value"
00320 << " may not be correct";
00321 print(larlight::MSG::WARNING,__FUNCTION__,msg.str());
00322 }
00323
00324
00325 if(temperature == 0.)
00326 temperature = Temperature();
00327
00328 if(temperature < 87.0 || temperature > 94.0) {
00329 std::ostringstream msg;
00330 msg << "DriftVelocity Warning! : Temperature value of "
00331 << temperature
00332 << " K is outside of range covered by drift velocity"
00333 << " parameterization. Returned value may not be"
00334 << " correct";
00335 print(larlight::MSG::WARNING,__FUNCTION__,msg.str());
00336 }
00337
00338 Double_t tshift = -87.203+temperature;
00339 Double_t xFit = 0.0938163-0.0052563*tshift-0.0001470*tshift*tshift;
00340 Double_t uFit = 5.18406+0.01448*tshift-0.003497*tshift*tshift-0.000516*tshift*tshift*tshift;
00341 Double_t vd;
00342
00343
00344
00345 Double_t P1 = -0.04640;
00346 Double_t P2 = 0.01712;
00347 Double_t P3 = 1.88125;
00348 Double_t P4 = 0.99408;
00349 Double_t P5 = 0.01172;
00350 Double_t P6 = 4.20214;
00351 Double_t T0 = 105.749;
00352
00353 Double_t P1W = -0.01481;
00354 Double_t P2W = -0.0075;
00355 Double_t P3W = 0.141;
00356 Double_t P4W = 12.4;
00357 Double_t P5W = 1.627;
00358 Double_t P6W = 0.317;
00359 Double_t T0W = 90.371;
00360
00361
00362
00363
00364 if (efield < xFit) vd=efield*uFit;
00365 else if (efield<0.619) {
00366 vd = ((P1*(temperature-T0)+1)
00367 *(P3*efield*std::log(1+P4/efield) + P5*std::pow(efield,P6))
00368 +P2*(temperature-T0));
00369 }
00370 else if (efield<0.699) {
00371 vd = 12.5*(efield-0.619)*((P1W*(temperature-T0W)+1)
00372 *(P3W*efield*std::log(1+P4W/efield) + P5W*std::pow(efield,P6W))
00373 +P2W*(temperature-T0W))+
00374 12.5*(0.699-efield)*((P1*(temperature-T0)+1)
00375 *(P3*efield*std::log(1+P4/efield) + P5*std::pow(efield,P6))
00376 +P2*(temperature-T0));
00377 }
00378 else {
00379 vd = ((P1W*(temperature-T0W)+1)
00380 *(P3W*efield*std::log(1+P4W/efield) + P5W*std::pow(efield,P6W))
00381 +P2W*(temperature-T0W));
00382 }
00383
00384 vd /= 10.;
00385
00386 return vd;
00387 }
00388
00389
00390
00391
00392
00393
00394
00395
00396 Double_t LArProperties::BirksCorrection(Double_t dQdx) const
00397 {
00398
00399
00400
00401 Double_t A3t = kRecombA;
00402 Double_t K3t = kRecombk;
00403 Double_t rho = this->Density();
00404 Double_t Wion = 1000./kGeVToElectrons;
00405 Double_t Efield = this->Efield();
00406 K3t /= rho;
00407 Double_t dEdx = dQdx/(A3t/Wion-K3t/Efield*dQdx);
00408
00409 return dEdx;
00410 }
00411
00412
00413 Double_t LArProperties::ModBoxCorrection(Double_t dQdx) const
00414 {
00415
00416
00417 Double_t rho = this->Density();
00418 Double_t Wion = 1000./kGeVToElectrons;
00419 Double_t Efield = this->Efield();
00420 Double_t Beta = kModBoxB / (rho * Efield);
00421 Double_t Alpha = kModBoxA;
00422 Double_t dEdx = (exp(Beta * Wion * dQdx ) - Alpha) / Beta;
00423
00424 return dEdx;
00425
00426 }
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 Double_t LArProperties::Eloss(Double_t mom, Double_t mass, Double_t tcut) const
00446 {
00447
00448
00449 Double_t K = 0.307075;
00450 Double_t me = 0.510998918;
00451
00452
00453
00454 Double_t bg = mom / mass;
00455 Double_t gamma = sqrt(1. + bg*bg);
00456 Double_t beta = bg / gamma;
00457 Double_t mer = 0.001 * me / mass;
00458 Double_t tmax = 2.*me* bg*bg / (1. + 2.*gamma*mer + mer*mer);
00459
00460
00461
00462 if(tcut == 0. || tcut > tmax)
00463 tcut = tmax;
00464
00465
00466
00467 Double_t x = std::log10(bg);
00468 Double_t delta = 0.;
00469 if(x >= fSx0) {
00470 delta = 2. * std::log(10.) * x - fScbar;
00471 if(x < fSx1)
00472 delta += fSa * std::pow(fSx1 - x, fSk);
00473 }
00474
00475
00476
00477 Double_t B = 0.5 * std::log(2.*me*bg*bg*tcut / (1.e-12 * fI*fI))
00478 - 0.5 * beta*beta * (1. + tcut / tmax) - 0.5 * delta;
00479
00480
00481
00482 if(B < 1.)
00483 B = 1.;
00484
00485
00486
00487 Double_t dedx = Density() * K*fZ*B / (fA * beta*beta);
00488
00489
00490
00491 return dedx;
00492 }
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503 Double_t LArProperties::ElossVar(Double_t mom, Double_t mass) const
00504 {
00505
00506
00507 Double_t K = 0.307075;
00508 Double_t me = 0.510998918;
00509
00510
00511
00512 Double_t bg = mom / mass;
00513 Double_t gamma2 = 1. + bg*bg;
00514 Double_t beta2 = bg*bg / gamma2;
00515
00516
00517
00518 Double_t result = gamma2 * (1. - 0.5 * beta2) * me * (fZ / fA) * K * Density();
00519 return result;
00520 }
00521
00522
00523 std::map<Double_t,Double_t> LArProperties::FastScintSpectrum() const
00524 {
00525 if(fFastScintSpectrum.size()!=fFastScintEnergies.size()){
00526 std::ostringstream msg;
00527 msg << "The vectors specifying the fast scintillation spectrum are "
00528 << " different sizes - " << fFastScintSpectrum.size()
00529 << " " << fFastScintEnergies.size();
00530 throw LArUtilException(msg.str());
00531 }
00532
00533 std::map<Double_t, Double_t> ToReturn;
00534 for(size_t i=0; i!=fFastScintSpectrum.size(); ++i)
00535 ToReturn[fFastScintEnergies.at(i)]=fFastScintSpectrum.at(i);
00536
00537 return ToReturn;
00538 }
00539
00540
00541 std::map<Double_t, Double_t> LArProperties::SlowScintSpectrum() const
00542 {
00543 if(fSlowScintSpectrum.size()!=fSlowScintEnergies.size()){
00544 std::ostringstream msg;
00545 msg << "The vectors specifying the slow scintillation spectrum are "
00546 << " different sizes - " << fFastScintSpectrum.size()
00547 << " " << fFastScintEnergies.size();
00548 throw LArUtilException(msg.str());
00549 }
00550
00551 std::map<Double_t, Double_t> ToReturn;
00552 for(size_t i=0; i!=fSlowScintSpectrum.size(); ++i)
00553 ToReturn[fSlowScintEnergies.at(i)]=fSlowScintSpectrum.at(i);
00554
00555 return ToReturn;
00556 }
00557
00558
00559 std::map<Double_t, Double_t> LArProperties::RIndexSpectrum() const
00560 {
00561 if(fRIndexSpectrum.size()!=fRIndexEnergies.size()){
00562 std::ostringstream msg;
00563 msg << "The vectors specifying the RIndex spectrum are "
00564 << " different sizes - " << fRIndexSpectrum.size()
00565 << " " << fRIndexEnergies.size();
00566 throw LArUtilException(msg.str());
00567 }
00568
00569 std::map<Double_t, Double_t> ToReturn;
00570 for(size_t i=0; i!=fRIndexSpectrum.size(); ++i)
00571 ToReturn[fRIndexEnergies.at(i)]=fRIndexSpectrum.at(i);
00572
00573 return ToReturn;
00574 }
00575
00576
00577
00578 std::map<Double_t, Double_t> LArProperties::AbsLengthSpectrum() const
00579 {
00580 if(fAbsLengthSpectrum.size()!=fAbsLengthEnergies.size()){
00581 std::ostringstream msg;
00582 msg << "The vectors specifying the Abs Length spectrum are "
00583 << " different sizes - " << fAbsLengthSpectrum.size()
00584 << " " << fAbsLengthEnergies.size();
00585 throw LArUtilException(msg.str());
00586 }
00587
00588 std::map<Double_t, Double_t> ToReturn;
00589 for(size_t i=0; i!=fAbsLengthSpectrum.size(); ++i)
00590 ToReturn[fAbsLengthEnergies.at(i)]=fAbsLengthSpectrum.at(i);
00591
00592 return ToReturn;
00593 }
00594
00595
00596 std::map<Double_t, Double_t> LArProperties::RayleighSpectrum() const
00597 {
00598 if(fRayleighSpectrum.size()!=fRayleighEnergies.size()){
00599 std::ostringstream msg;
00600 msg << "The vectors specifying the rayleigh spectrum are "
00601 << " different sizes - " << fRayleighSpectrum.size()
00602 << " " << fRayleighEnergies.size();
00603 throw LArUtilException(msg.str());
00604 }
00605
00606 std::map<Double_t, Double_t> ToReturn;
00607 for(size_t i=0; i!=fRayleighSpectrum.size(); ++i)
00608 ToReturn[fRayleighEnergies.at(i)]=fRayleighSpectrum.at(i);
00609
00610 return ToReturn;
00611 }
00612
00613
00614 std::map<std::string, std::map<Double_t,Double_t> > LArProperties::SurfaceReflectances() const
00615 {
00616 std::map<std::string, std::map<Double_t, Double_t> > ToReturn;
00617
00618 if(fReflectiveSurfaceNames.size()!=fReflectiveSurfaceReflectances.size()){
00619 std::ostringstream msg;
00620 msg << "The vectors specifying the surface reflectivities "
00621 << "do not have consistent sizes";
00622 throw LArUtilException(msg.str());
00623 }
00624 for(size_t i=0; i!=fReflectiveSurfaceNames.size(); ++i){
00625 if(fReflectiveSurfaceEnergies.size()!=fReflectiveSurfaceReflectances.at(i).size()){
00626 std::ostringstream msg;
00627 msg << "The vectors specifying the surface reflectivities do not have consistent sizes";
00628 throw LArUtilException(msg.str());
00629 }
00630 }
00631 for(size_t iName=0; iName!=fReflectiveSurfaceNames.size(); ++iName)
00632 for(size_t iEnergy=0; iEnergy!=fReflectiveSurfaceEnergies.size(); ++iEnergy)
00633 ToReturn[fReflectiveSurfaceNames.at(iName)][fReflectiveSurfaceEnergies.at(iEnergy)]=fReflectiveSurfaceReflectances.at(iName)[iEnergy];
00634
00635 return ToReturn;
00636
00637 }
00638
00639
00640 std::map<std::string, std::map<Double_t,Double_t> > LArProperties::SurfaceReflectanceDiffuseFractions() const
00641 {
00642 std::map<std::string, std::map<Double_t, Double_t> > ToReturn;
00643
00644 if(fReflectiveSurfaceNames.size()!=fReflectiveSurfaceDiffuseFractions.size()){
00645 std::ostringstream msg;
00646 msg << "The vectors specifying the surface reflectivities do not have consistent sizes";
00647 LArUtilException(msg.str());
00648 }
00649 for(size_t i=0; i!=fReflectiveSurfaceNames.size(); ++i){
00650 if(fReflectiveSurfaceEnergies.size()!=fReflectiveSurfaceDiffuseFractions.at(i).size()){
00651 std::ostringstream msg;
00652 msg << "The vectors specifying the surface reflectivities do not have consistent sizes";
00653 throw LArUtilException(msg.str());
00654 }
00655 }
00656 for(size_t iName=0; iName!=fReflectiveSurfaceNames.size(); ++iName)
00657 for(size_t iEnergy=0; iEnergy!=fReflectiveSurfaceEnergies.size(); ++iEnergy)
00658 ToReturn[fReflectiveSurfaceNames.at(iName)][fReflectiveSurfaceEnergies.at(iEnergy)]=fReflectiveSurfaceDiffuseFractions.at(iName)[iEnergy];
00659
00660 return ToReturn;
00661 }
00662
00663 }
00664 #endif