#include <GflashHadronShowerProfile.h>

Inheritance diagram for GflashHadronShowerProfile:

Public Member Functions
std::vector< GflashHit > &	getGflashHitList ()
GflashShowino *	getGflashShowino ()
	GflashHadronShowerProfile (edm::ParameterSet parSet)
void	hadronicParameterization ()
void	initialize (int showerType, double energy, double globalTime, double charge, Gflash3Vector &position, Gflash3Vector &momentum)
virtual void	loadParameters ()
virtual	~GflashHadronShowerProfile ()
Protected Member Functions
double	depthScale (double ssp, double ssp0, double length)
void	doCholeskyReduction (double cc, double vv, const int ndim)
double	fLnE1 (double einc, const double *par)
double	fTanh (double einc, const double *par)
double	gammaProfile (double alpha, double beta, double depth, double lengthUnit)
void	getFluctuationVector (double lowTriangle, double correlationVector)
int	getNumberOfSpots (Gflash::CalorimeterNumber kCalor)
double	hoProfile (double pathLength, double refDepth)
Gflash3Vector	locateHitPosition (GflashTrajectoryPoint &point, double lateralArm)
double	longitudinalProfile ()
double	medianLateralArm (double depth, Gflash::CalorimeterNumber kCalor)
void	setEnergyScale (double einc, Gflash3Vector ssp)
double	twoGammaProfile (double *par, double depth, Gflash::CalorimeterNumber kIndex)
Protected Attributes
double	averageSpotEnergy [Gflash::kNumberCalorimeter]
double	energyScale [Gflash::kNumberCalorimeter]
double	lateralPar [Gflash::kNumberCalorimeter][Gflash::Nrpar]
double	longEcal [Gflash::NPar]
double	longHcal [Gflash::NPar]
double	theBField
bool	theGflashHcalOuter
std::vector< GflashHit >	theGflashHitList
GflashHistogram *	theHisto
edm::ParameterSet	theParSet
GflashShowino *	theShowino

Detailed Description

Definition at line 15 of file GflashHadronShowerProfile.h.

Constructor & Destructor Documentation

GflashHadronShowerProfile::GflashHadronShowerProfile ( edm::ParameterSet parSet )

Definition at line 17 of file GflashHadronShowerProfile.cc.

References edm::ParameterSet::getParameter(), instance, theBField, theGflashHcalOuter, theHisto, and theShowino.

                                                                           : theParSet(parSet)
{
  theBField = parSet.getParameter<double>("bField");
  theGflashHcalOuter = parSet.getParameter<bool>("GflashHcalOuter");

  theShowino = new GflashShowino();
  theHisto = GflashHistogram::instance();
}

GflashHadronShowerProfile::~GflashHadronShowerProfile ( ) [virtual]

Definition at line 26 of file GflashHadronShowerProfile.cc.

References theShowino.

{
  if(theShowino) delete theShowino;
}

Member Function Documentation

double GflashHadronShowerProfile::depthScale	(	double	ssp,
		double	ssp0,
		double	length
	)		`[protected]`

Definition at line 497 of file GflashHadronShowerProfile.cc.

References funct::pow().

Referenced by GflashAntiProtonShowerProfile::loadParameters(), GflashProtonShowerProfile::loadParameters(), GflashKaonPlusShowerProfile::loadParameters(), GflashPiKShowerProfile::loadParameters(), and GflashKaonMinusShowerProfile::loadParameters().

                                                                                   {
  double func = 0.0;
  if(length>0.0) func = std::pow((ssp-ssp0)/length,2.0);
  return func;
}

void GflashHadronShowerProfile::doCholeskyReduction	(	double **	cc,
		double **	vv,
		const int	ndim
	)		`[protected]`

Definition at line 406 of file GflashHadronShowerProfile.cc.

References i, j, gen::k, and mathSSE::sqrt().

Referenced by getFluctuationVector().

                                                                                            {

  double sumCjkSquare;
  double vjjLess;
  double sumCikjk;

  cc[0][0] = std::sqrt(vv[0][0]);

  for(int j=1 ; j < ndim ; j++) {
    cc[j][0] = vv[j][0]/cc[0][0];
  }

  for(int j=1 ; j < ndim ; j++) {

    sumCjkSquare = 0.0;
    for (int k=0 ; k < j ; k++) sumCjkSquare += cc[j][k]*cc[j][k];

    vjjLess =  vv[j][j] - sumCjkSquare;

    //check for the case that vjjLess is negative
    cc[j][j] = std::sqrt(std::fabs(vjjLess));

    for (int i=j+1 ; i < ndim ; i++) {
      sumCikjk = 0.;
      for(int k=0 ; k < j ; k++) sumCikjk += cc[i][k]*cc[j][k];
      cc[i][j] = (vv[i][j] - sumCikjk)/cc[j][j];
    }
  }
}

