|
|
|
#include <EMShower.h>
Public Member Functions | |
| void | compute () |
| Compute the shower longitudinal and lateral development. | |
| EMShower (const RandomEngine *engine, GammaFunctionGenerator *gamma, EMECALShowerParametrization *const myParam, std::vector< const RawParticle * > *const myPart, EcalHitMaker *const myGrid=NULL, PreshowerHitMaker *const myPreshower=NULL) | |
| double | getMaximumOfShower () const |
| get the depth of the centre of gravity of the shower(s) | |
| void | prepareSteps () |
| Computes the steps before the real compute. | |
| void | setGrid (EcalHitMaker *const myGrid) |
| set the grid address | |
| void | setHcal (HcalHitMaker *const myHcal) |
| set the HCAL address | |
| void | setPreshower (PreshowerHitMaker *const myPresh) |
| set the preshower address | |
| virtual | ~EMShower () |
Private 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 |
Private Member Functions | |
| double | deposit (double t, double a, double b, double dt) |
| double | deposit (double a, double b, double t) |
| double | gam (double x, double a) const |
| void | setIntervals (unsigned icomp, RadialInterval &rad) |
Private Attributes | |
| std::vector< double > | a |
| std::vector< double > | aSpot |
| std::vector< double > | b |
| std::vector< double > | bSpot |
| std::vector< std::vector < double > > | depositedEnergy |
| std::vector< double > | E |
| std::vector< double > | Etot |
| double | globalMaximum |
| bool | hasPreshower |
| double | innerDepth |
| std::vector< double > | maximumOfShower |
| std::vector< double > | meanDepth |
| GammaFunctionGenerator * | myGammaGenerator |
| Genfun::IncompleteGamma | myIncompleteGamma |
| unsigned int | nPart |
| unsigned | nSteps |
| double | outerDepth |
| std::vector< double > | photos |
| const RandomEngine * | random |
| Steps | steps |
| bool | stepsCalculated |
| std::vector< double > | T |
| const ECALProperties * | theECAL |
| EcalHitMaker * | theGrid |
| const HCALProperties * | theHCAL |
| HcalHitMaker * | theHcalHitMaker |
| const PreshowerLayer1Properties * | theLayer1 |
| const PreshowerLayer2Properties * | theLayer2 |
| std::vector< double > | theNumberOfSpots |
| EMECALShowerParametrization *const | theParam |
| std::vector< const RawParticle * > *const | thePart |
| PreshowerHitMaker * | thePreshower |
| std::vector< double > | Ti |
| double | totalEnergy |
| std::vector< double > | TSpot |
Definition at line 26 of file EMShower.h.
typedef std::pair<XYZPoint,double> EMShower::Spot [private] |
Definition at line 31 of file EMShower.h.
typedef std::pair<unsigned int, double> EMShower::Step [private] |
Definition at line 32 of file EMShower.h.
typedef Steps::const_iterator EMShower::step_iterator [private] |
Definition at line 34 of file EMShower.h.
typedef std::vector<Step> EMShower::Steps [private] |
Definition at line 33 of file EMShower.h.
typedef math::XYZVector EMShower::XYZPoint [private] |
Definition at line 29 of file EMShower.h.
| EMShower::EMShower | ( | const RandomEngine * | engine, |
| GammaFunctionGenerator * | gamma, | ||
| EMECALShowerParametrization *const | myParam, | ||
| std::vector< const RawParticle * > *const | myPart, | ||
| EcalHitMaker *const | myGrid = NULL, |
||
| PreshowerHitMaker *const | myPreshower = NULL |
||
| ) |
Definition at line 17 of file EMShower.cc.
References a, aSpot, b, bSpot, EMECALShowerParametrization::correlationAlphaT(), ECALProperties::criticalEnergy(), ExpressReco_HICollisions_FallBack::e, E, EMECALShowerParametrization::ecalProperties(), Etot, funct::exp(), RandomEngine::gaussShoot(), globalMaximum, hasPreshower, EMECALShowerParametrization::hcalProperties(), i, EMECALShowerParametrization::layer1Properties(), EMECALShowerParametrization::layer2Properties(), ECALProperties::lightCollectionEfficiency(), funct::log(), maximumOfShower, EMECALShowerParametrization::meanAlpha(), EMECALShowerParametrization::meanAlphaSpot(), meanDepth, EMECALShowerParametrization::meanLnAlpha(), EMECALShowerParametrization::meanLnT(), EMECALShowerParametrization::meanT(), EMECALShowerParametrization::meanTSpot(), nPart, EMECALShowerParametrization::nSpots(), NULL, photos, ECALProperties::photoStatistics(), random, EMECALShowerParametrization::sigmaLnAlpha(), EMECALShowerParametrization::sigmaLnT(), mathSSE::sqrt(), stepsCalculated, T, theECAL, theHCAL, theLayer1, theLayer2, theNumberOfSpots, theParam, thePart, Ti, totalEnergy, and TSpot.
: theParam(myParam), thePart(myPart), theGrid(myGrid), thePreshower(myPresh), random(engine), myGammaGenerator(gamma) { // Get the Famos Histos pointer // myHistos = Histos::instance(); // myGammaGenerator = GammaFunctionGenerator::instance(); stepsCalculated=false; hasPreshower = myPresh!=NULL; theECAL = myParam->ecalProperties(); theHCAL = myParam->hcalProperties(); theLayer1 = myParam->layer1Properties(); theLayer2 = myParam->layer2Properties(); double fotos = theECAL->photoStatistics() * theECAL->lightCollectionEfficiency(); nPart = thePart->size(); totalEnergy = 0.; globalMaximum = 0.; double meanDepth=0.; // Initialize the shower parameters for each particle for ( unsigned int i=0; i<nPart; ++i ) { // std::cout << " AAA " << *(*thePart)[i] << std::endl; // The particle and the shower energy Etot.push_back(0.); E.push_back(((*thePart)[i])->e()); totalEnergy+=E[i]; double lny = std::log ( E[i] / theECAL->criticalEnergy() ); // Average and Sigma for T and alpha double theMeanT = myParam->meanT(lny); double theMeanAlpha = myParam->meanAlpha(lny); double theMeanLnT = myParam->meanLnT(lny); double theMeanLnAlpha = myParam->meanLnAlpha(lny); double theSigmaLnT = myParam->sigmaLnT(lny); double theSigmaLnAlpha = myParam->sigmaLnAlpha(lny); // The correlation matrix double theCorrelation = myParam->correlationAlphaT(lny); double rhop = std::sqrt( (1.+theCorrelation)/2. ); double rhom = std::sqrt( (1.-theCorrelation)/2. ); // The number of spots in ECAL / HCAL theNumberOfSpots.push_back(myParam->nSpots(E[i])); // theNumberOfSpots.push_back(myParam->nSpots(E[i])*spotFraction); //theNumberOfSpots = random->poissonShoot(myParam->nSpots(myPart->e())); // Photo-statistics photos.push_back(E[i] * fotos); // The longitudinal shower development parameters // Fluctuations of alpha, T and beta double z1=0.; double z2=0.; double aa=0.; // Protect against too large fluctuations (a < 1) for small energies while ( aa <= 1. ) { z1 = random->gaussShoot(0.,1.); z2 = random->gaussShoot(0.,1.); aa = std::exp(theMeanLnAlpha + theSigmaLnAlpha * (z1*rhop-z2*rhom)); } a.push_back(aa); T.push_back(std::exp(theMeanLnT + theSigmaLnT * (z1*rhop+z2*rhom))); b.push_back((a[i]-1.)/T[i]); maximumOfShower.push_back((a[i]-1.)/b[i]); globalMaximum += maximumOfShower[i]*E[i]; meanDepth += a[i]/b[i]*E[i]; // std::cout << " Adding max " << maximumOfShower[i] << " " << E[i] << " " <<maximumOfShower[i]*E[i] << std::endl; // std::cout << std::setw(8) << std::setprecision(5) << " a /b " << a[i] << " " << b[i] << std::endl; Ti.push_back( a[i]/b[i] * (std::exp(theMeanLnAlpha)-1.) / std::exp(theMeanLnAlpha)); // The parameters for the number of energy spots TSpot.push_back(theParam->meanTSpot(theMeanT)); aSpot.push_back(theParam->meanAlphaSpot(theMeanAlpha)); bSpot.push_back((aSpot[i]-1.)/TSpot[i]); // myHistos->fill("h7000",a[i]); // myHistos->fill("h7002",E[i],a[i]); } // std::cout << " PS1 : " << myGrid->ps1TotalX0() // << " PS2 : " << myGrid->ps2TotalX0() // << " ECAL : " << myGrid->ecalTotalX0() // << " HCAL : " << myGrid->hcalTotalX0() // << " Offset : " << myGrid->x0DepthOffset() // << std::endl; globalMaximum/=totalEnergy; meanDepth/=totalEnergy; // std::cout << " Total Energy " << totalEnergy << " Global max " << globalMaximum << std::endl; }
| virtual EMShower::~EMShower | ( | ) | [inline, virtual] |
Definition at line 45 of file EMShower.h.
{;}
| void EMShower::compute | ( | ) |
Compute the shower longitudinal and lateral development.
Definition at line 244 of file EMShower.cc.
References a, EcalHitMaker::addHit(), HcalHitMaker::addHit(), PreshowerHitMaker::addHit(), EcalHitMaker::addHitDepth(), RadialInterval::addInterval(), aSpot, b, bSpot, RadialInterval::compute(), deposit(), depositedEnergy, ExpressReco_HICollisions_FallBack::detector, dt, E, patCandidatesForDimuonsSequences_cff::ecal, Etot, RandomEngine::flatShoot(), gam(), RandomEngine::gaussShoot(), RadialInterval::getNumberOfSpots(), EcalHitMaker::getPads(), RadialInterval::getSpotEnergy(), RadialInterval::getUmax(), RadialInterval::getUmin(), EcalHitMaker::getX0back(), hasPreshower, patCandidatesForDimuonsSequences_cff::hcal, HCALProperties::hOverPi(), i, ECALProperties::lightCollectionUniformity(), M_PI, meanDepth, PreshowerLayer2Properties::mipsPerGeV(), PreshowerLayer1Properties::mipsPerGeV(), myGammaGenerator, RadialInterval::nIntervals(), nPart, nSteps, NULL, EMECALShowerParametrization::p(), phi, photos, RandomEngine::poissonShoot(), prepareSteps(), random, EMECALShowerParametrization::rC(), EMECALShowerParametrization::rT(), HcalHitMaker::setDepth(), setIntervals(), GammaFunctionGenerator::setParameters(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), PreshowerHitMaker::setSpotEnergy(), GammaFunctionGenerator::shoot(), HCALProperties::spotFraction(), mathSSE::sqrt(), ntuplemaker::status, steps, stepsCalculated, matplotRender::t, theECAL, theGrid, theHCAL, theHcalHitMaker, theLayer1, theLayer2, theNumberOfSpots, theParam, thePreshower, and Ti.
Referenced by CalorimetryManager::EMShowerSimulation().
{
double t = 0.;
double dt = 0.;
if(!stepsCalculated) prepareSteps();
// Prepare the grids in EcalHitMaker
// theGrid->setInnerAndOuterDepth(innerDepth,outerDepth);
float pstot=0.;
float ps2tot=0.;
float ps1tot=0.;
bool status=false;
// double E1 = 0.; // Energy layer 1
// double E2 = 0.; // Energy layer 2
// double n1 = 0.; // #mips layer 1
// double n2 = 0.; // #mips layer 2
// double E9 = 0.; // Energy ECAL
// Loop over all segments for the longitudinal development
for (unsigned iStep=0; iStep<nSteps; ++iStep ) {
// The length of the shower in this segment
dt = steps[iStep].second;
// std::cout << " Detector " << steps[iStep].first << " t " << t << " " << dt << std::endl;
// The elapsed length
t += dt;
// In what detector are we ?
unsigned detector=steps[iStep].first;
bool presh1 = detector==0;
bool presh2 = detector==1;
bool ecal = detector==2;
bool hcal = detector==3;
bool vfcal = detector==4;
bool gap = detector==5;
// Temporary. Will be removed
if ( theHCAL==NULL) hcal=false;
// Keep only ECAL for now
if ( vfcal ) continue;
// Nothing to do in the gap
if( gap ) continue;
// cout << " t = " << t << endl;
// Build the grid of crystals at this ECAL depth
// Actually, it might be useful to check if this grid is empty or not.
// If it is empty (because no crystal at this depth), it is of no use
// (and time consuming) to generate the spots
// middle of the step
double tt = t-0.5*dt;
double realTotalEnergy=0.;
for ( unsigned int i=0; i<nPart; ++i ) {
realTotalEnergy += depositedEnergy[iStep][i]*E[i];
}
// std::cout << " Step " << tt << std::endl;
// std::cout << "ecal " << ecal << " hcal " << hcal <<std::endl;
// If the amount of energy is greater than 1 MeV, make a new grid
// otherwise put in the previous one.
bool usePreviousGrid=(realTotalEnergy<0.001);
// If the amount of energy is greater than 1 MeV, make a new grid
// otherwise put in the previous one.
// If less than 1 kEV. Just skip
if(iStep>2&&realTotalEnergy<0.000001) continue;
if (ecal && !usePreviousGrid)
{
status=theGrid->getPads(meanDepth[iStep]);
}
if (hcal)
{
status=theHcalHitMaker->setDepth(tt);
}
if((ecal ||hcal) && !status) continue;
bool detailedShowerTail=false;
// check if a detailed treatment of the rear leakage should be applied
if(ecal && !usePreviousGrid)
{
detailedShowerTail=(t-dt > theGrid->getX0back());
}
// The particles of the shower are processed in parallel
for ( unsigned int i=0; i<nPart; ++i ) {
// double Edepo=deposit(t,a[i],b[i],dt);
// integration of the shower profile between t-dt and t
double dE = (!hcal)? depositedEnergy[iStep][i]:1.-deposit(a[i],b[i],t-dt);
// no need to do the full machinery if there is ~nothing to distribute)
if(dE*E[i]<0.000001) continue;
if(detailedShowerTail)
{
myGammaGenerator->setParameters(floor(a[i]+0.5),b[i],t-dt);
}
// The number of energy spots (or mips)
double nS = 0;
// ECAL case : Account for photostatistics and long'al non-uniformity
if (ecal) {
dE = random->poissonShoot(dE*photos[i])/photos[i];
double z0 = random->gaussShoot(0.,1.);
dE *= 1. + z0*theECAL->lightCollectionUniformity();
// Expected spot number
nS = ( theNumberOfSpots[i] * gam(bSpot[i]*tt,aSpot[i])
* bSpot[i] * dt
/ tgamma(aSpot[i]) );
// Preshower : Expected number of mips + fluctuation
}
else if ( hcal ) {
nS = ( theNumberOfSpots[i] * gam(bSpot[i]*tt,aSpot[i])
* bSpot[i] * dt
/ tgamma(aSpot[i]))* theHCAL->spotFraction();
double nSo = nS ;
nS = random->poissonShoot(nS);
// 'Quick and dirty' fix (but this line should be better removed):
if( nSo > 0. && nS/nSo < 10.) dE *= nS/nSo;
// if(true)
// {
// std::cout << " theHCAL->spotFraction = " <<theHCAL->spotFraction() <<std::endl;
// std::cout << " nSpot Ecal : " << nSo/theHCAL->spotFraction() << " Final " << nS << std::endl;
// }
}
else if ( presh1 ) {
nS = random->poissonShoot(dE*E[i]*theLayer1->mipsPerGeV());
// std::cout << " dE *E[i] (1)" << dE*E[i] << " " << dE*E[i]*theLayer1->mipsPerGeV() << " "<< nS << std::endl;
pstot+=dE*E[i];
ps1tot+=dE*E[i];
dE = nS/(E[i]*theLayer1->mipsPerGeV());
// E1 += dE*E[i];
// n1 += nS;
// if (presh2) { E2 += SpotEnergy; ++n2; }
} else if ( presh2 ) {
nS = random->poissonShoot(dE*E[i]*theLayer2->mipsPerGeV());
// std::cout << " dE *E[i] (2) " << dE*E[i] << " " << dE*E[i]*theLayer2->mipsPerGeV() << " "<< nS << std::endl;
pstot+=dE*E[i];
ps2tot+=dE*E[i];
dE = nS/(E[i]*theLayer2->mipsPerGeV());
// E2 += dE*E[i];
// n2 += nS;
}
// myHistos->fill("h100",t,dE);
// The lateral development parameters
// Energy of the spots
double eSpot = (nS>0.) ? dE/nS : 0.;
double SpotEnergy=eSpot*E[i];
if(hasPreshower&&(presh1||presh2)) thePreshower->setSpotEnergy(0.00009);
if(hcal)
{
SpotEnergy*=theHCAL->hOverPi();
theHcalHitMaker->setSpotEnergy(SpotEnergy);
}
// Poissonian fluctuations for the number of spots
// int nSpot = random->poissonShoot(nS);
int nSpot = (int)(nS+0.5);
// Fig. 11 (right) *** Does not match.
// myHistos->fill("h101",t,(double)nSpot/theNumberOfSpots);
//double taui = t/T;
double taui = tt/Ti[i];
double proba = theParam->p(taui,E[i]);
double theRC = theParam->rC(taui,E[i]);
double theRT = theParam->rT(taui,E[i]);
// Fig. 10
// myHistos->fill("h300",taui,theRC);
// myHistos->fill("h301",taui,theRT);
// myHistos->fill("h302",taui,proba);
double dSpotsCore =
random->gaussShoot(proba*nSpot,std::sqrt(proba*(1.-proba)*nSpot));
if(dSpotsCore<0) dSpotsCore=0;
unsigned nSpots_core = (unsigned)(dSpotsCore+0.5);
unsigned nSpots_tail = ((unsigned)nSpot>nSpots_core) ? nSpot-nSpots_core : 0;
for(unsigned icomp=0;icomp<2;++icomp)
{
double theR=(icomp==0) ? theRC : theRT ;
unsigned ncompspots=(icomp==0) ? nSpots_core : nSpots_tail;
RadialInterval radInterval(theR,ncompspots,SpotEnergy,random);
if(ecal)
{
if(icomp==0)
{
setIntervals(icomp,radInterval);
}
else
{
setIntervals(icomp,radInterval);
}
}
else
{
radInterval.addInterval(100.,1.);// 100% of the spots
}
radInterval.compute();
// irad = 0 : central circle; irad=1 : outside
unsigned nrad=radInterval.nIntervals();
for(unsigned irad=0;irad<nrad;++irad)
{
double spote=radInterval.getSpotEnergy(irad);
if(ecal) theGrid->setSpotEnergy(spote);
if(hcal) theHcalHitMaker->setSpotEnergy(spote);
unsigned nradspots=radInterval.getNumberOfSpots(irad);
double umin=radInterval.getUmin(irad);
double umax=radInterval.getUmax(irad);
// Go for the lateral development
for ( unsigned ispot=0; ispot<nradspots; ++ispot )
{
double z3=random->flatShoot(umin,umax);
double ri=theR * std::sqrt(z3/(1.-z3)) ;
//Fig. 12
/*
if ( 2. < t && t < 3. )
myHistos->fill("h401",ri,1./1000.*eSpot/dE/0.2);
if ( 6. < t && t < 7. )
myHistos->fill("h402",ri,1./1000.*eSpot/dE/0.2);
if ( 19. < t && t < 20. )
myHistos->fill("h403",ri,1./1000.*eSpot/dE/0.2);
*/
// Fig. 13 (top)
// myHistos->fill("h400",ri,1./1000.*eSpot/0.2);
// Generate phi
double phi = 2.*M_PI*random->flatShoot();
// Add the hit in the crystal
// if( ecal ) theGrid->addHit(ri*theECAL->moliereRadius(),phi);
// Now the *moliereRadius is done in EcalHitMaker
if ( ecal )
{
if(detailedShowerTail)
{
// std::cout << "About to call addHitDepth " << std::endl;
double depth;
do
{
depth=myGammaGenerator->shoot();
}
while(depth>t);
theGrid->addHitDepth(ri,phi,depth);
// std::cout << " Done " << std::endl;
}
else
theGrid->addHit(ri,phi);
}
else if (hasPreshower&&presh1) thePreshower->addHit(ri,phi,1);
else if (hasPreshower&&presh2) thePreshower->addHit(ri,phi,2);
else if (hcal)
{
// std::cout << " About to add a spot in the HCAL" << status << std::endl;
theHcalHitMaker->addHit(ri,phi);
// std::cout << " Added a spot in the HCAL" << status << std::endl;
}
// if (ecal) E9 += SpotEnergy;
// if (presh1) { E1 += SpotEnergy; ++n1; }
// if (presh2) { E2 += SpotEnergy; ++n2; }
Etot[i] += spote;
}
}
}
// std::cout << " Done with the step " << std::endl;
// The shower!
//myHistos->fill("h500",theSpot.z(),theSpot.perp());
}
// std::cout << " nPart " << nPart << std::endl;
}
// std::cout << " Finshed the step loop " << std::endl;
// myHistos->fill("h500",E1+0.7*E2,E9);
// myHistos->fill("h501",n1+0.7*n2,E9);
// myHistos->fill("h400",n1);
// myHistos->fill("h401",n2);
// myHistos->fill("h402",E9+E1+0.7*E2);
// if(!standalone)theGrid->printGrid();
double Etotal=0.;
for(unsigned i=0;i<nPart;++i)
{
// myHistos->fill("h10",Etot[i]);
Etotal+=Etot[i];
}
// myHistos->fill("h20",Etotal);
// if(thePreshower)
// std::cout << " PS " << thePreshower->layer1Calibrated() << " " << thePreshower->layer2Calibrated() << " " << thePreshower->totalCalibrated() << " " << ps1tot << " " <<ps2tot << " " << pstot << std::endl;
}
| double EMShower::deposit | ( | double | t, |
| double | a, | ||
| double | b, | ||
| double | dt | ||
| ) | [private] |
Definition at line 607 of file EMShower.cc.
References dt, myIncompleteGamma, query::result, and matplotRender::t.
Referenced by compute(), and prepareSteps().
{
myIncompleteGamma.a().setValue(a);
double b1=b*(t-dt);
double b2=b*t;
double result = 0.;
double rb1=(b1!=0.) ? myIncompleteGamma(b1) : 0.;
double rb2=(b2!=0.) ? myIncompleteGamma(b2) : 0.;
result = (rb2-rb1);
return result;
}
| double EMShower::deposit | ( | double | a, |
| double | b, | ||
| double | t | ||
| ) | [private] |
Definition at line 648 of file EMShower.cc.
References ExpressReco_HICollisions_FallBack::e, myIncompleteGamma, query::result, and matplotRender::t.
{
// std::cout << " Deposit " << std::endl;
myIncompleteGamma.a().setValue(a);
double b2=b*t;
double result = 0.;
if(fabs(b2) < 1.e-9 ) b2 = 1.e-9;
result=myIncompleteGamma(b2);
// std::cout << " deposit t = " << t << " " << result <<std::endl;
return result;
}
| double EMShower::gam | ( | double | x, |
| double | a | ||
| ) | const [private] |
Definition at line 571 of file EMShower.cc.
References funct::exp(), and funct::pow().
Referenced by compute().
| double EMShower::getMaximumOfShower | ( | ) | const [inline] |
get the depth of the centre of gravity of the shower(s)
get the depth of the maximum of the shower
Definition at line 57 of file EMShower.h.
References globalMaximum.
Referenced by CalorimetryManager::EMShowerSimulation().
{return globalMaximum;}
| void EMShower::prepareSteps | ( | ) |
Computes the steps before the real compute.
Definition at line 123 of file EMShower.cc.
References a, b, deposit(), depositedEnergy, dt, E, EcalHitMaker::ecalTotalX0(), EcalHitMaker::hcalTotalX0(), i, innerDepth, meanDepth, nPart, nSteps, evf::evtn::offset(), outerDepth, EcalHitMaker::ps1TotalX0(), EcalHitMaker::ps2eeTotalX0(), EcalHitMaker::ps2TotalX0(), ExpressReco_HICollisions_FallBack::step, steps, stepsCalculated, matplotRender::t, theGrid, EcalHitMaker::totalX0(), and EcalHitMaker::x0DepthOffset().
Referenced by compute().
{
// TimeMe theT("EMShower::compute");
// Determine the longitudinal intervals
// std::cout << " EMShower compute" << std::endl;
double dt;
double radlen;
int stps;
int first_Ecal_step=0;
int last_Ecal_step=0;
// The maximum is in principe 8 (with 5X0 steps in the ECAL)
steps.reserve(20);
// std::cout << " PS1 : " << theGrid->ps1TotalX0()
// << " PS2 : " << theGrid->ps2TotalX0()
// << " ECAL : " << theGrid->ecalTotalX0()
// << " HCAL : " << theGrid->hcalTotalX0()
// << " Offset : " << theGrid->x0DepthOffset()
// << std::endl;
radlen = -theGrid->x0DepthOffset();
// Preshower Layer 1
radlen += theGrid->ps1TotalX0();
if ( radlen > 0. ) {
steps.push_back(Step(0,radlen));
radlen = 0.;
}
// Preshower Layer 2
radlen += theGrid->ps2TotalX0();
if ( radlen > 0. ) {
steps.push_back(Step(1,radlen));
radlen = 0.;
}
// add a step between preshower and ee
radlen += theGrid->ps2eeTotalX0();
if ( radlen > 0.) {
steps.push_back(Step(5,radlen));
radlen = 0.;
}
// ECAL
radlen += theGrid->ecalTotalX0();
if ( radlen > 0. ) {
stps=(int)((radlen+2.5)/5.);
// stps=(int)((radlen+.5)/1.);
if ( stps == 0 ) stps = 1;
dt = radlen/(double)stps;
Step step(2,dt);
first_Ecal_step=steps.size();
for ( int ist=0; ist<stps; ++ist )
steps.push_back(step);
last_Ecal_step=steps.size()-1;
radlen = 0.;
}
// I should had a gap here !
// HCAL
radlen += theGrid->hcalTotalX0();
if ( radlen > 0. ) {
double dtFrontHcal=theGrid->totalX0()-theGrid->hcalTotalX0();
// One single step for the full HCAL
if(dtFrontHcal<30.)
{
dt=30.-dtFrontHcal;
Step step(3,dt);
steps.push_back(step);
}
}
nSteps=steps.size();
if(nSteps==0) return;
double ESliceTot=0.;
double MeanDepth=0.;
depositedEnergy.resize(nSteps);
meanDepth.resize(nSteps);
double t=0.;
int offset=0;
for(unsigned iStep=0;iStep<nSteps;++iStep)
{
ESliceTot=0.;
MeanDepth=0.;
double realTotalEnergy=0;
dt=steps[iStep].second;
t+=dt;
for ( unsigned int i=0; i<nPart; ++i ) {
depositedEnergy[iStep].push_back(deposit(t,a[i],b[i],dt));
ESliceTot +=depositedEnergy[iStep][i];
MeanDepth += deposit(t,a[i]+1.,b[i],dt)/b[i]*a[i];
realTotalEnergy+=depositedEnergy[iStep][i]*E[i];
}
if( ESliceTot > 0. ) // can happen for the shower tails; this depth will be skipped anyway
MeanDepth/=ESliceTot;
else
MeanDepth=t-dt;
meanDepth[iStep]=MeanDepth;
if(realTotalEnergy<0.001)
{
offset-=1;
}
}
innerDepth=meanDepth[first_Ecal_step];
if(last_Ecal_step+offset>=0)
outerDepth=meanDepth[last_Ecal_step+offset];
else
outerDepth=innerDepth;
stepsCalculated=true;
}
| void EMShower::setGrid | ( | EcalHitMaker *const | myGrid | ) | [inline] |
set the grid address
Definition at line 60 of file EMShower.h.
References theGrid.
Referenced by CalorimetryManager::EMShowerSimulation().
{ theGrid=myGrid;}
| void EMShower::setHcal | ( | HcalHitMaker *const | myHcal | ) |
set the HCAL address
Definition at line 642 of file EMShower.cc.
References theHcalHitMaker.
Referenced by CalorimetryManager::EMShowerSimulation().
{
theHcalHitMaker = myHcal;
}
| void EMShower::setIntervals | ( | unsigned | icomp, |
| RadialInterval & | rad | ||
| ) | [private] |
Definition at line 619 of file EMShower.cc.
References RadialInterval::addInterval(), EMECALShowerParametrization::getCoreIntervals(), EMECALShowerParametrization::getTailIntervals(), and theParam.
Referenced by compute().
{
// std::cout << " Got the pointer " << std::endl;
const std::vector<double>& myValues((icomp)?theParam->getTailIntervals():theParam->getCoreIntervals());
// std::cout << " Got the vector " << myValues.size () << std::endl;
unsigned nvals=myValues.size()/2;
for(unsigned iv=0;iv<nvals;++iv)
{
// std::cout << myValues[2*iv] << " " << myValues[2*iv+1] <<std::endl;
rad.addInterval(myValues[2*iv],myValues[2*iv+1]);
}
}
| void EMShower::setPreshower | ( | PreshowerHitMaker *const | myPresh | ) |
set the preshower address
Definition at line 632 of file EMShower.cc.
References hasPreshower, NULL, and thePreshower.
Referenced by CalorimetryManager::EMShowerSimulation().
{
if(myPresh!=NULL)
{
thePreshower = myPresh;
hasPreshower=true;
}
}
std::vector<double> EMShower::a [private] |
Definition at line 101 of file EMShower.h.
Referenced by compute(), EMShower(), and prepareSteps().
std::vector<double> EMShower::aSpot [private] |
Definition at line 105 of file EMShower.h.
Referenced by compute(), and EMShower().
std::vector<double> EMShower::b [private] |
Definition at line 102 of file EMShower.h.
Referenced by compute(), EMShower(), and prepareSteps().
std::vector<double> EMShower::bSpot [private] |
Definition at line 106 of file EMShower.h.
Referenced by compute(), and EMShower().
std::vector<std::vector<double> > EMShower::depositedEnergy [private] |
Definition at line 112 of file EMShower.h.
Referenced by compute(), and prepareSteps().
std::vector<double> EMShower::E [private] |
Definition at line 98 of file EMShower.h.
Referenced by compute(), EMShower(), and prepareSteps().
std::vector<double> EMShower::Etot [private] |
Definition at line 97 of file EMShower.h.
Referenced by compute(), and EMShower().
double EMShower::globalMaximum [private] |
Definition at line 116 of file EMShower.h.
Referenced by EMShower(), and getMaximumOfShower().
bool EMShower::hasPreshower [private] |
Definition at line 135 of file EMShower.h.
Referenced by compute(), EMShower(), and setPreshower().
double EMShower::innerDepth [private] |
Definition at line 114 of file EMShower.h.
Referenced by prepareSteps().
std::vector<double> EMShower::maximumOfShower [private] |
Definition at line 111 of file EMShower.h.
Referenced by EMShower().
std::vector<double> EMShower::meanDepth [private] |
Definition at line 113 of file EMShower.h.
Referenced by compute(), EMShower(), and prepareSteps().
GammaFunctionGenerator* EMShower::myGammaGenerator [private] |
Definition at line 145 of file EMShower.h.
Referenced by compute().
Genfun::IncompleteGamma EMShower::myIncompleteGamma [private] |
Definition at line 139 of file EMShower.h.
Referenced by deposit().
unsigned int EMShower::nPart [private] |
Definition at line 93 of file EMShower.h.
Referenced by compute(), EMShower(), and prepareSteps().
unsigned EMShower::nSteps [private] |
Definition at line 122 of file EMShower.h.
Referenced by compute(), and prepareSteps().
double EMShower::outerDepth [private] |
Definition at line 114 of file EMShower.h.
Referenced by prepareSteps().
std::vector<double> EMShower::photos [private] |
Definition at line 99 of file EMShower.h.
Referenced by compute(), and EMShower().
const RandomEngine* EMShower::random [private] |
Definition at line 142 of file EMShower.h.
Referenced by compute(), and EMShower().
Steps EMShower::steps [private] |
Definition at line 121 of file EMShower.h.
Referenced by compute(), and prepareSteps().
bool EMShower::stepsCalculated [private] |
Definition at line 123 of file EMShower.h.
Referenced by compute(), EMShower(), and prepareSteps().
std::vector<double> EMShower::T [private] |
Definition at line 100 of file EMShower.h.
Referenced by EMShower().
const ECALProperties* EMShower::theECAL [private] |
Definition at line 86 of file EMShower.h.
Referenced by compute(), and EMShower().
EcalHitMaker* EMShower::theGrid [private] |
Definition at line 126 of file EMShower.h.
Referenced by compute(), prepareSteps(), and setGrid().
const HCALProperties* EMShower::theHCAL [private] |
Definition at line 87 of file EMShower.h.
Referenced by compute(), and EMShower().
HcalHitMaker* EMShower::theHcalHitMaker [private] |
Definition at line 132 of file EMShower.h.
const PreshowerLayer1Properties* EMShower::theLayer1 [private] |
Definition at line 88 of file EMShower.h.
Referenced by compute(), and EMShower().
const PreshowerLayer2Properties* EMShower::theLayer2 [private] |
Definition at line 89 of file EMShower.h.
Referenced by compute(), and EMShower().
std::vector<double> EMShower::theNumberOfSpots [private] |
Definition at line 96 of file EMShower.h.
Referenced by compute(), and EMShower().
EMECALShowerParametrization* const EMShower::theParam [private] |
Definition at line 83 of file EMShower.h.
Referenced by compute(), EMShower(), and setIntervals().
std::vector<const RawParticle*>* const EMShower::thePart [private] |
Definition at line 92 of file EMShower.h.
Referenced by EMShower().
PreshowerHitMaker* EMShower::thePreshower [private] |
Definition at line 129 of file EMShower.h.
Referenced by compute(), and setPreshower().
std::vector<double> EMShower::Ti [private] |
Definition at line 103 of file EMShower.h.
Referenced by compute(), and EMShower().
double EMShower::totalEnergy [private] |
Definition at line 118 of file EMShower.h.
Referenced by EMShower().
std::vector<double> EMShower::TSpot [private] |
Definition at line 104 of file EMShower.h.
Referenced by EMShower().
1.7.3