|
|
|
#include <SimG4Core/GFlash/interface/GflashHadronShowerProfile.h>
Definition at line 18 of file GflashHadronShowerProfile.h.
| GflashHadronShowerProfile::GflashHadronShowerProfile | ( | G4Region * | envelope | ) |
Definition at line 17 of file GflashHadronShowerProfile.cc.
References Exception, fillFluctuationVector(), GflashHistogram::instance(), edm::Service< T >::isAvailable(), jCalorimeter, Gflash::kNULL, showerType, theHelix, theHisto, theRandGamma, and theRandGauss.
00018 { 00019 showerType = 0; 00020 jCalorimeter = Gflash::kNULL; 00021 theHelix = new GflashTrajectory; 00022 theHisto = GflashHistogram::instance(); 00023 00024 edm::Service<edm::RandomNumberGenerator> rng; 00025 if ( ! rng.isAvailable()) { 00026 throw cms::Exception("Configuration") 00027 << "GflashHadronShowerProfile requires RandomNumberGeneratorService\n" 00028 << "which is not present in the configuration file. " 00029 << "You must add the service\n in the configuration file or " 00030 << "remove the modules that require it."; 00031 } 00032 theRandGauss = new CLHEP::RandGaussQ(rng->getEngine()); 00033 theRandGamma = new CLHEP::RandGamma(rng->getEngine()); 00034 00035 //correllation and fluctuation matrix 00036 fillFluctuationVector(); 00037 }
| GflashHadronShowerProfile::~GflashHadronShowerProfile | ( | ) |
Definition at line 39 of file GflashHadronShowerProfile.cc.
References theHelix, theRandGamma, and theRandGauss.
00040 { 00041 // delete theGflashStep; 00042 delete theHelix; 00043 delete theRandGauss; 00044 delete theRandGamma; 00045 }
| void GflashHadronShowerProfile::clearSpotList | ( | ) | [inline] |
Definition at line 28 of file GflashHadronShowerProfile.h.
References aEnergySpotList.
Referenced by GflashHadronShowerModel::DoIt().
00028 { aEnergySpotList.clear(); }
| void GflashHadronShowerProfile::doCholeskyReduction | ( | G4double ** | cc, | |
| G4double ** | vv, | |||
| const G4int | ndim | |||
| ) | [private] |
Referenced by fillFluctuationVector().
| void GflashHadronShowerProfile::fillFluctuationVector | ( | ) | [private] |
Definition at line 547 of file GflashHadronShowerProfile.cc.
References correlationVector, GenMuonPlsPt100GeV_cfg::cout, doCholeskyReduction(), lat::endl(), i, j, k, Gflash::NDim, Gflash::NRegion, Gflash::NStart, and Gflash::rho.
Referenced by GflashHadronShowerProfile().
00547 { 00548 // G4double RMX[186]; //21*6 = 186 00549 00550 for(G4int k = 0 ; k < Gflash::NRegion ; k++) { 00551 const G4int dim = Gflash::NDim[k]; 00552 G4double **xr = new G4double *[dim]; 00553 G4double **xrho = new G4double *[dim]; 00554 00555 for(G4int j=0;j<dim;j++) { 00556 xr[j] = new G4double [dim]; 00557 xrho[j] = new G4double [dim]; 00558 } 00559 00560 for(G4int i = 0; i < dim; i++) { 00561 for(G4int j = 0; j < i+1 ; j++) { 00562 xrho[i][j] = Gflash::rho[i+Gflash::NRegion*k][j]; 00563 xrho[i][j] = Gflash::rho[i][j]; 00564 xrho[j][i] = xrho[i][j]; 00565 } 00566 } 00567 00568 doCholeskyReduction(xrho,xr,dim); 00569 00570 for(G4int i = 0 ; i < dim ; i++) { 00571 for (G4int j = 0 ; j < i+1 ; j++){ 00572 correlationVector[Gflash::NStart[k]+i*(i+1)/2 + j] = xr[i][j]; 00573 } 00574 } 00575 00576 std::cout << "this should be calcuated at constructor" << std::endl; 00577 for(int i = 0; i < 21 ; i++) std::cout << correlationVector[i] << std::endl; 00578 00579 for(G4int j=0;j<dim;j++) delete [] xr[j]; 00580 delete [] xr; 00581 for(G4int j=0;j<dim;j++) delete [] xrho[j]; 00582 delete [] xrho; 00583 } 00584 }
| Gflash::CalorimeterNumber GflashHadronShowerProfile::getCalorimeterNumber | ( | ) | [inline, private] |
Definition at line 36 of file GflashHadronShowerProfile.h.
References jCalorimeter.
Referenced by hadronicParameterization().
00036 {return jCalorimeter;}
| Gflash::CalorimeterNumber GflashHadronShowerProfile::getCalorimeterNumber | ( | const G4ThreeVector | position | ) |
Definition at line 47 of file GflashHadronShowerProfile.cc.
References eta, Gflash::EtaMax, Gflash::EtaMin, index, Gflash::kENCA, Gflash::kESPM, Gflash::kHB, Gflash::kHE, Gflash::kNULL, Gflash::Rmax, Gflash::Rmin, Gflash::Zmax, and Gflash::Zmin.
00048 { 00049 Gflash::CalorimeterNumber index = Gflash::kNULL; 00050 G4double eta = position.getEta(); 00051 00052 //central 00053 if (fabs(eta) < Gflash::EtaMax[Gflash::kESPM] || fabs(eta) < Gflash::EtaMax[Gflash::kHB]) { 00054 if(position.getRho() > Gflash::Rmin[Gflash::kESPM] && 00055 position.getRho() < Gflash::Rmax[Gflash::kESPM] ) { 00056 index = Gflash::kESPM; 00057 } 00058 if(position.getRho() > Gflash::Rmin[Gflash::kHB] && 00059 position.getRho() < Gflash::Rmax[Gflash::kHB]) { 00060 index = Gflash::kHB; 00061 } 00062 } 00063 //forward 00064 else if (fabs(eta) > Gflash::EtaMin[Gflash::kENCA] || fabs(eta) > Gflash::EtaMin[Gflash::kHE]) { 00065 if( fabs(position.getZ()) > Gflash::Zmin[Gflash::kENCA] && 00066 fabs(position.getZ()) < Gflash::Zmax[Gflash::kENCA] ) { 00067 index = Gflash::kENCA; 00068 } 00069 if( fabs(position.getZ()) > Gflash::Zmin[Gflash::kHE] && 00070 fabs(position.getZ()) < Gflash::Zmax[Gflash::kHE] ) { 00071 index = Gflash::kHE; 00072 } 00073 } 00074 return index; 00075 }
| std::vector<GflashEnergySpot>& GflashHadronShowerProfile::getEnergySpotList | ( | ) | [inline] |
Definition at line 30 of file GflashHadronShowerProfile.h.
References aEnergySpotList.
Referenced by GflashHadronShowerModel::DoIt().
00030 {return aEnergySpotList;};
| void GflashHadronShowerProfile::hadronicParameterization | ( | const G4FastTrack & | fastTrack | ) |
Definition at line 77 of file GflashHadronShowerProfile.cc.
References aEnergySpotList, funct::cos(), energyToDeposit, funct::exp(), getCalorimeterNumber(), GflashTrajectoryPoint::getCrossUnitVector(), GflashEnergySpot::getEnergy(), GflashTrajectory::getGflashTrajectoryPoint(), GflashTrajectoryPoint::getOrthogonalUnitVector(), GflashTrajectory::getPathLengthAtRhoEquals(), GflashTrajectory::getPathLengthAtZ(), GflashTrajectoryPoint::getPosition(), GflashHistogram::getStoreFlag(), GflashHistogram::gfhlongProfile, GflashTrajectory::initializeTrajectory(), jCalorimeter, Gflash::kENCA, Gflash::kESPM, Gflash::kHB, Gflash::kHE, lateralPar, GflashHistogram::lateralx, GflashHistogram::lateraly, loadParameters(), funct::log(), longitudinalProfile(), max, min, funct::pow(), Gflash::Rmax, Gflash::Rmin, rMoliere, GflashHistogram::rshower, GflashEnergySpot::setEnergy(), GflashEnergySpot::setPosition(), showerType, funct::sin(), funct::sqrt(), theHelix, theHisto, theRandGauss, Gflash::Zmax, and Gflash::Zmin.
Referenced by GflashHadronShowerModel::DoIt().
00078 { 00079 // The skeleton of this method is based on the fortran code gfshow.F originally written 00080 // by S. Peters and G. Grindhammer (also see NIM A290 (1990) 469-488), but longitudinal 00081 // parameterizations of hadron showers are significantly modified for the CMS calorimeter 00082 00083 // unit convention: energy in [GeV] and length in [cm] 00084 00085 // maximum number of energy spots 00086 const G4int maxNumberOfSpots = 10000; 00087 00088 // low energy cutoff (unit in GeV) 00089 // const G4double energyCutoff = 0.01; 00090 00091 // intrinsic properties of hadronic showers (lateral shower profile) 00092 const G4double maxShowerDepthforR50 = 10.0; 00093 00094 G4double rShower = 0.; 00095 G4double rGauss = theRandGauss->fire(); 00096 00097 // The shower starting point is the PostStepPoint of Hadronic Inelestic interaction; 00098 // see GflashHadronShowerModel::ModelTrigger 00099 00100 G4double einc = fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV; 00101 G4ThreeVector positionShower = fastTrack.GetPrimaryTrack()->GetPosition()/cm; 00102 G4ThreeVector momentumShower = fastTrack.GetPrimaryTrack()->GetMomentum()/GeV; 00103 00104 //find the calorimeter at the shower starting point 00105 jCalorimeter = getCalorimeterNumber(positionShower); 00106 00107 //get all necessary parameters for hadronic shower profiles including energyToDeposit 00108 loadParameters(fastTrack); 00109 00110 // The direction of shower is assumed to be along the showino trajectory 00111 // inside the magnetic field; 00112 00113 const G4double bField = 4.0*tesla; 00114 00115 double charge = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge(); 00116 00117 theHelix->initializeTrajectory(momentumShower,positionShower,charge,bField/tesla); 00118 00119 //path Length from the origin to the shower starting point in cm 00120 00121 G4double pathLength0 = 0; 00122 G4double transDepth = 0; 00123 00124 // The step length left is the total path length from the shower starting point to 00125 // the maximum distance inside paramerized envelopes 00126 00127 //distance to the end of HB/HB now 00128 //@@@extend the trajectory outside bField and HO later 00129 G4double stepLengthLeft = 0.0; 00130 00131 if(jCalorimeter == Gflash::kESPM || jCalorimeter == Gflash::kHB ) { 00132 pathLength0 = theHelix->getPathLengthAtRhoEquals(positionShower.getRho()); 00133 stepLengthLeft = theHelix->getPathLengthAtRhoEquals(Gflash::Rmax[Gflash::kHB]) 00134 - theHelix->getPathLengthAtRhoEquals(positionShower.getRho()); 00135 if(showerType == 3 ) { 00136 transDepth = theHelix->getPathLengthAtRhoEquals(Gflash::Rmin[Gflash::kHB]) - pathLength0; 00137 } 00138 } 00139 else if (jCalorimeter == Gflash::kENCA || jCalorimeter == Gflash::kHE ) { 00140 pathLength0 = theHelix->getPathLengthAtZ(positionShower.getZ()); 00141 stepLengthLeft = theHelix->getPathLengthAtRhoEquals(Gflash::Zmax[Gflash::kHE]) 00142 - theHelix->getPathLengthAtRhoEquals(positionShower.getZ()); 00143 if ( showerType ==7 ) { 00144 transDepth = theHelix->getPathLengthAtZ(Gflash::Zmin[Gflash::kHE]) - pathLength0; 00145 } 00146 } 00147 else { 00148 //@@@extend for HF later 00149 stepLengthLeft = 200.0; 00150 } 00151 00152 G4double pathLength = pathLength0; // this will grow along the shower development 00153 00154 // Limit number of spots to maxNumberOfSpots 00155 G4int numberOfSpots = std::max( 50, static_cast<int>(800.*std::log(einc)+50.)); 00156 numberOfSpots = std::min(numberOfSpots,maxNumberOfSpots); 00157 00158 // Spot energy to simulate sampling fluctuations (SampleEnergySpot) and 00159 // number of spots needed to fullfill geometry constraints(TotalNumberOfSpots): 00160 00161 G4double spotEnergy = energyToDeposit/numberOfSpots; 00162 00163 // The step size of showino along the helix trajectory in cm unit 00164 const G4double divisionStep = 1.0; 00165 G4double deltaStep = 0.0; 00166 G4double showerDepth = 0.0; 00167 00168 00169 G4int totalNumberOfSpots = 0; 00170 00171 //empty energy spot vector for a new track 00172 aEnergySpotList.clear(); 00173 00174 double scaleLateral = 0.0; 00175 const double rMoliere = 2.19; //Moliere Radius in [cm] 00176 00177 while(stepLengthLeft > 0.0) { 00178 00179 // update shower depth and stepLengthLeft 00180 if ( stepLengthLeft < divisionStep ) { 00181 deltaStep = stepLengthLeft; 00182 stepLengthLeft = 0.0; 00183 } 00184 else { 00185 deltaStep = divisionStep; 00186 stepLengthLeft -= deltaStep; 00187 } 00188 00189 showerDepth += deltaStep; 00190 pathLength += deltaStep; 00191 00192 // energy in this deltaStep along the longitudinal shower profile 00193 double deltaEnergy = 0.; 00194 00195 double heightProfile = longitudinalProfile(showerDepth,pathLength,transDepth); 00196 // deltaEnergy = longitudinalProfile(showerDepth)*divisionStep*energyToDeposit; 00197 deltaEnergy = heightProfile*divisionStep*energyToDeposit; 00198 00199 //@@@ When depthShower is inside Hcal, the sampling fluctuation for deposited 00200 // energy will be treated in SD. However we should put some scale factor 00201 // to relate the spot energy to the energy deposited in each geant4 step. 00202 00203 double hadronicFraction = 1.0; 00204 G4double fluctuatedEnergy = deltaEnergy; 00205 G4int nSpotsInStep = std::max(1,static_cast<int>(fluctuatedEnergy/spotEnergy)); 00206 G4double sampleSpotEnergy = hadronicFraction*fluctuatedEnergy/nSpotsInStep; 00207 00208 // Sampling fluctuations determine the number of spots: 00209 // if (insideSampling(positionShower)) samplingFluctuation(fluctuatedEnergy,einc); 00210 // G4double hadSpotEnergy = std::max(0.,sampleSpotEnergy * hadronicFraction); 00211 // G4double emSpotEnergy = std::max(0.,sampleSpotEnergy - hadSpotEnergy); 00212 // hadSpotEnergy *= Gflash::PBYMIP[jCalorimeter]; 00213 00214 00215 //@@@this part of code may not be not need, but leave it for further consideration 00216 00217 // Lateral shape and fluctuations 00218 00219 double showerDepthR50 = std::min(showerDepth/20.7394, maxShowerDepthforR50); 00220 00221 double R50 = lateralPar[0] + lateralPar[1] * showerDepthR50; 00222 double varinanceR50 = std::pow((lateralPar[2] + lateralPar[3] * showerDepthR50) * R50, 2); 00223 00224 // Simulation of lognormal distribution 00225 00226 double sigmaSq = std::log(varinanceR50/(R50*R50)+1.0); 00227 double sigmaR50 = std::sqrt(sigmaSq); 00228 double meanR50 = std::log(R50) - (sigmaSq/2.); 00229 00230 R50 = std::exp(rGauss*sigmaR50 + meanR50); 00231 00232 // Averaging lateral scale in terms of Moliere radius 00233 // const G4double rMoliere = Gflash::RLTHAD[jCalorimeter]; 00234 00235 //@@@this should be each spot basis 00236 if(showerType == 4 || showerType == 8) { 00237 scaleLateral = (3.5+1.0*showerDepth)*rMoliere; 00238 } 00239 else { 00240 //@@@need better division for showerDepth arosse the Hcal front face 00241 if(showerDepthR50 < 2.0 ) { 00242 scaleLateral = (5.5-0.4*std::log(einc))*rMoliere; 00243 } 00244 else { 00245 scaleLateral = ( 14-1.5*std::log(einc))*rMoliere; 00246 } 00247 } 00248 // region0 && inside Ecal: scaleLateral = (5.5-0.4*logEinc)*rMoliere; 00249 // region0 && inside Hcal: scaleLateral = (14-1.5*logEinc)*rMoliere; 00250 // region1 00251 00252 R50 *= scaleLateral; 00253 00254 GflashEnergySpot eSpot; 00255 00256 for (G4int ispot = 0 ; ispot < nSpotsInStep ; ispot++) { 00257 00258 totalNumberOfSpots++; 00259 00260 // Smearing in r according to f(r)= 2.*r*R50**2/(r**2+R50**2)**2 00261 G4double rnunif = G4UniformRand(); 00262 G4double rxPDF = std::sqrt(rnunif/(1.-rnunif)); 00263 rShower = R50 * rxPDF; 00264 00265 // Uniform smearing in phi, for 66% of lateral containm. 00266 G4double azimuthalAngle = 0.0; 00267 00268 azimuthalAngle = twopi*G4UniformRand(); 00269 00270 // Compute global position of generated spots with taking into account magnetic field 00271 // Divide deltaStep into nSpotsInStep and give a spot a global position 00272 G4double incrementPath = (deltaStep/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep); 00273 00274 // trajectoryPoint give a spot an imaginary point along the shower development 00275 GflashTrajectoryPoint trajectoryPoint; 00276 theHelix->getGflashTrajectoryPoint(trajectoryPoint,pathLength+incrementPath); 00277 00278 // actual spot position by adding a radial vector to a trajectoryPoint 00279 G4ThreeVector SpotPosition = trajectoryPoint.getPosition() + 00280 rShower*std::cos(azimuthalAngle)*trajectoryPoint.getOrthogonalUnitVector() + 00281 rShower*std::sin(azimuthalAngle)*trajectoryPoint.getCrossUnitVector(); 00282 00283 //convert unit of energy to geant4 default MeV 00284 // eSpot.setEnergy((hadSpotEnergy+emSpotEnergy)*GeV); 00285 eSpot.setEnergy(sampleSpotEnergy*GeV); 00286 eSpot.setPosition(SpotPosition*cm); 00287 aEnergySpotList.push_back(eSpot); 00288 00289 //@@@debugging histograms 00290 if(theHisto->getStoreFlag()) { 00291 theHisto->rshower->Fill(rShower); 00292 theHisto->lateralx->Fill(rShower*std::cos(azimuthalAngle)); 00293 theHisto->lateraly->Fill(rShower*std::sin(azimuthalAngle)); 00294 theHisto->gfhlongProfile->Fill(pathLength+incrementPath-pathLength0,positionShower.getRho(),eSpot.getEnergy()*GeV); 00295 } 00296 } 00297 } 00298 00299 }
| G4bool GflashHadronShowerProfile::insideSampling | ( | const G4ThreeVector | pos | ) | [private] |
Definition at line 538 of file GflashHadronShowerProfile.cc.
References jCalorimeter, Gflash::kENCA, Gflash::kESPM, Gflash::kHB, and Gflash::kHE.
00538 { 00539 G4bool issampling = false; 00540 00541 if((jCalorimeter == Gflash::kHB) || (jCalorimeter == Gflash::kHE) || 00542 ((jCalorimeter == Gflash::kESPM) && ((pos.rho()/cm - 177.5) > 0)) || 00543 ((jCalorimeter == Gflash::kENCA) && ( fabs(pos.z()/cm - 391.95) > 0 ))) issampling = true; 00544 return issampling; 00545 }
| void GflashHadronShowerProfile::loadParameters | ( | const G4FastTrack & | fastTrack | ) | [private] |
Definition at line 301 of file GflashHadronShowerProfile.cc.
References correlationVector, e, energyToDeposit, i, j, jCalorimeter, k, Gflash::kENCA, Gflash::kESPM, Gflash::kHB, Gflash::kHE, lateralPar, funct::log(), longPar, longSigma, lv, max, Gflash::NRegion, Gflash::NStart, Gflash::NxN, showerType, funct::sqrt(), and theRandGauss.
Referenced by hadronicParameterization().
00302 { 00303 // Initialization of longitudinal and lateral parameters for 00304 // hadronic showers. Simulation of the intrinsic fluctuations 00305 00306 G4double einc = fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV; 00307 00308 // type of hadron showers subject to the shower starting point (ssp) 00309 00310 // showerType = 0 : default (invalid) 00311 // showerType = 1 : ssp before EB 00312 // showerType = 2 : ssp inside EB 00313 // showerType = 3 : ssp after EB before HB 00314 // showerType = 4 : ssp inside HB 00315 // showerType = 5 : ssp before EE 00316 // showerType = 6 : ssp inside EE 00317 // showerType = 7 : ssp after EE before HE 00318 // showerType = 8 : ssp inside HE 00319 00320 G4TouchableHistory* touch = (G4TouchableHistory*)(fastTrack.GetPrimaryTrack()->GetTouchable()); 00321 G4LogicalVolume* lv = touch->GetVolume()->GetLogicalVolume(); 00322 std::size_t pos1 = lv->GetName().find("EBRY"); 00323 std::size_t pos2 = lv->GetName().find("EFRY"); 00324 00325 G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition()/cm; 00326 00327 showerType = 0; 00328 00329 G4double correctionAsDepth = 0.0; 00330 00331 //central 00332 if (jCalorimeter == Gflash::kESPM || jCalorimeter == Gflash::kHB ) { 00333 00334 G4double posRho = position.getRho(); 00335 00336 if(pos1 != std::string::npos) { 00337 showerType = 2; 00338 } 00339 else { 00340 if(jCalorimeter == Gflash::kESPM) { 00341 showerType = 3; 00342 if( posRho < 129.0 ) showerType = 1; 00343 } 00344 else showerType = 4; 00345 } 00346 00347 if ( posRho < 150.0 ) { 00348 correctionAsDepth = 0.01-2.0/((posRho-150.)*(posRho-150.) +5.4*5.4); 00349 } 00350 else { 00351 correctionAsDepth = 0.03-2.0/((posRho-150.)*(posRho-150.) +4.7*4.7); 00352 } 00353 } 00354 //forward 00355 else if (jCalorimeter == Gflash::kENCA || jCalorimeter == Gflash::kHE) { 00356 if(pos2 != std::string::npos) { 00357 showerType = 6; 00358 } 00359 else { 00360 if(jCalorimeter == Gflash::kENCA) { 00361 showerType = 7; 00362 if(fabs(position.getZ()) < 330.0 ) showerType = 5; 00363 } 00364 else showerType = 8; 00365 } 00366 //@@@need z-dependent correction on the mean energy reponse 00367 } 00368 00369 00370 // total energy to deposite 00371 //@@@ need additional parameterization by the shower starting point 00372 G4double fractionEnergy = 1.0; 00373 G4double sigmaEnergy = 0.0; 00374 00375 if( showerType == 4 || showerType == 8) { 00376 //Mip-like particle 00377 fractionEnergy = 0.7125 + 0.0812*std::tanh(0.9040*(std::log(einc) - 2.6307)); 00378 sigmaEnergy = 0.0257/std::sqrt(einc) + 0.0734; 00379 } 00380 else { 00381 fractionEnergy = 0.7125 + 0.0812*std::tanh(0.9040*(std::log(einc) - 2.6307)); 00382 sigmaEnergy = 0.0844/std::sqrt(einc) + 0.0592; 00383 } 00384 00385 energyToDeposit = fractionEnergy*(1.0+correctionAsDepth)*einc*(1.0+sigmaEnergy*theRandGauss->fire()); 00386 energyToDeposit = std::max(0.0,energyToDeposit); 00387 00388 // parameters for the longitudinal profiles 00389 //@@@check longitudinal profiles of endcaps for possible varitations 00390 //@@@need to add fluctuation and correlation for individual shower 00391 00392 longPar[0][0] = 1.41*std::max(0.0,-5.96481e-03 + 0.18231*std::tanh(0.55451*(std::log(einc)-0.458775))) ; 00393 longPar[0][1] = std::max(0.0,2.01611 + 1.77483 * std::tanh(0.75719*(std::log(einc) - 2.58172))); 00394 longPar[0][2] = std::max(0.0,0.21261 + 0.24168 * std::tanh(0.76962*(std::log(einc) - 2.11936))); 00395 longPar[0][3] = std::max(0.0,1.05577e-02 + 1.00807 * std::tanh(-6.31044e-04*(std::log(einc) - 4.60658))); 00396 longPar[0][4] = 0.87*std::max(0.0,1.19845e-01 + 6.87070e-02 * std::tanh(-8.23888e-01*(std::log(einc) - 2.90178))); 00397 longPar[0][5] = std::max(0.0,2.49694e+01 + 1.10258e+01 * std::tanh(6.16435e-01*(std::log(einc) - 3.56012))); 00398 00399 longSigma[0][0] = 0.02; 00400 longSigma[0][1] = 0.16; 00401 longSigma[0][2] = 0.02; 00402 longSigma[0][3] = 0.01; 00403 longSigma[0][4] = 0.03; 00404 longSigma[0][5] = 2.50; 00405 00406 longPar[1][0] = 0.1126; 00407 longPar[1][1] = 1.3857; 00408 longPar[1][2] = std::max(0.0,1.1353 + 0.4997*std::tanh(-0.6382*(std::log(einc) - 2.0035))); 00409 longPar[1][3] = 0.2300; 00410 longPar[1][4] = 3.5018; 00411 longPar[1][5] = std::max(0.0,0.6151 - 0.0561*std::log(einc)); 00412 00413 longSigma[1][0] = 0.01; 00414 longSigma[1][1] = 0.44; 00415 longSigma[1][2] = 0.01; 00416 longSigma[1][3] = 0.01; 00417 longSigma[1][4] = 0.20; 00418 longSigma[1][5] = 0.04; 00419 00420 longPar[2][0] = std::max(0.0,-1.55624e+01+1.56831e+01*std::tanh(5.93651e-01*(std::log(einc) + 4.89902))); 00421 longPar[2][1] = std::max(0.0,7.28995e-01+ 7.71148e-01*std::tanh(4.77898e-01*(std::log(einc) - 1.69087))); 00422 longPar[2][2] = std::max(0.0,1.23387+ 7.34778e-01*std::tanh(-3.14958e-01*(std::log(einc) - 0.529206))); 00423 longPar[2][3] = std::max(0.0,1.02070e+02+1.01873e+02*std::tanh(-4.99805e-01*(std::log(einc) + 5.04012))); 00424 longPar[2][4] = std::max(0.0,3.59765+8.53358e-01*std::tanh( 8.47277e-01*(std::log(einc) - 3.36548))); 00425 longPar[2][5] = std::max(0.0,4.27294e-01+1.62535e-02*std::tanh(-2.26278*(std::log(einc) - 1.81308))); 00426 00427 longSigma[2][0] = 0.01; 00428 longSigma[2][1] = 0.44; 00429 longSigma[2][2] = 0.01; 00430 longSigma[2][3] = 0.01; 00431 longSigma[2][4] = 0.20; 00432 longSigma[2][5] = 0.04; 00433 00434 double normalZ[Gflash::NxN]; 00435 for (int i = 0; i < Gflash::NxN ; i++) normalZ[i] = theRandGauss->fire(); 00436 00437 for(int k = 0 ; k < Gflash::NRegion ; k++) { 00438 for(int i = 0 ; i < Gflash::NxN ; i++) { 00439 double correlationSum = 0.0; 00440 for(int j = 0 ; j < Gflash::NxN ; j++) { 00441 correlationSum += correlationVector[Gflash::NStart[Gflash::NRegion]+(i+1)/2+j]*normalZ[i]; 00442 } 00443 longPar[k][i] = std::max(0.0,longPar[k][i]+longSigma[k][i]*correlationSum); 00444 } 00445 } 00446 00447 // parameters for the lateral profile 00448 00449 lateralPar[0] = 0.20; 00450 lateralPar[1] = std::max(0.0,0.40 -0.06*std::log(einc)); 00451 lateralPar[2] = 0.70 - 0.05*std::max(0.,std::log(einc)); 00452 lateralPar[3] = 0.20 * lateralPar[2]; 00453 }
| G4double GflashHadronShowerProfile::longitudinalProfile | ( | G4double | showerDepth, | |
| G4double | pathLength, | |||
| G4double | transDepth | |||
| ) | [private] |
Definition at line 455 of file GflashHadronShowerProfile.cc.
References funct::abs(), funct::exp(), GflashTrajectory::getGflashTrajectoryPoint(), GflashTrajectoryPoint::getPosition(), Gflash::kENCA, Gflash::kHB, Gflash::kHE, longPar, max, funct::pow(), Gflash::Rmin, Gflash::ScaleSensitive, showerType, theHelix, x, and Gflash::Zmin.
Referenced by hadronicParameterization().
00455 { 00456 00457 G4double heightProfile = 0; 00458 00459 // Energy in a delta step (dz) = (energy to deposite)*[Gamma(z+dz)-Gamma(z)]*dz 00460 // where the incomplete Gamma function gives an intergrate probability of the longitudinal 00461 // shower u[ to the shower depth (z). 00462 // Instead, we use approximated energy; energy in dz = (energy to deposite)*gamma(z)*dz 00463 // where gamma is the Gamma-distributed probability function 00464 00465 Genfun::LogGamma lgam; 00466 GflashTrajectoryPoint tempPoint; 00467 theHelix->getGflashTrajectoryPoint(tempPoint,pathLength); 00468 00469 double x = 0.0; 00470 //get parameters 00471 if(showerType == 1 || showerType == 2 ) { 00472 // std::cout << " pathLength tempPoint.getPosition().getRho()= " << pathLength << " " << tempPoint.getPosition().getRho() << std::endl; 00473 if(tempPoint.getPosition().getRho() < 150.0 ) { 00474 x = showerDepth*longPar[0][2]; 00475 heightProfile = longPar[0][0]*std::pow(x,longPar[0][1]-1.0)*std::exp(-x)/std::exp(lgam(longPar[0][1]))+longPar[0][3]; 00476 } 00477 else if (tempPoint.getPosition().getRho() > Gflash::Rmin[Gflash::kHB] ){ 00478 x = showerDepth; 00479 heightProfile = longPar[0][4]*std::exp(-x/longPar[0][5]); 00480 heightProfile *= Gflash::ScaleSensitive; 00481 } 00482 else heightProfile = 0.; 00483 } 00484 else if(showerType == 5 || showerType == 6){ 00485 //@@@use new parameterization for EE/HE 00486 if(std::abs(tempPoint.getPosition().getZ()) < Gflash::Zmin[Gflash::kENCA]+23.0 ) { 00487 x = showerDepth*longPar[0][2]; 00488 heightProfile = longPar[0][0]*std::pow(x,longPar[0][1]-1.0)*std::exp(-x)/std::exp(lgam(longPar[0][1]))+longPar[0][3]; 00489 } 00490 else if (std::abs(tempPoint.getPosition().getZ()) > Gflash::Rmin[Gflash::kHE] ){ 00491 x = showerDepth; 00492 heightProfile = longPar[0][4]*std::exp(-x/longPar[0][5]); 00493 heightProfile *= Gflash::ScaleSensitive; 00494 } 00495 else heightProfile = 0.; 00496 } 00497 else if (showerType == 3 || showerType == 7 ) { 00498 //two gammas between crystal and Hcal 00499 if((showerDepth - transDepth) > 0.0) { 00500 double x1 = (showerDepth-transDepth)*longPar[1][2]/16.42; 00501 double x2 = (showerDepth-transDepth)*longPar[1][5]/1.49; 00502 00503 heightProfile = longPar[1][3]*std::pow(x1,longPar[1][1]-1.0)*std::exp(-x1)/std::exp(lgam(longPar[1][1])) 00504 + (1.0-longPar[1][3])*std::pow(x2,longPar[1][4]-1.0)*std::exp(-x2)/std::exp(lgam(longPar[1][4])); 00505 heightProfile = std::max(0.0,longPar[1][0]*heightProfile); 00506 heightProfile *= Gflash::ScaleSensitive; 00507 } 00508 else heightProfile = 0.; 00509 } 00510 else if (showerType == 4 || showerType == 8 ) { 00511 //two gammas inside Hcal 00512 double x1 = showerDepth*longPar[2][2]/16.42; 00513 double x2 = showerDepth*longPar[2][5]/1.49; 00514 heightProfile = longPar[2][3]*std::pow(x1,longPar[2][1]-1.0)*std::exp(-x1)/std::exp(lgam(longPar[2][1])) 00515 + (1.0-longPar[2][3])*std::pow(x2,longPar[2][4]-1.0)*std::exp(-x2)/std::exp(lgam(longPar[2][4])); 00516 heightProfile = std::max(0.0,longPar[2][0]*heightProfile); 00517 heightProfile *= Gflash::ScaleSensitive; 00518 } 00519 00520 return heightProfile; 00521 }
| void GflashHadronShowerProfile::samplingFluctuation | ( | G4double & | de, | |
| G4double | einc | |||
| ) | [private] |
Definition at line 523 of file GflashHadronShowerProfile.cc.
References energyToDeposit, i, jCalorimeter, Gflash::NDET, funct::pow(), Gflash::SAMHAD, and theRandGamma.
00523 { 00524 00525 G4double spot[Gflash::NDET]; 00526 00527 for(G4int i = 0 ; i < Gflash::NDET ; i++) { 00528 spot[i] = std::pow(Gflash::SAMHAD[0][i],2) // resolution 00529 + std::pow(Gflash::SAMHAD[1][i],2)/einc // noisy 00530 + std::pow(Gflash::SAMHAD[2][i],2)*einc; // constant 00531 } 00532 G4double ein = de * (energyToDeposit/einc); 00533 00534 de = (ein > 0 ) ? 00535 theRandGamma->fire(ein/spot[jCalorimeter],1.0)*spot[jCalorimeter] : ein; 00536 }
std::vector<GflashEnergySpot> GflashHadronShowerProfile::aEnergySpotList [private] |
Definition at line 45 of file GflashHadronShowerProfile.h.
Referenced by clearSpotList(), getEnergySpotList(), and hadronicParameterization().
G4double GflashHadronShowerProfile::correlationVector[Gflash::NRegion *Gflash::NxN *(Gflash::NxN+1)/2] [private] |
Definition at line 49 of file GflashHadronShowerProfile.h.
Referenced by fillFluctuationVector(), and loadParameters().
G4double GflashHadronShowerProfile::energyToDeposit [private] |
Definition at line 47 of file GflashHadronShowerProfile.h.
Referenced by hadronicParameterization(), loadParameters(), and samplingFluctuation().
Definition at line 44 of file GflashHadronShowerProfile.h.
Referenced by getCalorimeterNumber(), GflashHadronShowerProfile(), hadronicParameterization(), insideSampling(), loadParameters(), and samplingFluctuation().
G4double GflashHadronShowerProfile::lateralPar[4] [private] |
Definition at line 53 of file GflashHadronShowerProfile.h.
Referenced by hadronicParameterization(), and loadParameters().
G4double GflashHadronShowerProfile::longPar[Gflash::NRegion][Gflash::NxN] [private] |
Definition at line 50 of file GflashHadronShowerProfile.h.
Referenced by loadParameters(), and longitudinalProfile().
G4double GflashHadronShowerProfile::longSigma[Gflash::NRegion][Gflash::NxN] [private] |
G4int GflashHadronShowerProfile::showerType [private] |
Definition at line 43 of file GflashHadronShowerProfile.h.
Referenced by GflashHadronShowerProfile(), hadronicParameterization(), loadParameters(), and longitudinalProfile().
Definition at line 60 of file GflashHadronShowerProfile.h.
Referenced by GflashHadronShowerProfile(), hadronicParameterization(), longitudinalProfile(), and ~GflashHadronShowerProfile().
Definition at line 59 of file GflashHadronShowerProfile.h.
Referenced by GflashHadronShowerProfile(), and hadronicParameterization().
CLHEP::RandGamma* GflashHadronShowerProfile::theRandGamma [private] |
Definition at line 63 of file GflashHadronShowerProfile.h.
Referenced by GflashHadronShowerProfile(), samplingFluctuation(), and ~GflashHadronShowerProfile().
CLHEP::RandGaussQ* GflashHadronShowerProfile::theRandGauss [private] |
Definition at line 62 of file GflashHadronShowerProfile.h.
Referenced by GflashHadronShowerProfile(), hadronicParameterization(), loadParameters(), and ~GflashHadronShowerProfile().
1.5.4