double GflashHadronShowerProfile::fLnE1	(	double	einc,
		const double *	par
	)		`[protected]`

Definition at line 491 of file GflashHadronShowerProfile.cc.

References create_public_lumi_plots::log.

                                                                      {
  double func = 0.0;
  if(einc>0.0) func = par[0]+par[1]*std::log(einc);
  return func;
}

double GflashHadronShowerProfile::fTanh	(	double	einc,
		const double *	par
	)		`[protected]`

Definition at line 485 of file GflashHadronShowerProfile.cc.

References create_public_lumi_plots::log.

Referenced by hadronicParameterization(), GflashAntiProtonShowerProfile::loadParameters(), GflashProtonShowerProfile::loadParameters(), GflashKaonPlusShowerProfile::loadParameters(), GflashPiKShowerProfile::loadParameters(), and GflashKaonMinusShowerProfile::loadParameters().

                                                                      {
  double func = 0.0;
  if(einc>0.0) func = par[0]+par[1]*std::tanh(par[2]*(std::log(einc)-par[3])) + par[4]*std::log(einc);
  return func;
}

double GflashHadronShowerProfile::gammaProfile	(	double	alpha,
		double	beta,
		double	depth,
		double	lengthUnit
	)		`[protected]`

Definition at line 503 of file GflashHadronShowerProfile.cc.

References create_public_lumi_plots::exp, funct::pow(), and x.

Referenced by twoGammaProfile().

                                                                                                               {
  double gamma = 0.0;
  //  if(alpha > 0 && beta > 0 && lengthUnit > 0) {
  if(showerDepth>0.0) {
    Genfun::LogGamma lgam;
    double x = showerDepth*(beta/lengthUnit);
    gamma = (beta/lengthUnit)*std::pow(x,alpha-1.0)*std::exp(-x)/std::exp(lgam(alpha));
  }
  return gamma;
}

void GflashHadronShowerProfile::getFluctuationVector	(	double *	lowTriangle,
		double *	correlationVector
	)		`[protected]`

Definition at line 370 of file GflashHadronShowerProfile.cc.

References doCholeskyReduction(), i, j, and Gflash::NPar.

Referenced by GflashAntiProtonShowerProfile::loadParameters(), GflashProtonShowerProfile::loadParameters(), GflashKaonPlusShowerProfile::loadParameters(), GflashPiKShowerProfile::loadParameters(), and GflashKaonMinusShowerProfile::loadParameters().

                                                                                                   {

    const int dim = Gflash::NPar;

    double **xr   = new double *[dim];
    double **xrho = new double *[dim];
    
    for(int j=0;j<dim;j++) {
      xr[j]   = new double [dim];
      xrho[j] = new double [dim];
    }
    
    for(int i = 0; i < dim; i++) {
      for(int j = 0; j < i+1 ; j++) {
        if(j==i) xrho[i][j] = 1.0;
        else {
          xrho[i][j] = lowTriangle[i*(i-1)/2 + j];
          xrho[j][i] = xrho[i][j];
        }
      }
    }

    doCholeskyReduction(xrho,xr,dim);

    for(int i = 0 ; i < dim ; i++) {
      for (int j = 0 ; j < i+1 ; j++){
        correlationVector[i*(i+1)/2 + j] = xr[i][j];
      }
    }

    for(int j=0;j<dim;j++) delete [] xr[j];
    delete [] xr;
    for(int j=0;j<dim;j++) delete [] xrho[j];
    delete [] xrho;
}

std::vector<GflashHit>& GflashHadronShowerProfile::getGflashHitList ( ) [inline]

Definition at line 30 of file GflashHadronShowerProfile.h.

References theGflashHitList.

Referenced by CalorimetryManager::HDShowerSimulation(), and GflashHadronShowerModel::makeHits().

{return theGflashHitList;};

GflashShowino* GflashHadronShowerProfile::getGflashShowino ( ) [inline]

