#include <HFShower.h>
Classes | |
struct | Hit |
Public Types | |
typedef std::pair< XYZPoint, double > | Spot |
typedef std::pair< unsigned int, double > | Step |
typedef Steps::const_iterator | step_iterator |
typedef std::vector< Step > | Steps |
typedef math::XYZVector | XYZPoint |
Public Member Functions | |
bool | compute () |
Compute the shower longitudinal and lateral development. | |
std::vector< Hit > | getHits (G4Step *aStep, bool forLibrary) |
std::vector< Hit > | getHits (G4Step *aStep) |
HFShower (std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p, int chk=0) | |
HFShower (const RandomEngine *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart) | |
void | initRun (G4ParticleTable *) |
virtual | ~HFShower () |
virtual | ~HFShower () |
Private Member Functions | |
double | gam (double x, double a) const |
std::vector< double > | getDDDArray (const std::string &, const DDsvalues_type &, int &) |
int | indexFinder (double x, const std::vector< double > &Fhist) |
void | makeSteps (int nsteps) |
double | transProb (double factor, double R, double r) |
Private Attributes | |
double | aloge |
double | alpEM |
double | alpHD |
bool | applyFidCut |
double | balanceEH |
double | betEM |
double | betHD |
HFCherenkov * | cherenkov |
int | chkFibre |
double | criticalEnergy |
double | depthStart |
double | depthStep |
std::vector< int > | detector |
double | e |
double | eSpotSize |
std::vector< double > | eStep |
HFFibre * | fibre |
std::vector< double > | gpar |
double | hcalDepthFactor |
int | infinity |
double | lambdaEM |
double | lambdaHD |
std::vector< double > | lamcurr |
std::vector< double > | lamdepth |
std::vector< double > | lamstep |
std::vector< double > | lamtotal |
int | lossesOpt |
double | maxTRfactor |
int | nDepthSteps |
std::vector< int > | nspots |
int | nTRsteps |
int | onEcal |
double | part |
double | probMax |
const RandomEngine * | random |
std::vector< double > | rlamStep |
double | tgamEM |
double | tgamHD |
const ECALProperties * | theECALproperties |
EcalHitMaker * | theGrid |
HcalHitMaker * | theHcalHitMaker |
const HCALProperties * | theHCALproperties |
HDShowerParametrization * | theParam |
double | theR1 |
double | theR2 |
double | theR3 |
double | transFactor |
double | transParam |
std::vector< double > | x0curr |
std::vector< double > | x0depth |
double | x0EM |
double | x0HD |
Definition at line 22 of file HFShower.h.
typedef std::pair<XYZPoint,double> HFShower::Spot |
Definition at line 29 of file HFShower.h.
typedef std::pair<unsigned int, double> HFShower::Step |
Definition at line 30 of file HFShower.h.
typedef Steps::const_iterator HFShower::step_iterator |
Definition at line 32 of file HFShower.h.
typedef std::vector<Step> HFShower::Steps |
Definition at line 31 of file HFShower.h.
typedef math::XYZVector HFShower::XYZPoint |
Definition at line 27 of file HFShower.h.
HFShower::HFShower | ( | const RandomEngine * | engine, |
HDShowerParametrization * | myParam, | ||
EcalHitMaker * | myGrid, | ||
HcalHitMaker * | myHcalHitMaker, | ||
int | onECAL, | ||
double | epart | ||
) |
Definition at line 30 of file HFShower.cc.
References aloge, HDShowerParametrization::alpe1(), HDShowerParametrization::alpe2(), alpEM, HDShowerParametrization::alph1(), HDShowerParametrization::alph2(), alpHD, balanceEH, HDShowerParametrization::bete1(), HDShowerParametrization::bete2(), betEM, HDShowerParametrization::beth1(), HDShowerParametrization::beth2(), betHD, criticalEnergy, debug, depthStart, depthStep, detector, e, HDShowerParametrization::e1(), HDShowerParametrization::e2(), EcalHitMaker::ecalHcalGapTotalL0(), EcalHitMaker::ecalHcalGapTotalX0(), HDShowerParametrization::ecalProperties(), EcalHitMaker::ecalTotalL0(), HDShowerParametrization::emax(), HDShowerParametrization::emid(), HDShowerParametrization::emin(), eSpotSize, RandomEngine::flatShoot(), HSParameters::getHDbalanceEH(), HSParameters::getHDcriticalEnergy(), HSParameters::getHDdepthStep(), HSParameters::getHDeSpotSize(), HSParameters::getHDhcalDepthFactor(), HSParameters::getHDlossesOpt(), HSParameters::getHDmaxTRfactor(), HSParameters::getHDnDepthSteps(), HSParameters::getHDnTRsteps(), HSParameters::getHDtransParam(), hcalDepthFactor, HDShowerParametrization::hcalProperties(), EcalHitMaker::hcalTotalL0(), HDShowerParametrization::hsParameters(), i, HCALProperties::interactionLength(), ECALProperties::interactionLength(), lambdaEM, lambdaHD, lamcurr, lamdepth, lamstep, lamtotal, create_public_lumi_plots::log, LogDebug, lossesOpt, makeSteps(), maxTRfactor, nDepthSteps, nTRsteps, onEcal, HDShowerParametrization::part1(), HDShowerParametrization::part2(), HDShowerParametrization::r1(), HDShowerParametrization::r2(), HDShowerParametrization::r3(), ECALProperties::radLenIncm(), HCALProperties::radLenIncm(), random, HDShowerParametrization::setCase(), tgamEM, tgamHD, theECALproperties, theGrid, theHCALproperties, theParam, theR1, theR2, theR3, transFactor, transParam, x0curr, x0depth, x0EM, and x0HD.
: theParam(myParam), theGrid(myGrid), theHcalHitMaker(myHcalHitMaker), onEcal(onECAL), e(epart), random(engine) { // To get an access to constants read in FASTCalorimeter // FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance(); // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc lossesOpt = myParam->hsParameters()->getHDlossesOpt(); nDepthSteps = myParam->hsParameters()->getHDnDepthSteps(); nTRsteps = myParam->hsParameters()->getHDnTRsteps(); transParam = myParam->hsParameters()->getHDtransParam(); eSpotSize = myParam->hsParameters()->getHDeSpotSize(); depthStep = myParam->hsParameters()->getHDdepthStep(); criticalEnergy = myParam->hsParameters()->getHDcriticalEnergy(); maxTRfactor = myParam->hsParameters()->getHDmaxTRfactor(); balanceEH = myParam->hsParameters()->getHDbalanceEH(); hcalDepthFactor = myParam->hsParameters()->getHDhcalDepthFactor(); // Special tr.size fluctuations transParam *= (1. + random->flatShoot()); // Special ad hoc long. extension + some fluctuations double depthExt; if (e < 50.) depthExt = 0.8 * (50. - e) / 50. + 0.3; else { if (e < 500.) depthExt = (500. - e) / 500. * 0.4 - 0.1; else depthExt = -0.1; } hcalDepthFactor += depthExt + 0.05 * (2.* random->flatShoot() - 1.); // normally 1, in HF - might be smaller to take into account // a narrowness of the HF shower (Cherenkov light) if( e < 50. ) transFactor = 0.5 - (50. - e ) / 50. * 0.2; else transFactor = 0.7 - (1000. - e ) /1000. * 0.2; // simple protection ... if(e < 0) e = 0.; // Get the Famos Histos pointer // myHistos = FamosHistos::instance(); // std::cout << " Hello FamosShower " << std::endl; theECALproperties = theParam->ecalProperties(); theHCALproperties = theParam->hcalProperties(); double emax = theParam->emax(); double emid = theParam->emid(); double emin = theParam->emin(); double effective = e; if( e < emid ) { theParam->setCase(1); // avoid "underflow" below Emin (for parameters calculation only) if(e < emin) effective = emin; } else theParam->setCase(2); // A bit coarse espot size for HF... eSpotSize *= 2.5; if(effective > 0.5 * emax) { eSpotSize *= 2.; if(effective > emax) { effective = emax; eSpotSize *= 3.; depthStep *= 2.; } } if(debug == 2 ) LogDebug("FastCalorimetry") << " HFShower : " << std::endl << " Energy " << e << std::endl << " lossesOpt " << lossesOpt << std::endl << " nDepthSteps " << nDepthSteps << std::endl << " nTRsteps " << nTRsteps << std::endl << " transParam " << transParam << std::endl << " eSpotSize " << eSpotSize << std::endl << " criticalEnergy " << criticalEnergy << std::endl << " maxTRfactor " << maxTRfactor << std::endl << " balanceEH " << balanceEH << std::endl << "hcalDepthFactor " << hcalDepthFactor << std::endl; double alpEM1 = theParam->alpe1(); double alpEM2 = theParam->alpe2(); double betEM1 = theParam->bete1(); double betEM2 = theParam->bete2(); double alpHD1 = theParam->alph1(); double alpHD2 = theParam->alph2(); double betHD1 = theParam->beth1(); double betHD2 = theParam->beth2(); double part1 = theParam->part1(); double part2 = theParam->part2(); aloge = std::log(effective); double edpar = (theParam->e1() + aloge * theParam->e2()) * effective; double aedep = std::log(edpar); if(debug == 2) LogDebug("FastCalorimetry") << " HFShower : " << std::endl << " edpar " << edpar << " aedep " << aedep << std::endl << " alpEM1 " << alpEM1 << std::endl << " alpEM2 " << alpEM2 << std::endl << " betEM1 " << betEM1 << std::endl << " betEM2 " << betEM2 << std::endl << " alpHD1 " << alpHD1 << std::endl << " alpHD2 " << alpHD2 << std::endl << " betHD1 " << betHD1 << std::endl << " betHD2 " << betHD2 << std::endl << " part1 " << part1 << std::endl << " part2 " << part2 << std::endl; // private members to set theR1 = theParam->r1(); theR2 = theParam->r2(); theR3 = theParam->r3(); alpEM = alpEM1 + alpEM2 * aedep; tgamEM = tgamma(alpEM); betEM = betEM1 - betEM2 * aedep; alpHD = alpHD1 + alpHD2 * aedep; tgamHD = tgamma(alpHD); betHD = betHD1 - betHD2 * aedep; part = part1 - part2 * aedep; if(part > 1.) part = 1.; // protection - just in case of if(debug == 2 ) LogDebug("FastCalorimetry") << " HFShower : " << std::endl << " alpEM " << alpEM << std::endl << " tgamEM " << tgamEM << std::endl << " betEM " << betEM << std::endl << " alpHD " << alpHD << std::endl << " tgamHD " << tgamHD << std::endl << " betHD " << betHD << std::endl << " part " << part << std::endl; if(onECAL){ lambdaEM = theParam->ecalProperties()->interactionLength(); x0EM = theParam->ecalProperties()->radLenIncm(); } else { lambdaEM = 0.; x0EM = 0.; } lambdaHD = theParam->hcalProperties()->interactionLength(); x0HD = theParam->hcalProperties()->radLenIncm(); if(debug == 2) LogDebug("FastCalorimetry") << " HFShower e " << e << std::endl << " x0EM = " << x0EM << std::endl << " x0HD = " << x0HD << std::endl << " lamEM = " << lambdaEM << std::endl << " lamHD = " << lambdaHD << std::endl; // Starting point of the shower // try first with ECAL lambda double sum1 = 0.; // lambda path from the ECAL/HF entrance; double sum2 = 0.; // lambda path from the interaction point; double sum3 = 0.; // x0 path from the interaction point; int nsteps = 0; // full number of longitudinal steps (counter); int nmoresteps; // how many longitudinal steps in addition to // one (if interaction happens there) in ECAL if(e < criticalEnergy ) nmoresteps = 1; else nmoresteps = nDepthSteps; double depthECAL = 0.; double depthGAP = 0.; double depthGAPx0 = 0.; if(onECAL ) { depthECAL = theGrid->ecalTotalL0(); // ECAL depth segment depthGAP = theGrid->ecalHcalGapTotalL0(); // GAP depth segment depthGAPx0 = theGrid->ecalHcalGapTotalX0(); // GAP depth x0 } double depthHCAL = theGrid->hcalTotalL0(); // HCAL depth segment double depthToHCAL = depthECAL + depthGAP; //--------------------------------------------------------------------------- // Depth simulation & various protections, among them // if too deep - get flat random in the allowed region // if no HCAL material behind - force to deposit in ECAL double maxDepth = depthToHCAL + depthHCAL - 1.1 * depthStep; double depthStart = std::log(1./random->flatShoot()); // starting point lambda unts if(e < emin) { if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : e <emin -> depthStart = 0" << std::endl; depthStart = 0.; } if(depthStart > maxDepth) { if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : depthStart too big ... = " << depthStart << std::endl; depthStart = maxDepth * random->flatShoot(); if( depthStart < 0.) depthStart = 0.; if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : depthStart re-calculated = " << depthStart << std::endl; } if( onECAL && e < emid ) { if((depthECAL - depthStart)/depthECAL > 0.2 && depthECAL > depthStep ) { depthStart = 0.5 * depthECAL * random->flatShoot(); if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : small energy, " << " depthStart reduced to = " << depthStart << std::endl; } } if( depthHCAL < depthStep) { if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : depthHCAL too small ... = " << depthHCAL << " depthStart -> forced to 0 !!!" << std::endl; depthStart = 0.; nmoresteps = 0; if(depthECAL < depthStep) { nsteps = -1; LogInfo("FastCalorimetry") << " FamosHFShower : too small ECAL and HCAL depths - " << " particle is lost !!! " << std::endl; } } if(debug) LogDebug("FastCalorimetry") << " FamosHFShower depths(lam) - " << std::endl << " ECAL = " << depthECAL << std::endl << " GAP = " << depthGAP << std::endl << " HCAL = " << depthHCAL << std::endl << " starting point = " << depthStart << std::endl; if( onEcal ) { if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : onECAL" << std::endl; if(depthStart < depthECAL) { if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : depthStart < depthECAL" << std::endl; if((depthECAL - depthStart)/depthECAL > 0.25 && depthECAL > depthStep) { if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : enough space to make ECAL step" << std::endl; // ECAL - one step nsteps++; sum1 += depthECAL; // at the end of step sum2 += depthECAL-depthStart; sum3 += sum2 * lambdaEM / x0EM; lamtotal.push_back(sum1); lamdepth.push_back(sum2); lamcurr.push_back(lambdaEM); lamstep.push_back(depthECAL-depthStart); x0depth.push_back(sum3); x0curr.push_back(x0EM); detector.push_back(1); if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : " << " in ECAL sum1, sum2 " << sum1 << " " << sum2 << std::endl; // // Gap - no additional step after ECAL // // just move further to HCAL over the gap sum1 += depthGAP; sum2 += depthGAP; sum3 += depthGAPx0; } // Just shift starting point to HCAL else { // cout << " FamosHFShower : not enough space to make ECAL step" << std::endl; if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl; depthStart = depthToHCAL; sum1 += depthStart; } } else { // GAP or HCAL if(depthStart >= depthECAL && depthStart < depthToHCAL ) { depthStart = depthToHCAL; // just a shift to HCAL for simplicity } sum1 += depthStart; if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl; } } else { // Forward if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : forward" << std::endl; sum1 += depthStart; } for (int i = 0; i < nmoresteps ; i++) { sum1 += depthStep; if (sum1 > (depthECAL + depthGAP + depthHCAL)) break; sum2 += depthStep; sum3 += sum2 * lambdaHD / x0HD; lamtotal.push_back(sum1); lamdepth.push_back(sum2); lamcurr.push_back(lambdaHD); lamstep.push_back(depthStep); x0depth.push_back(sum3); x0curr.push_back(x0HD); detector.push_back(3); nsteps++; } // Make fractions of energy and transverse radii at each step // PV // std::cout << "HFShower::HFShower() : Nsteps = " << nsteps << std::endl; if(nsteps > 0) { makeSteps(nsteps); } }
HFShower::~HFShower | ( | ) | [inline, virtual] |
Definition at line 41 of file HFShower.h.
{;}
HFShower::HFShower | ( | std::string & | name, |
const DDCompactView & | cpv, | ||
edm::ParameterSet const & | p, | ||
int | chk = 0 |
||
) |
Definition at line 25 of file HFShower.cc.
References DDFilteredView::addFilter(), applyFidCut, cherenkov, chkFibre, DDSpecificsFilter::equals, Exception, fibre, alcazmumu_cfi::filter, DDFilteredView::firstChild(), getDDDArray(), edm::ParameterSet::getParameter(), gpar, DDFilteredView::mergedSpecifics(), mergeVDriftHistosByStation::name, probMax, DDSpecificsFilter::setCriteria(), and relativeConstraints::value.
: cherenkov(0), fibre(0), chkFibre(chk) { edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower"); applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut"); probMax = m_HF.getParameter<double>("ProbMax"); edm::LogInfo("HFShower") << "HFShower:: Maximum probability cut off " << probMax << " Check flag " << chkFibre; G4String attribute = "ReadOutName"; G4String value = name; DDSpecificsFilter filter; DDValue ddv(attribute,value,0); filter.setCriteria(ddv,DDSpecificsFilter::equals); DDFilteredView fv(cpv); fv.addFilter(filter); bool dodet = fv.firstChild(); if ( dodet ) { DDsvalues_type sv(fv.mergedSpecifics()); //Special Geometry parameters int ngpar = 7; gpar = getDDDArray("gparHF", sv,ngpar); edm::LogInfo("HFShower") << "HFShower: " << ngpar << " gpar (cm)"; for (int ig=0; ig<ngpar; ig++) edm::LogInfo("HFShower") << "HFShower: gpar[" << ig << "] = " << gpar[ig]/cm << " cm"; } else { edm::LogError("HFShower") << "HFShower: cannot get filtered " << " view for " << attribute << " matching " << name; throw cms::Exception("Unknown", "HFShower") << "cannot match " << attribute << " to " << name <<"\n"; } cherenkov = new HFCherenkov(m_HF); fibre = new HFFibre(name, cpv, p); }
virtual HFShower::~HFShower | ( | ) | [virtual] |
bool HFShower::compute | ( | ) |
Compute the shower longitudinal and lateral development.
Definition at line 475 of file HFShower.cc.
References EcalHitMaker::addHit(), HcalHitMaker::addHit(), prof2calltree::count, debug, detector, patCandidatesForDimuonsSequences_cff::ecal, eStep, RandomEngine::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, i, getHLTprescales::index, indexFinder(), EcalCondDB::inf, infinity, j, lamstep, lamtotal, LogDebug, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, CosmicsPD_Skims::radius, random, query::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), ntuplemaker::status, theGrid, theHcalHitMaker, and transProb().
Referenced by CalorimetryManager::HDShowerSimulation().
{ // TimeMe theT("FamosHFShower::compute"); bool status = false; int numLongit = eStep.size(); if(debug) LogDebug("FastCalorimetry") << " FamosHFShower::compute - " << " N_long.steps required : " << numLongit << std::endl; if(numLongit > 0) { status = true; // Prepare the trsanverse probability function std::vector<double> Fhist; std::vector<double> rhist; for (int j = 0; j < nTRsteps + 1; j++) { rhist.push_back(maxTRfactor * j / nTRsteps ); Fhist.push_back(transProb(maxTRfactor,1.,rhist[j])); if(debug == 3) LogDebug("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl; } //================================================================ // Longitudinal steps //================================================================ for (int i = 0; i < numLongit ; i++) { double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i]; // vary the longitudinal profile if needed if(detector[i] != 1) currentDepthL0 *= hcalDepthFactor; if(debug) LogDebug("FastCalorimetry") << " FamosHFShower::compute - detector = " << detector[i] << " currentDepthL0 = " << currentDepthL0 << std::endl; double maxTRsize = maxTRfactor * rlamStep[i]; // in lambda units double rbinsize = maxTRsize / nTRsteps; double espot = eStep[i] / (double)nspots[i]; // re-adjust espot if( espot > 4. || espot < 0. ) LogDebug("FastCalorimetry") << " FamosHFShower::compute - unphysical espot = " << espot << std::endl; int ecal = 0; if(detector[i] != 1) { bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0); if(debug) LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of " << " theHcalHitMaker->setDepth(currentDepthL0) is " << setHDdepth << std::endl; if(!setHDdepth) { currentDepthL0 -= lamstep[i]; setHDdepth = theHcalHitMaker->setDepth(currentDepthL0); } if(!setHDdepth) continue; theHcalHitMaker->setSpotEnergy(espot); } else { ecal = 1; bool status = theGrid->getPads(currentDepthL0); if(debug) LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of Grid = " << status << std::endl; if(!status) continue; theGrid->setSpotEnergy(espot); } //------------------------------------------------------------ // Transverse distribution //------------------------------------------------------------ int nok = 0; // counter of OK int count = 0; int inf = infinity; if(lossesOpt) inf = nspots[i]; // losses are enabled, otherwise // only OK points are counted ... // Total energy in this layer double eremaining = eStep[i]; bool converged = false; while (eremaining > 0. && !converged && count<inf ) { ++count; // energy spot (HFL) double newespot = espot; // We need to know a priori if this energy spot if for // a long (1) or short (2) fiber unsigned layer = 1; if( currentDepthL0 < 1.3 ) // first 22 cm = 1.3 lambda - only HFL layer = 1; else layer = random->flatShoot() < 0.5 ? 1 : 2; if ( layer == 2 ) newespot = 2. * espot; if ( eremaining - newespot < 0. ) newespot = eremaining; // process transverse distribution double prob = random->flatShoot(); int index = indexFinder(prob,Fhist); double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize; // in-bin double phi = 2.*M_PI*random->flatShoot(); if(debug == 2) LogDebug("FastCalorimetry") << std::endl << " FamosHFShower::compute " << " r = " << radius << " phi = " << phi << std::endl; // add hit theHcalHitMaker->setSpotEnergy(newespot); theGrid->setSpotEnergy(newespot); bool result; if(ecal) { result = theGrid->addHit(radius,phi,0); // shouldn't get here ! if(debug == 2) LogDebug("FastCalorimetry") << " FamosHFShower::compute - " << " theGrid->addHit result = " << result << std::endl; } else { // PV assign espot to long/short fibers result = theHcalHitMaker->addHit(radius,phi,layer); if(debug == 2) LogDebug("FastCalorimetry") << " FamosHFShower::compute - " << " theHcalHitMaker->addHit result = " << result << std::endl; } if (result) { ++nok; eremaining -= newespot; if ( eremaining <= 0. ) converged = true; // std::cout << "Transverse : " << nok << " " // << " , E= " << newespot // << " , Erem = " << eremaining // << std::endl; } else { // std::cout << "WARNING : hit not added" << std::endl; } } // end of tranverse simulation //----------------------------------------------------- if(count == infinity) { status = false; if(debug) LogDebug("FastCalorimetry") << "*** FamosHFShower::compute " << " maximum number of" << " transverse points " << count << " is used !!!" << std::endl; break; } if(debug) LogDebug("FastCalorimetry") << " FamosHFShower::compute " << " long.step No." << i << " Ntry, Nok = " << count << " " << nok << std::endl; } // end of longitudinal steps } // end of no steps return status; }
double HFShower::gam | ( | double | x, |
double | a | ||
) | const [inline, private] |
Definition at line 49 of file HFShower.h.
References create_public_lumi_plots::exp, and funct::pow().
Referenced by makeSteps().
std::vector< double > HFShower::getDDDArray | ( | const std::string & | str, |
const DDsvalues_type & | sv, | ||
int & | nmin | ||
) | [private] |
Definition at line 301 of file HFShower.cc.
References DDfetch(), DDValue::doubles(), Exception, LogDebug, and relativeConstraints::value.
Referenced by HFShower().
{ #ifdef DebugLog LogDebug("HFShower") << "HFShower:getDDDArray called for " << str << " with nMin " << nmin; #endif DDValue value(str); if (DDfetch(&sv, value)) { #ifdef DebugLog LogDebug("HFShower") << value; #endif const std::vector<double> & fvec = value.doubles(); int nval = fvec.size(); if (nmin > 0) { if (nval < nmin) { edm::LogError("HFShower") << "HFShower : # of " << str << " bins " << nval << " < " << nmin << " ==> illegal"; throw cms::Exception("Unknown", "HFShower") << "nval < nmin for array " << str << "\n"; } } else { if (nval < 2) { edm::LogError("HFShower") << "HFShower : # of " << str << " bins " << nval << " < 2 ==> illegal" << " (nmin=" << nmin << ")"; throw cms::Exception("Unknown", "HFShower") << "nval < 2 for array " << str << "\n"; } } nmin = nval; return fvec; } else { edm::LogError("HFShower") << "HFShower : cannot get array " << str; throw cms::Exception("Unknown", "HFShower") << "cannot get array " << str << "\n"; } }
std::vector< HFShower::Hit > HFShower::getHits | ( | G4Step * | aStep, |
bool | forLibrary | ||
) |
Definition at line 202 of file HFShower.cc.
References abs, applyFidCut, cherenkov, chkFibre, HFCherenkov::computeNPE(), relval_parameters_module::energy, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, LogDebug, HFShower::Hit::momentum, mergeVDriftHistosByStation::name, convertSQLiteXML::ok, HFFibreFiducial::PMTNumber(), HFShower::Hit::position, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), v, w(), and HFShower::Hit::wavelength.
{ std::vector<HFShower::Hit> hits; int nHit = 0; double edep = aStep->GetTotalEnergyDeposit(); double stepl = 0.; if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.) stepl = aStep->GetStepLength(); if ((edep == 0.) || (stepl == 0.)) { #ifdef DebugLog LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit; #endif return hits; } G4Track *aTrack = aStep->GetTrack(); const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle(); HFShower::Hit hit; double energy = aParticle->GetTotalEnergy(); double momentum = aParticle->GetTotalMomentum(); double pBeta = momentum / energy; double dose = 0.; int npeDose = 0; G4ThreeVector momentumDir = aParticle->GetMomentumDirection(); G4ParticleDefinition *particleDef = aTrack->GetDefinition(); G4StepPoint * preStepPoint = aStep->GetPreStepPoint(); G4ThreeVector globalPos = preStepPoint->GetPosition(); G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName(); double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1]; G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(), zv); G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()-> GetTopTransform().TransformAxis(momentumDir); // @@ Here the depth should be changed (Fibers are all long in Geometry!) int depth = 1; int npmt = 0; bool ok = true; if (zv < 0 || zv > gpar[1]) { ok = false; // beyond the fiber in Z } if (ok && applyFidCut) { npmt = HFFibreFiducial:: PMTNumber(globalPos); if (npmt <= 0) { ok = false; } else if (npmt > 24) { // a short fibre if (zv > gpar[0]) { depth = 2; } else { ok = false; } } #ifdef DebugLog edm::LogInfo("HFShower") << "HFShower:getHits(SL): npmt " << npmt << " zv " << std::abs(pe_effect.z()) << ":" << gpar[4] << ":" << zv << ":" << gpar[0] << " ok " << ok << " depth " << depth; #endif } else { depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG! } G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation(); double u = localMom.x(); double v = localMom.y(); double w = localMom.z(); double zCoor = localPos.z(); double zFibre = (0.5*gpar[1]-zCoor-translation.z()); double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime()); double time = fibre->tShift(localPos, depth, chkFibre); #ifdef DebugLog LogDebug("HFShower") << "HFShower::getHits(SL): in " << name << " Z " <<zCoor << "(" << globalPos.z() << ") " << zFibre << " Trans " << translation << " Time " << tSlice << " " << time << "\n Direction " << momentumDir << " Local " << localMom; #endif // npe should be 0 int npe = 0; if(ok) npe = cherenkov->computeNPE(particleDef,pBeta,u,v,w, stepl,zFibre, dose, npeDose); std::vector<double> wavelength = cherenkov->getWL(); std::vector<double> momz = cherenkov->getMom(); for (int i = 0; i<npe; ++i) { hit.time = tSlice+time; hit.wavelength = wavelength[i]; hit.momentum = momz[i]; hit.position = globalPos; hits.push_back(hit); nHit++; } return hits; }
std::vector< HFShower::Hit > HFShower::getHits | ( | G4Step * | aStep | ) |
Definition at line 70 of file HFShower.cc.
References abs, applyFidCut, HFFibre::attLength(), cherenkov, chkFibre, HFCherenkov::computeNPE(), HFShower::Hit::depth, relval_parameters_module::energy, create_public_lumi_plots::exp, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, LogDebug, HFShower::Hit::momentum, mergeVDriftHistosByStation::name, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, HFFibreFiducial::PMTNumber(), HFShower::Hit::position, probMax, diffTwoXMLs::r1, diffTwoXMLs::r2, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), v, w(), and HFShower::Hit::wavelength.
Referenced by HCalSD::hitForFibre(), and FiberSD::ProcessHits().
{ std::vector<HFShower::Hit> hits; int nHit = 0; double edep = aStep->GetTotalEnergyDeposit(); double stepl = 0.; if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.) stepl = aStep->GetStepLength(); if ((edep == 0.) || (stepl == 0.)) { #ifdef DebugLog LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit; #endif return hits; } G4Track *aTrack = aStep->GetTrack(); const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle(); HFShower::Hit hit; double energy = aParticle->GetTotalEnergy(); double momentum = aParticle->GetTotalMomentum(); double pBeta = momentum / energy; double dose = 0.; int npeDose = 0; G4ThreeVector momentumDir = aParticle->GetMomentumDirection(); G4ParticleDefinition *particleDef = aTrack->GetDefinition(); G4StepPoint * preStepPoint = aStep->GetPreStepPoint(); G4ThreeVector globalPos = preStepPoint->GetPosition(); G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName(); //double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1]; double zv = std::abs(globalPos.z()) - gpar[4]; G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(), zv); G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()-> GetTopTransform().TransformAxis(momentumDir); // @@ Here the depth should be changed (Fibers all long in Geometry!) int depth = 1; int npmt = 0; bool ok = true; if (zv < 0. || zv > gpar[1]) { ok = false; // beyond the fiber in Z } if (ok && applyFidCut) { npmt = HFFibreFiducial::PMTNumber(globalPos); if (npmt <= 0) { ok = false; } else if (npmt > 24) { // a short fibre if (zv > gpar[0]) { depth = 2; } else { ok = false; } } #ifdef DebugLog edm::LogInfo("HFShower") << "HFShower: npmt " << npmt << " zv " << std::abs(pe_effect.z()) << ":" << gpar[4] << ":" << zv << ":" << gpar[0] << " ok " << ok << " depth " << depth; #endif } else { depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG! } G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation(); double u = localMom.x(); double v = localMom.y(); double w = localMom.z(); double zCoor = localPos.z(); double zFibre = (0.5*gpar[1]-zCoor-translation.z()); double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime()); double time = fibre->tShift(localPos, depth, chkFibre); #ifdef DebugLog LogDebug("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") " << zFibre << " Trans " << translation << " Time " << tSlice << " " << time << "\n Direction " << momentumDir << " Local " << localMom; #endif // Here npe should be 0 if there is no fiber (@@ M.K.) int npe = 1; std::vector<double> wavelength; std::vector<double> momz; if (!applyFidCut) { // _____ Tmp close of the cherenkov function if (ok) npe = cherenkov->computeNPE(particleDef,pBeta,u,v,w,stepl,zFibre,dose, npeDose); wavelength = cherenkov->getWL(); momz = cherenkov->getMom(); } // ^^^^^ End of Tmp close of the cherenkov function if(ok && npe>0) { for (int i = 0; i<npe; ++i) { double p=1.; if (!applyFidCut) p = fibre->attLength(wavelength[i]); double r1 = G4UniformRand(); double r2 = G4UniformRand(); #ifdef DebugLog LogDebug("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1 <<":" << exp(-p*zFibre) << " r2 " << r2 << ":" << probMax << " Survive: " << (r1 <= exp(-p*zFibre) && r2 <= probMax); #endif if (applyFidCut || chkFibre < 0 || (r1 <= exp(-p*zFibre) && r2 <= probMax)) { hit.depth = depth; hit.time = tSlice+time; if (!applyFidCut) { hit.wavelength = wavelength[i]; // Tmp hit.momentum = momz[i]; // Tmp } else { hit.wavelength = 300.; // Tmp hit.momentum = 1.; // Tmp } hit.position = globalPos; hits.push_back(hit); nHit++; } } } #ifdef DebugLog LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit; for (int i=0; i<nHit; ++i) LogDebug("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time " << hits[i].time << " Momentum " << hits[i].momentum << " Position " << hits[i].position; #endif return hits; }
int HFShower::indexFinder | ( | double | x, |
const std::vector< double > & | Fhist | ||
) | [private] |
Definition at line 659 of file HFShower.cc.
References debug, LogDebug, findQualityFiles::size, and launcher::step.
Referenced by compute().
{ // binary search in the vector of doubles int size = Fhist.size(); int curr = size / 2; int step = size / 4; int iter; int prevdir = 0; int actudir = 0; for (iter = 0; iter < size ; iter++) { if( curr >= size || curr < 1 ) LogWarning("FastCalorimetry") << " FamosHFShower::indexFinder - wrong current index = " << curr << " !!!" << std::endl; if ((x <= Fhist[curr]) && (x > Fhist[curr-1])) break; prevdir = actudir; if(x > Fhist[curr]) {actudir = 1;} else {actudir = -1;} if(prevdir * actudir < 0) { if(step > 1) step /= 2;} curr += actudir * step; if(curr > size) curr = size; else { if(curr < 1) {curr = 1;}} if(debug == 3) LogDebug("FastCalorimetry") << " indexFinder - end of iter." << iter << " curr, F[curr-1], F[curr] = " << curr << " " << Fhist[curr-1] << " " << Fhist[curr] << std::endl; } if(debug == 3) LogDebug("FastCalorimetry") << " indexFinder x = " << x << " found index = " << curr-1 << std::endl; return curr-1; }
void HFShower::initRun | ( | G4ParticleTable * | theParticleTable | ) |
Definition at line 340 of file HFShower.cc.
Referenced by HCalSD::initRun().
{// Define PDG codes /* emPDG = theParticleTable->FindParticle("e-")->GetPDGEncoding(); epPDG = theParticleTable->FindParticle("e+")->GetPDGEncoding(); gammaPDG = theParticleTable->FindParticle("gamma")->GetPDGEncoding(); #ifdef DebugLog edm::LogInfo("HFShower") << "HFShower: Particle codes" << " for e- = " << emPDG << " for e+ = " << epPDG << " for gamma = " << gammaPDG; #endif */ }
void HFShower::makeSteps | ( | int | nsteps | ) | [private] |
Definition at line 366 of file HFShower.cc.
References aloge, alpEM, alpHD, balanceEH, betEM, betHD, prof2calltree::count, debug, detector, e, eSpotSize, eStep, RandomEngine::flatShoot(), gam(), i, infinity, lamcurr, lamdepth, lamstep, LogDebug, nspots, random, rlamStep, groupFilesInBlocks::temp, tgamEM, tgamHD, theR1, theR2, theR3, transFactor, transParam, x, x0curr, x0depth, and detailsBasic3DVector::y.
Referenced by HFShower().
{ double sumes = 0.; double sum = 0.; std::vector<double> temp; if(debug) LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - " << " nsteps required : " << nsteps << std::endl; int count = 0; for (int i = 0; i < nsteps; i++) { double deplam = lamdepth[i] - 0.5 * lamstep[i]; double depx0 = x0depth[i] - 0.5 * lamstep[i] / x0curr[i]; double x = betEM * depx0; double y = betHD * deplam; if(debug == 2) LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps " << " - step " << i << " depx0, x = " << depx0 << ", " << x << " deplam, y = " << deplam << ", " << y << std::endl; double est = (part * betEM * gam(x,alpEM) * lamcurr[i] / (x0curr[i] * tgamEM) + (1.-part) * betHD * gam(y,alpHD) / tgamHD) * lamstep[i]; // protection ... if(est < 0.) { LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps " << " - negative step energy !!!" << std::endl; est = 0.; break ; } // for estimates only sum += est; int nPest = (int) (est * e / sum / eSpotSize) ; if(debug == 2) LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - nPoints estimate = " << nPest << std::endl; if(nPest <= 1 && count !=0 ) break; // good step - to proceed temp.push_back(est); sumes += est; rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge)) * deplam * transFactor); count ++; } // fluctuations in ECAL and re-distribution of remaining energy in HCAL if(detector[0] == 1 && count > 1) { double oldECALenergy = temp[0]; double oldHCALenergy = sumes - oldECALenergy ; double newECALenergy = 2. * sumes; for (int i = 0; newECALenergy > sumes && i < infinity; i++) newECALenergy = 2.* balanceEH * random->flatShoot() * oldECALenergy; if(debug == 2) LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps " << " ECAL fraction : old/new - " << oldECALenergy/sumes << "/" << newECALenergy/sumes << std::endl; temp[0] = newECALenergy; double newHCALenergy = sumes - newECALenergy; double newHCALreweight = newHCALenergy / oldHCALenergy; for (int i = 1; i < count; i++) { temp[i] *= newHCALreweight; } } // final re-normalization of the energy fractions double etot = 0.; for (int i = 0; i < count ; i++) { eStep.push_back(temp[i] * e / sumes ); nspots.push_back((int)(eStep[i]/eSpotSize)+1); etot += eStep[i]; if(debug) LogDebug("FastCalorimetry") << i << " xO and lamdepth at the end of step = " << x0depth[i] << " " << lamdepth[i] << " Estep func = " << eStep[i] << " Rstep = " << rlamStep[i] << " Nspots = " << nspots[i] << std::endl; } // The only step is in ECAL - let's make the size bigger ... if(count == 1 and detector[0] == 1) rlamStep[0] *= 2.; if(debug) { if(eStep[0] > 0.95 * e && detector[0] == 1) LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - " << "ECAL energy = " << eStep[0] << " out of total = " << e << std::endl; } }
double HFShower::transProb | ( | double | factor, |
double | R, | ||
double | r | ||
) | [inline, private] |
double HFShower::aloge [private] |
Definition at line 77 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::alpEM [private] |
Definition at line 72 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::alpHD [private] |
Definition at line 72 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
bool HFShower::applyFidCut [private] |
Definition at line 47 of file HFShower.h.
Referenced by getHits(), and HFShower().
double HFShower::balanceEH [private] |
Definition at line 117 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::betEM [private] |
Definition at line 72 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::betHD [private] |
Definition at line 72 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
HFCherenkov* HFShower::cherenkov [private] |
Definition at line 51 of file HFShower.h.
Referenced by getHits(), and HFShower().
int HFShower::chkFibre [private] |
Definition at line 54 of file HFShower.h.
Referenced by getHits(), and HFShower().
double HFShower::criticalEnergy [private] |
Definition at line 113 of file HFShower.h.
Referenced by HFShower().
double HFShower::depthStart [private] |
Definition at line 76 of file HFShower.h.
Referenced by HFShower().
double HFShower::depthStep [private] |
Definition at line 111 of file HFShower.h.
Referenced by HFShower().
std::vector<int> HFShower::detector [private] |
Definition at line 79 of file HFShower.h.
Referenced by compute(), HFShower(), and makeSteps().
double HFShower::e [private] |
Definition at line 96 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::eSpotSize [private] |
Definition at line 109 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
std::vector<double> HFShower::eStep [private] |
Definition at line 80 of file HFShower.h.
Referenced by compute(), and makeSteps().
HFFibre* HFShower::fibre [private] |
Definition at line 52 of file HFShower.h.
Referenced by getHits(), and HFShower().
std::vector<double> HFShower::gpar [private] |
Definition at line 56 of file HFShower.h.
Referenced by getHits(), and HFShower().
double HFShower::hcalDepthFactor [private] |
Definition at line 119 of file HFShower.h.
Referenced by compute(), and HFShower().
int HFShower::infinity [private] |
Definition at line 84 of file HFShower.h.
Referenced by compute(), and makeSteps().
double HFShower::lambdaEM [private] |
Definition at line 75 of file HFShower.h.
Referenced by HFShower().
double HFShower::lambdaHD [private] |
Definition at line 75 of file HFShower.h.
Referenced by HFShower().
std::vector<double> HFShower::lamcurr [private] |
Definition at line 82 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
std::vector<double> HFShower::lamdepth [private] |
Definition at line 82 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
std::vector<double> HFShower::lamstep [private] |
Definition at line 82 of file HFShower.h.
Referenced by compute(), HFShower(), and makeSteps().
std::vector<double> HFShower::lamtotal [private] |
Definition at line 82 of file HFShower.h.
Referenced by compute(), and HFShower().
int HFShower::lossesOpt [private] |
Definition at line 99 of file HFShower.h.
Referenced by compute(), and HFShower().
double HFShower::maxTRfactor [private] |
Definition at line 115 of file HFShower.h.
Referenced by compute(), and HFShower().
int HFShower::nDepthSteps [private] |
Definition at line 101 of file HFShower.h.
Referenced by HFShower().
std::vector<int> HFShower::nspots [private] |
Definition at line 79 of file HFShower.h.
Referenced by compute(), and makeSteps().
int HFShower::nTRsteps [private] |
Definition at line 103 of file HFShower.h.
Referenced by compute(), and HFShower().
int HFShower::onEcal [private] |
Definition at line 93 of file HFShower.h.
Referenced by HFShower().
double HFShower::part [private] |
Definition at line 72 of file HFShower.h.
double HFShower::probMax [private] |
Definition at line 55 of file HFShower.h.
Referenced by getHits(), and HFShower().
const RandomEngine* HFShower::random [private] |
Definition at line 122 of file HFShower.h.
Referenced by compute(), HFShower(), and makeSteps().
std::vector<double> HFShower::rlamStep [private] |
Definition at line 80 of file HFShower.h.
Referenced by compute(), and makeSteps().
double HFShower::tgamEM [private] |
Definition at line 72 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::tgamHD [private] |
Definition at line 72 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
const ECALProperties* HFShower::theECALproperties [private] |
Definition at line 67 of file HFShower.h.
Referenced by HFShower().
EcalHitMaker* HFShower::theGrid [private] |
Definition at line 87 of file HFShower.h.
Referenced by compute(), and HFShower().
HcalHitMaker* HFShower::theHcalHitMaker [private] |
Definition at line 90 of file HFShower.h.
Referenced by compute().
const HCALProperties* HFShower::theHCALproperties [private] |
Definition at line 68 of file HFShower.h.
Referenced by HFShower().
HDShowerParametrization* HFShower::theParam [private] |
Definition at line 64 of file HFShower.h.
Referenced by HFShower().
double HFShower::theR1 [private] |
Definition at line 71 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::theR2 [private] |
Definition at line 71 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::theR3 [private] |
Definition at line 71 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::transFactor [private] |
Definition at line 107 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::transParam [private] |
Definition at line 105 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
std::vector<double> HFShower::x0curr [private] |
Definition at line 81 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
std::vector<double> HFShower::x0depth [private] |
Definition at line 81 of file HFShower.h.
Referenced by HFShower(), and makeSteps().
double HFShower::x0EM [private] |
Definition at line 75 of file HFShower.h.
Referenced by HFShower().
double HFShower::x0HD [private] |
Definition at line 75 of file HFShower.h.
Referenced by HFShower().