Definition at line 29 of file GflashHadronShowerProfile.h.

References theShowino.

Referenced by CalorimetryManager::HDShowerSimulation().

{ return theShowino; }

int GflashHadronShowerProfile::getNumberOfSpots ( Gflash::CalorimeterNumber kCalor ) [protected]

Definition at line 436 of file GflashHadronShowerProfile.cc.

References GflashShowino::getEnergy(), GflashShowino::getShowerType(), Gflash::kENCA, Gflash::kESPM, Gflash::kHB, Gflash::kHE, create_public_lumi_plots::log, max(), and theShowino.

Referenced by hadronicParameterization().

                                                                              {
  //generator number of spots: energy dependent Gamma distribution of Nspots based on Geant4
  //replacing old parameterization of H1,
  //int numberOfSpots = std::max( 50, static_cast<int>(80.*std::log(einc)+50.));

  double einc = theShowino->getEnergy();
  int showerType = theShowino->getShowerType();

  int numberOfSpots = 0;
  double nmean  = 0.0;
  double nsigma = 0.0;

  if(showerType == 0 || showerType == 1 || showerType == 4 || showerType == 5 ) {
    if(kCalor == Gflash::kESPM || kCalor == Gflash::kENCA) {
      nmean = 10000 + 5000*log(einc);
      nsigma = 1000;
    }
    if(kCalor == Gflash::kHB || kCalor == Gflash::kHE) {
      nmean =  5000 + 2500*log(einc);
      nsigma =  500;
    }
  }
  else if (showerType == 2 || showerType == 3 || showerType == 6 || showerType == 7 ) {
    if(kCalor == Gflash::kHB || kCalor == Gflash::kHE) {
      nmean =  5000 + 2500*log(einc);
      nsigma =  500;
    }
    else {
      nmean = 10000;
      nsigma = 1000;
    }
  }
  //@@@need correlation and individual fluctuation on alphaNspots and betaNspots here:
  //evaluating covariance should be straight forward since the distribution is 'one' Gamma

  numberOfSpots = std::max(500,static_cast<int> (nmean+nsigma* CLHEP::RandGaussQ::shoot()));

  //until we optimize the reduction scale in the number of Nspots
      
  if( kCalor == Gflash::kESPM ||  kCalor == Gflash::kENCA) {
    numberOfSpots = static_cast<int>(numberOfSpots/100);
  }
  else {
    numberOfSpots = static_cast<int>(numberOfSpots/3.0);
  }

  return numberOfSpots;
}

void GflashHadronShowerProfile::hadronicParameterization ( )

Definition at line 40 of file GflashHadronShowerProfile.cc.

References GflashShowino::addEnergyDeposited(), Gflash::divisionStep, alignCSCRings::e, energyScale, Gflash::EtaMax, fTanh(), Gflash::getCalorimeterNumber(), GflashShowino::getEnergy(), GflashShowino::getEnergyDeposited(), GflashTrajectory::getGflashTrajectoryPoint(), GflashShowino::getGlobalTime(), GflashShowino::getHelix(), getNumberOfSpots(), GflashShowino::getPathLength(), GflashTrajectory::getPathLengthAtRhoEquals(), GflashShowino::getPathLengthAtShower(), GflashShowino::getPosition(), GflashShowino::getPositionAtShower(), GflashShowino::getStepLengthToOut(), Gflash::ho_nonzero, hoProfile(), Gflash::kHB, Gflash::kHE, Gflash::kHO, Gflash::kNULL, locateHitPosition(), create_public_lumi_plots::log, longitudinalProfile(), max(), medianLateralArm(), Gflash::MinEnergyCutOffForHO, Gflash::Rmax, Gflash::Rmin, GflashHit::setEnergy(), GflashShowino::setPathLength(), GflashHit::setPosition(), GflashHit::setTime(), theGflashHcalOuter, theGflashHitList, theShowino, GflashShowino::updateShowino(), and Gflash::Zmax.

Referenced by GflashHadronShowerModel::DoIt(), and CalorimetryManager::HDShowerSimulation().

{
  // The skeleton of this method is based on the fortran code gfshow.F originally written  
  // by S. Peters and G. Grindhammer (also see NIM A290 (1990) 469-488), but longitudinal
  // parameterizations of hadron showers are significantly modified for the CMS calorimeter  

  // unit convention: energy in [GeV] and length in [cm]
  // intrinsic properties of hadronic showers (lateral shower profile)
  
  // The step size of showino along the helix trajectory in cm unit
  double showerDepthR50 = 0.0;
  bool firstHcalHit = true;

  //initial valuses that will be changed as the shower developes
  double stepLengthLeft = theShowino->getStepLengthToOut();

  double deltaStep = 0.0;
  double showerDepth = 0.0;

  Gflash::CalorimeterNumber whichCalor = Gflash::kNULL;

  theGflashHitList.clear();

  GflashHit aHit;

  while(stepLengthLeft > 0.0) {

    // update shower depth and stepLengthLeft
    if ( stepLengthLeft < Gflash::divisionStep ) {
      deltaStep = stepLengthLeft;
      stepLengthLeft  = 0.0;
    }
    else {
      deltaStep = Gflash::divisionStep;
      stepLengthLeft -= deltaStep;
    }

    showerDepth += deltaStep;
    showerDepthR50 += deltaStep;

    //update GflashShowino
    theShowino->updateShowino(deltaStep);

    //evaluate energy in this deltaStep along the longitudinal shower profile
    double heightProfile = 0.; 
    double deltaEnergy = 0.;

    whichCalor = Gflash::getCalorimeterNumber(theShowino->getPosition());

    //skip if Showino is outside envelopes
    if(whichCalor == Gflash::kNULL ) continue;

    heightProfile = longitudinalProfile();

    //skip if the delta energy for this step will be very small
    if(heightProfile < 1.00e-08 ) continue;

    //get energy deposition for this step 
    deltaEnergy =  heightProfile*Gflash::divisionStep*energyScale[whichCalor];    
    theShowino->addEnergyDeposited(deltaEnergy);

    //apply sampling fluctuation if showino is inside the sampling calorimeter
    double hadronicFraction = 1.0;
    double fluctuatedEnergy = deltaEnergy;

    int nSpotsInStep = std::max(1,static_cast<int>(getNumberOfSpots(whichCalor)*(deltaEnergy/energyScale[whichCalor]) ));
    double sampleSpotEnergy = hadronicFraction*fluctuatedEnergy/nSpotsInStep;

    // Lateral shape and fluctuations

    if((whichCalor==Gflash::kHB || whichCalor==Gflash::kHE) && firstHcalHit) {
      firstHcalHit = false;
      //reset the showerDepth used in the lateral parameterization inside Hcal
      showerDepthR50 = Gflash::divisionStep;
    }

    //evaluate the fluctuated median of the lateral distribution, R50
    double R50 = medianLateralArm(showerDepthR50,whichCalor);

    double hitEnergy = sampleSpotEnergy*CLHEP::GeV;
    double hitTime = theShowino->getGlobalTime()*CLHEP::nanosecond;

    Gflash::CalorimeterNumber hitCalor = Gflash::kNULL;

    for (int ispot = 0 ;  ispot < nSpotsInStep ; ispot++) {
        
      // Compute global position of generated spots with taking into account magnetic field
      // Divide deltaStep into nSpotsInStep and give a spot a global position
      double incrementPath = theShowino->getPathLength() 
        + (deltaStep/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
        
      // trajectoryPoint give a spot an imaginary point along the shower development
      GflashTrajectoryPoint trajectoryPoint;
      theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint,incrementPath);
        
      Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint,R50);   
        
      hitCalor = Gflash::getCalorimeterNumber(hitPosition);       

      if( hitCalor == Gflash::kNULL) continue;

      hitPosition *= CLHEP::cm;

      aHit.setTime(hitTime);
      aHit.setEnergy(hitEnergy);
      aHit.setPosition(hitPosition);
      theGflashHitList.push_back(aHit);
        
    } // end of for spot iteration
  } // end of while for longitudinal integration

  //HO parameterization

  if(theShowino->getEnergy()> Gflash::MinEnergyCutOffForHO &&
     fabs(theShowino->getPositionAtShower().pseudoRapidity()) < Gflash::EtaMax[Gflash::kHO] && 
     theGflashHcalOuter) {
    
    //non zero ho fraction to simulate based on geant4
    double nonzeroProb = 0.7*fTanh(theShowino->getEnergy(),Gflash::ho_nonzero);
    double r0 = CLHEP::HepUniformRand();
    double leftoverE =  theShowino->getEnergy() - theShowino->getEnergyDeposited();

    //@@@ nonzeroProb is not random - need further correlation for non-zero HO energy

    if( r0 < nonzeroProb && leftoverE > 0.0) {
            
      //starting path Length and stepLengthLeft
      double pathLength = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHO]-10);
      stepLengthLeft = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHO]+10) - pathLength;
      showerDepth = pathLength - theShowino->getPathLengthAtShower();
      
      theShowino->setPathLength(pathLength);
      
      double pathLengthx = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHO]);
      double pathLengthy = theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
      
      while(stepLengthLeft > 0.0) {
        
        // update shower depth and stepLengthLeft
        if ( stepLengthLeft < Gflash::divisionStep ) {
          deltaStep = stepLengthLeft;
          stepLengthLeft  = 0.0;
        }
        else {
          deltaStep = Gflash::divisionStep;
          stepLengthLeft -= deltaStep;
        }
        
        showerDepth += deltaStep;
        
        //update GflashShowino
        theShowino->updateShowino(deltaStep);
        
        //evaluate energy in this deltaStep along the longitudinal shower profile
        double heightProfile = 0.; 
        double deltaEnergy = 0.;
        
        double hoScale  = leftoverE*(pathLengthx-pathLengthy)/(pathLengthx- theShowino->getPathLengthAtShower());
        double refDepth = theShowino->getPathLength() 
          - theShowino->getHelix()->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]);
        
        if( refDepth > 0) {
          heightProfile = hoProfile(theShowino->getPathLength(),refDepth);
          deltaEnergy = heightProfile*Gflash::divisionStep*hoScale;
        }
        
        int nSpotsInStep = std::max(50,static_cast<int>((160.+40* CLHEP::RandGaussQ::shoot())*std::log(theShowino->getEnergy())+50.));
        
        double hoFraction = 1.00;
        double poissonProb = CLHEP::RandPoissonQ::shoot(1.0);
        
        double fluctuatedEnergy = deltaEnergy*poissonProb;
        double sampleSpotEnergy = hoFraction*fluctuatedEnergy/nSpotsInStep;
        
        // Lateral shape and fluctuations
        
        //evaluate the fluctuated median of the lateral distribution, R50
        double R50 = medianLateralArm(showerDepth,Gflash::kHB);

        double hitEnergy = sampleSpotEnergy*CLHEP::GeV;
        double hitTime = theShowino->getGlobalTime()*CLHEP::nanosecond;
        
        for (int ispot = 0 ;  ispot < nSpotsInStep ; ispot++) {
          
          double incrementPath = theShowino->getPathLength() 
            + (deltaStep/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
          
          // trajectoryPoint give a spot an imaginary point along the shower development
          GflashTrajectoryPoint trajectoryPoint;
          theShowino->getHelix()->getGflashTrajectoryPoint(trajectoryPoint,incrementPath);
          
          Gflash3Vector hitPosition = locateHitPosition(trajectoryPoint,R50);   
          hitPosition *= CLHEP::cm;

          if(std::fabs(hitPosition.getZ()/CLHEP::cm) > Gflash::Zmax[Gflash::kHO] ) continue;

          aHit.setTime(hitTime);
          aHit.setEnergy(hitEnergy);
          aHit.setPosition(hitPosition);
          theGflashHitList.push_back(aHit);
          
        } // end of for HO spot iteration
      } // end of while for HO longitudinal integration
    }
  }

  //  delete theGflashNavigator;

}

double GflashHadronShowerProfile::hoProfile	(	double	pathLength,
		double	refDepth
	)		`[protected]`

Definition at line 357 of file GflashHadronShowerProfile.cc.

References create_public_lumi_plots::exp, GflashTrajectory::getGflashTrajectoryPoint(), GflashShowino::getHelix(), GflashTrajectoryPoint::getPosition(), Gflash::intLength, Gflash::kHO, funct::sin(), and theShowino.

Referenced by hadronicParameterization().

                                                                              {

  double heightProfile = 0;

  GflashTrajectoryPoint tempPoint;
  theShowino->getHelix()->getGflashTrajectoryPoint(tempPoint,pathLength);

  double dint = 1.4*Gflash::intLength[Gflash::kHO]*std::sin(tempPoint.getPosition().getTheta());
  heightProfile = std::exp(-1.0*refDepth/dint);

  return heightProfile;
}

void GflashHadronShowerProfile::initialize	(	int	showerType,
		double	energy,
		double	globalTime,
		double	charge,
		Gflash3Vector &	position,
		Gflash3Vector &	momentum
	)

Definition at line 31 of file GflashHadronShowerProfile.cc.

References GflashShowino::initialize(), theBField, and theShowino.

Referenced by GflashHadronShowerModel::DoIt(), and CalorimetryManager::HDShowerSimulation().

                                                                                            { 

  //initialize GflashShowino for this track
  theShowino->initialize(showerType, energy, globalTime, charge, 
                         position, momentum, theBField);

}

void GflashHadronShowerProfile::loadParameters ( ) [virtual]

Reimplemented in GflashAntiProtonShowerProfile, GflashKaonMinusShowerProfile, GflashKaonPlusShowerProfile, GflashPiKShowerProfile, and GflashProtonShowerProfile.

Definition at line 250 of file GflashHadronShowerProfile.cc.

Referenced by GflashHadronShowerModel::DoIt(), and CalorimetryManager::HDShowerSimulation().

{
  edm::LogInfo("SimGeneralGFlash") << "GflashHadronShowerProfile::loadParameters() "
                                   << "should be implimented for each particle type";
}

Gflash3Vector GflashHadronShowerProfile::locateHitPosition	(	GflashTrajectoryPoint &	point,
		double	lateralArm
	)		`[protected]`

Definition at line 278 of file GflashHadronShowerProfile.cc.

References funct::cos(), GflashTrajectoryPoint::getCrossUnitVector(), GflashTrajectoryPoint::getOrthogonalUnitVector(), GflashTrajectoryPoint::getPosition(), GflashHistogram::getStoreFlag(), GflashHistogram::lateralx, GflashHistogram::lateraly, Gflash::maxLateralArmforR50, min, position, GflashHistogram::rshower, funct::sin(), mathSSE::sqrt(), and theHisto.

Referenced by hadronicParameterization().

{
  // Smearing in r according to f(r)= 2.*r*R50**2/(r**2+R50**2)**2
  double rnunif = CLHEP::HepUniformRand();
  double rxPDF = std::sqrt(rnunif/(1.-rnunif));
  double rShower = lateralArm*rxPDF;
  
  //rShower within maxLateralArmforR50
  rShower = std::min(Gflash::maxLateralArmforR50,rShower);

  // Uniform smearing in phi
  double azimuthalAngle = CLHEP::twopi*CLHEP::HepUniformRand(); 

  // actual spot position by adding a radial vector to a trajectoryPoint
  Gflash3Vector position = point.getPosition() +
    rShower*std::cos(azimuthalAngle)*point.getOrthogonalUnitVector() +
    rShower*std::sin(azimuthalAngle)*point.getCrossUnitVector();
  
  //@@@debugging histograms
  if(theHisto->getStoreFlag()) {
    theHisto->rshower->Fill(rShower);
    theHisto->lateralx->Fill(rShower*std::cos(azimuthalAngle));
    theHisto->lateraly->Fill(rShower*std::sin(azimuthalAngle));
  }
  return position;
}

double GflashHadronShowerProfile::longitudinalProfile ( ) [protected]

Definition at line 306 of file GflashHadronShowerProfile.cc.

References Gflash::getCalorimeterNumber(), GflashShowino::getDepth(), GflashShowino::getPathLengthAtShower(), GflashShowino::getPathLengthOnEcal(), GflashShowino::getPosition(), GflashShowino::getShowerType(), GflashShowino::getStepLengthToHcal(), Gflash::kENCA, Gflash::kESPM, Gflash::kHB, Gflash::kHE, longEcal, longHcal, pos, theShowino, and twoGammaProfile().

Referenced by hadronicParameterization().

                                                      {

  double heightProfile = 0.0;

  Gflash3Vector pos = theShowino->getPosition();
  int showerType = theShowino->getShowerType();
  double showerDepth = theShowino->getDepth();
  double transDepth = theShowino->getStepLengthToHcal();

  // Energy in a delta step (dz) = (energy to deposite)*[Gamma(z+dz)-Gamma(z)]*dz
  // where the incomplete Gamma function gives an intergrate probability of the longitudinal 
  // shower up to the shower depth (z).
  // Instead, we use approximated energy; energy in dz = (energy to deposite)*gamma(z)*dz
  // where gamma is the Gamma-distributed probability function

  Gflash::CalorimeterNumber whichCalor = Gflash::getCalorimeterNumber(pos);

  if(showerType == 0 || showerType == 4 ) {
    double shiftDepth = theShowino->getPathLengthOnEcal() - theShowino->getPathLengthAtShower();
    if(shiftDepth > 0 ) {
      heightProfile = twoGammaProfile(longEcal,showerDepth-shiftDepth,whichCalor);
    }
    else  {
      heightProfile = 0.;
      //      std::cout << "negative shiftDepth for showerType 0 " << shiftDepth << std::endl;
    }
  }  
  else if(showerType == 1 || showerType == 5 ) {
    if(whichCalor == Gflash::kESPM || whichCalor == Gflash::kENCA ) {
      heightProfile = twoGammaProfile(longEcal,showerDepth,whichCalor);
    }
    else if(whichCalor == Gflash::kHB || whichCalor == Gflash::kHE) {
      heightProfile = twoGammaProfile(longHcal,showerDepth-transDepth,whichCalor);
    }
    else  heightProfile = 0.;
  }  
  else if (showerType == 2 || showerType == 6 ) {
    //two gammas between crystal and Hcal
    if((showerDepth - transDepth) > 0.0) {
      heightProfile = twoGammaProfile(longHcal,showerDepth-transDepth,Gflash::kHB);
    }
    else heightProfile = 0.;
  }
  else if (showerType == 3 || showerType == 7 ) {
    //two gammas inside Hcal
    heightProfile = twoGammaProfile(longHcal,showerDepth,Gflash::kHB);
  }

  return heightProfile;
}

double GflashHadronShowerProfile::medianLateralArm	(	double	depth,
		Gflash::CalorimeterNumber	kCalor
	)		`[protected]`

Definition at line 256 of file GflashHadronShowerProfile.cc.

References create_public_lumi_plots::exp, Gflash::kNULL, lateralPar, create_public_lumi_plots::log, max(), Gflash::maxShowerDepthforR50, min, funct::pow(), and mathSSE::sqrt().

Referenced by hadronicParameterization().

{
  double lateralArm = 0.0;
  if(kCalor != Gflash::kNULL) {
    
    double showerDepthR50X = std::min(showerDepthR50/22.4, Gflash::maxShowerDepthforR50);
    double R50          = lateralPar[kCalor][0] + std::max(0.0,lateralPar[kCalor][1]) * showerDepthR50X;
    double varinanceR50 = std::pow((lateralPar[kCalor][2] + lateralPar[kCalor][3] * showerDepthR50X) * R50, 2);
    
    // Simulation of lognormal distribution
    
    if(R50>0) {
      double sigmaSq  = std::log(varinanceR50/(R50*R50)+1.0);
      double sigmaR50 = std::sqrt(sigmaSq);
      double meanR50  = std::log(R50) - (sigmaSq/2.);
      
      lateralArm = std::exp(meanR50 + sigmaR50* CLHEP::RandGaussQ::shoot());
    }
  }
  return lateralArm;
}

void GflashHadronShowerProfile::setEnergyScale	(	double	einc,
		Gflash3Vector	ssp
	)		`[protected]`

double GflashHadronShowerProfile::twoGammaProfile	(	double *	par,
		double	depth,
		Gflash::CalorimeterNumber	kIndex
	)		`[protected]`

Definition at line 514 of file GflashHadronShowerProfile.cc.

References create_public_lumi_plots::exp, gammaProfile(), Gflash::intLength, max(), min, and Gflash::radLength.

Referenced by longitudinalProfile().

                                                                                                               {
  double twoGamma = 0.0;

  longPar[0] = std::min(1.0,longPar[0]);
  longPar[0] = std::max(0.0,longPar[0]);

  if(longPar[3] > 4.0 || longPar[4] > 4.0) {
    double rfactor = 2.0/std::max(longPar[3],longPar[4]);
    longPar[3] = rfactor*(longPar[3]+1.0);  
    longPar[4] *= rfactor;  
  }

  twoGamma  = longPar[0]*gammaProfile(exp(longPar[1]),exp(longPar[2]),depth,Gflash::radLength[kIndex])
         +(1-longPar[0])*gammaProfile(exp(longPar[3]),exp(longPar[4]),depth,Gflash::intLength[kIndex]);
  return twoGamma;
}