Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id$ 26 // 27 // 27 // 28 // 28 // ------------------------------------------- 29 // ------------------------------------------------------------ 29 // GEANT 4 class implementation 30 // GEANT 4 class implementation 30 // 31 // 31 // ------- GFlashSamplingShowerParameteri 32 // ------- GFlashSamplingShowerParameterisation ------- 32 // 33 // 33 // Authors: E.Barberio & Joanna Weng - 11.2005 34 // Authors: E.Barberio & Joanna Weng - 11.2005 34 // ------------------------------------------- 35 // ------------------------------------------------------------ 35 36 36 #include <cmath> 37 #include <cmath> 37 38 38 #include "GFlashSamplingShowerParameterisation 39 #include "GFlashSamplingShowerParameterisation.hh" 39 #include "GVFlashShowerParameterisation.hh" 40 #include "GVFlashShowerParameterisation.hh" 40 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh" 41 #include "Randomize.hh" 42 #include "Randomize.hh" 42 #include "G4ios.hh" 43 #include "G4ios.hh" 43 #include "G4Material.hh" 44 #include "G4Material.hh" 44 #include "G4MaterialTable.hh" 45 #include "G4MaterialTable.hh" 45 46 46 GFlashSamplingShowerParameterisation::GFlashSa << 47 GFlashSamplingShowerParameterisation:: 47 G4Material* aMat1, G4Material* aMat2, G4doub << 48 GFlashSamplingShowerParameterisation(G4Material* aMat1, G4Material* aMat2, 48 GFlashSamplingShowerTuning* aPar) << 49 G4double dd1, G4double dd2, >> 50 GFlashSamplingShowerTuning* aPar) 49 : GVFlashShowerParameterisation(), 51 : GVFlashShowerParameterisation(), 50 ParAveT2(0.), << 52 ParAveT2(0.), ParSigLogT1(0.), ParSigLogT2(0.), 51 ParSigLogT1(0.), << 53 ParSigLogA1(0.), ParSigLogA2(0.), ParRho1(0.), ParRho2(0.), ParsAveA2(0.), 52 ParSigLogT2(0.), << 54 AveLogAlphah(0.), AveLogTmaxh(0.), SigmaLogAlphah(0.), SigmaLogTmaxh(0.), 53 ParSigLogA1(0.), << 55 Rhoh(0.), Alphah(0.), Tmaxh(0.), Betah(0.), AveLogAlpha(0.), AveLogTmax(0.), 54 ParSigLogA2(0.), << 56 SigmaLogAlpha(0.), SigmaLogTmax(0.), Rho(0.), Alpha(0.), Tmax(0.), Beta(0.) 55 ParRho1(0.), << 57 { 56 ParRho2(0.), << 58 if(!aPar) { thePar = new GFlashSamplingShowerTuning; owning = true; } 57 ParsAveA2(0.), << 59 else { thePar = aPar; owning = false; } 58 AveLogAlphah(0.), << 60 59 AveLogTmaxh(0.), << 61 SetMaterial(aMat1,aMat2 ); 60 SigmaLogAlphah(0.), << 62 d1=dd1; 61 SigmaLogTmaxh(0.), << 63 d2=dd2; 62 Rhoh(0.), << 63 Alphah(0.), << 64 Tmaxh(0.), << 65 Betah(0.), << 66 AveLogAlpha(0.), << 67 AveLogTmax(0.), << 68 SigmaLogAlpha(0.), << 69 SigmaLogTmax(0.), << 70 Rho(0.), << 71 Alpha(0.), << 72 Tmax(0.), << 73 Beta(0.) << 74 { << 75 if (!aPar) { << 76 thePar = new GFlashSamplingShowerTuning; << 77 } << 78 else { << 79 thePar = aPar; << 80 } << 81 << 82 SetMaterial(aMat1, aMat2); << 83 d1 = dd1; << 84 d2 = dd2; << 85 64 86 // Longitudinal Coefficients for a homogenio 65 // Longitudinal Coefficients for a homogenious calo 87 66 88 // shower max 67 // shower max 89 ParAveT1 = thePar->ParAveT1(); // ln (ln y << 68 ParAveT1 = thePar->ParAveT1(); // ln (ln y -0.812) 90 ParAveA1 = thePar->ParAveA1(); // ln a (0.8 << 69 ParAveA1 = thePar->ParAveA1(); // ln a (0.81 + (0.458 + 2.26/Z)ln y) 91 ParAveA2 = thePar->ParAveA2(); << 70 ParAveA2 = thePar->ParAveA2(); 92 ParAveA3 = thePar->ParAveA3(); << 71 ParAveA3 = thePar->ParAveA3(); 93 // Variance of shower max sampling << 94 ParSigLogT1 = << 95 thePar << 96 ->ParsSigLogT1(); // Sigma T1 (-1.4 + 1 << 97 ParSigLogT2 = thePar->ParsSigLogT2(); // le << 98 // variance of 'alpha' << 99 ParSigLogA1 = << 100 thePar << 101 ->ParSigLogA1(); // Sigma a (-0.58 + 0. << 102 ParSigLogA2 = thePar->ParSigLogA2(); // lea << 103 // correlation alpha%T << 104 ParRho1 = thePar->ParRho1(); // Rho = 0.705 << 105 ParRho2 = thePar->ParRho2(); // leaving Par << 106 // Sampling 72 // Sampling 107 ParsAveT1 = thePar->ParsAveT1(); // T_sam = << 73 ParsAveT1 = thePar->ParsAveT1(); // T_sam = log(exp( log T_hom) + t1*Fs-1 + t2*(1-ehat)); 108 ParsAveT2 = thePar->ParsAveT2(); << 74 ParsAveT2 = thePar->ParsAveT2(); 109 ParsAveA1 = thePar->ParsAveA1(); << 75 ParsAveA1 = thePar->ParsAveA1(); 110 // Variance of shower max sampling << 76 // Variance of shower max sampling 111 ParsSigLogT1 = thePar->ParsSigLogT1(); // S << 77 ParsSigLogT1 = thePar->ParSigLogT1(); // Sigma T1 (-2.5 + 1.25 ln y)**-1 112 // w << 78 ParsSigLogT2 = thePar->ParSigLogT2(); 113 ParsSigLogT2 = thePar->ParsSigLogT2(); << 114 // variance of 'alpha' 79 // variance of 'alpha' 115 ParsSigLogA1 = thePar->ParsSigLogA1(); // S << 80 ParsSigLogA1 = thePar->ParSigLogA1(); // Sigma a (-0.82 + 0.79 ln y)**-1 116 // w << 81 ParsSigLogA2 = thePar->ParSigLogA2(); 117 ParsSigLogA2 = thePar->ParsSigLogA2(); << 118 // correlation alpha%T 82 // correlation alpha%T 119 ParsRho1 = << 83 ParsRho1 = thePar->ParRho1(); // Rho = 0.784 -0.023 ln y 120 thePar->ParsRho1(); // Rho = 0.784 -0.023 << 84 ParsRho2 = thePar->ParRho2(); 121 ParsRho2 = thePar->ParsRho2(); << 122 85 123 // Radial Coefficients 86 // Radial Coefficients 124 // r_C (tau)= z_1 +z_2 tau 87 // r_C (tau)= z_1 +z_2 tau 125 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+st 88 // r_t (tau)= k1 (std::exp (k3(tau -k2 ))+std::exp (k_4 (tau- k_2)))) 126 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 << 89 ParRC1 = thePar->ParRC1(); // z_1 = 0.0251 + 0.00319 ln E 127 ParRC2 = thePar->ParRC2(); << 90 ParRC2 = thePar->ParRC2(); 128 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 << 91 ParRC3 = thePar->ParRC3(); // z_2 = 0.1162 + - 0.000381 Z 129 ParRC4 = thePar->ParRC4(); << 92 ParRC4 = thePar->ParRC4(); 130 93 131 ParWC1 = thePar->ParWC1(); 94 ParWC1 = thePar->ParWC1(); 132 ParWC2 = thePar->ParWC2(); 95 ParWC2 = thePar->ParWC2(); 133 ParWC3 = thePar->ParWC3(); 96 ParWC3 = thePar->ParWC3(); 134 ParWC4 = thePar->ParWC4(); 97 ParWC4 = thePar->ParWC4(); 135 ParWC5 = thePar->ParWC5(); << 98 ParWC5 = thePar->ParWC5(); 136 ParWC6 = thePar->ParWC6(); 99 ParWC6 = thePar->ParWC6(); 137 ParRT1 = thePar->ParRT1(); 100 ParRT1 = thePar->ParRT1(); 138 ParRT2 = thePar->ParRT2(); 101 ParRT2 = thePar->ParRT2(); 139 ParRT3 = thePar->ParRT3(); 102 ParRT3 = thePar->ParRT3(); 140 ParRT4 = thePar->ParRT4(); << 103 ParRT4 = thePar->ParRT4(); 141 ParRT5 = thePar->ParRT5(); 104 ParRT5 = thePar->ParRT5(); 142 ParRT6 = thePar->ParRT6(); 105 ParRT6 = thePar->ParRT6(); 143 106 144 // additional sampling parameter << 107 //additional sampling parameter 145 ParsRC1 = thePar->ParsRC1(); << 108 ParsRC1= thePar->ParsRC1(); 146 ParsRC2 = thePar->ParsRC2(); << 109 ParsRC2= thePar->ParsRC2(); 147 ParsWC1 = thePar->ParsWC1(); << 110 ParsWC1= thePar->ParsWC1(); 148 ParsWC2 = thePar->ParsWC2(); << 111 ParsWC2= thePar->ParsWC2(); 149 ParsRT1 = thePar->ParsRT1(); << 112 ParsRT1= thePar->ParsRT1(); 150 ParsRT2 = thePar->ParsRT2(); << 113 ParsRT2= thePar->ParsRT2(); 151 114 152 // Coeff for fluctuedted radial profiles fo 115 // Coeff for fluctuedted radial profiles for a sampling media 153 ParsSpotT1 = thePar->ParSpotT1(); // T_spot << 116 ParsSpotT1 = thePar->ParSpotT1(); // T_spot = T_hom =(0.698 + 0.00212) 154 ParsSpotT2 = thePar->ParSpotT2(); << 117 ParsSpotT2 = thePar->ParSpotT2(); 155 ParsSpotA1 = thePar->ParSpotA1(); // a_spot << 118 ParsSpotA1 = thePar->ParSpotA1(); // a_spot= a_hom (0.639 + 0.00334) 156 ParsSpotA2 = thePar->ParSpotA2(); << 119 ParsSpotA2 = thePar->ParSpotA2(); 157 ParsSpotN1 = thePar->ParSpotN1(); // N_Spot << 120 ParsSpotN1 = thePar->ParSpotN1(); // N_Spot 93 * ln(Z) E ** 0.876 158 ParsSpotN2 = thePar->ParSpotN2(); << 121 ParsSpotN2 = thePar->ParSpotN2(); 159 SamplingResolution = thePar->SamplingResolut << 122 SamplingResolution = thePar->SamplingResolution(); 160 ConstantResolution = thePar->ConstantResolut << 123 ConstantResolution = thePar->ConstantResolution(); 161 NoiseResolution = thePar->NoiseResolution(); << 124 NoiseResolution = thePar->NoiseResolution(); 162 125 163 // Inits 126 // Inits 164 NSpot = 0.00; << 127 NSpot = 0.00; 165 AlphaNSpot = 0.00; << 128 AlphaNSpot = 0.00; 166 TNSpot = 0.00; << 129 TNSpot = 0.00; 167 BetaNSpot = 0.00; << 130 BetaNSpot = 0.00; 168 RadiusCore = 0.00; << 131 RadiusCore = 0.00; 169 WeightCore = 0.00; << 132 WeightCore = 0.00; 170 RadiusTail = 0.00; << 133 RadiusTail = 0.00; 171 ComputeZAX0EFFetc(); 134 ComputeZAX0EFFetc(); 172 135 173 G4cout << "/******************************** 136 G4cout << "/********************************************/ " << G4endl; 174 G4cout << " - GFlashSamplingShowerParameter 137 G4cout << " - GFlashSamplingShowerParameterisation::Constructor - " << G4endl; 175 G4cout << "/******************************** << 138 G4cout << "/********************************************/ " << G4endl; 176 } 139 } 177 140 178 // ------------------------------------------- 141 // ------------------------------------------------------------ 179 142 180 GFlashSamplingShowerParameterisation::~GFlashS 143 GFlashSamplingShowerParameterisation::~GFlashSamplingShowerParameterisation() 181 { 144 { 182 delete thePar; << 145 if(owning) { delete thePar; } 183 } 146 } 184 147 185 // ------------------------------------------- 148 // ------------------------------------------------------------ 186 149 187 void GFlashSamplingShowerParameterisation::Set << 150 void GFlashSamplingShowerParameterisation:: >> 151 SetMaterial(G4Material *mat1, G4Material *mat2) 188 { 152 { 189 const G4double Es = 21.2 * MeV; << 153 G4double Es = 21*MeV; 190 material1 = mat1; << 154 material1= mat1; 191 Z1 = GetEffZ(material1); 155 Z1 = GetEffZ(material1); 192 A1 = GetEffA(material1); 156 A1 = GetEffA(material1); 193 density1 = material1->GetDensity(); 157 density1 = material1->GetDensity(); 194 X01 = material1->GetRadlen(); << 158 X01 = material1->GetRadlen(); 195 Ec1 = 2.66 * std::pow((X01 * Z1 / A1), 1.1); << 159 Ec1 = 2.66 * std::pow((X01 * Z1 / A1),1.1); 196 // Ec1 = 610.0 * MeV / (Z1 + 1.24); << 160 Rm1 = X01*Es/Ec1; 197 Rm1 = X01 * Es / Ec1; << 198 161 199 material2 = mat2; << 162 material2= mat2; 200 Z2 = GetEffZ(material2); 163 Z2 = GetEffZ(material2); 201 A2 = GetEffA(material2); 164 A2 = GetEffA(material2); 202 density2 = material2->GetDensity(); 165 density2 = material2->GetDensity(); 203 X02 = material2->GetRadlen(); << 166 X02 = material2->GetRadlen(); 204 Ec2 = 2.66 * std::pow((X02 * Z2 / A2), 1.1); << 167 Ec2 = 2.66 * std::pow((X02 * Z2 / A2),1.1); 205 // Ec2 = 610.0 * MeV / (Z2 + 1.24); << 168 Rm2 = X02*Es/Ec2; 206 Rm2 = X02 * Es / Ec2; << 169 // PrintMaterial(); 207 // PrintMaterial(); << 208 } 170 } 209 171 210 // ------------------------------------------- 172 // ------------------------------------------------------------ 211 173 212 void GFlashSamplingShowerParameterisation::Com 174 void GFlashSamplingShowerParameterisation::ComputeZAX0EFFetc() 213 { 175 { 214 G4cout << "/************ ComputeZAX0EFFetc * 176 G4cout << "/************ ComputeZAX0EFFetc ************/" << G4endl; 215 G4cout << " - GFlashSamplingShowerParameter 177 G4cout << " - GFlashSamplingShowerParameterisation::Material - " << G4endl; 216 178 217 const G4double Es = 21.2 * MeV; << 179 G4double Es = 21*MeV; //constant 218 180 219 // material and geometry parameters for a sa 181 // material and geometry parameters for a sampling calorimeter 220 G4double denominator = (d1 * density1 + d2 * << 182 G4double denominator = (d1*density1 + d2*density2); 221 G4double W1 = (d1 * density1) / denominator; << 183 G4double W1 = (d1*density1) / denominator; 222 G4double W2 = (d2 * density2) / denominator; << 184 G4double W2 = (d2*density2)/denominator; 223 Zeff = (W1 * Z1) + (W2 * Z2); // X0*Es/Ec; << 185 Zeff = ( W1*Z2 ) + (W2*Z1); //X0*Es/Ec; 224 Aeff = (W1 * A1) + (W2 * A2); << 186 Aeff = ( W1*A1 ) + (W2*A2); 225 Rhoeff = ((d1 * density1) + (d2 * density2)) << 187 X0eff =(1/ (( W1 / X01) +( W2 / X02))); 226 X0eff = (W1 * Rhoeff) / (X01 * density1) + ( << 188 Rhoeff = ( (d1 *density1 ) + (d2 * density2 ))/G4double (d2 + d1 ); 227 X0eff = 1. / X0eff; << 189 Rmeff = 1/ ((((W1*Ec1)/ X01) + ((W2* Ec2)/ X02) ) / Es ) ; 228 Rmeff = 1. / ((((W1 * Ec1) / X01) + ((W2 * E << 190 Eceff = X0eff *((W1*Ec1)/ X01 + (W2* Ec2)/ X02 ); 229 Eceff = X0eff * ((W1 * Ec1) / X01 + (W2 * Ec << 191 Fs = X0eff/G4double ((d1/mm )+(d2/mm) ); 230 Fs = X0eff / (d1 + d2); // --> was G4double << 192 ehat = (1. / (1+ 0.007*(Z1- Z2))); 231 // mm makes sense.. << 232 ehat = (1. / (1 + 0.007 * (Z1 - Z2))); << 233 193 234 G4cout << "W1= " << W1 << G4endl; << 194 G4cout << "W1= " << W1 << G4endl; 235 G4cout << "W2= " << W2 << G4endl; 195 G4cout << "W2= " << W2 << G4endl; 236 G4cout << "effective quantities Zeff = " << << 196 G4cout << "effective quantities Zeff = "<<Zeff<< G4endl; 237 G4cout << "effective quantities Aeff = " << << 197 G4cout << "effective quantities Aeff = "<<Aeff<< G4endl; 238 G4cout << "effective quantities Rhoeff = " < << 198 G4cout << "effective quantities Rhoeff = "<<Rhoeff/g *cm3<<" g/cm3" << G4endl; 239 G4cout << "effective quantities X0eff = " << << 199 G4cout << "effective quantities X0eff = "<<X0eff/cm <<" cm" << G4endl; 240 200 241 X0eff = X0eff * Rhoeff; 201 X0eff = X0eff * Rhoeff; 242 202 243 G4cout << "effective quantities X0eff = " << << 203 G4cout << "effective quantities X0eff = "<<X0eff/g*cm2 <<" g/cm2" << G4endl; 244 X0eff = X0eff / Rhoeff; << 204 X0eff = X0eff /Rhoeff; 245 G4cout << "effective quantities RMeff = " << << 205 G4cout << "effective quantities RMeff = "<<Rmeff/cm<<" cm" << G4endl; 246 Rmeff = Rmeff * Rhoeff; << 206 Rmeff = Rmeff* Rhoeff; 247 G4cout << "effective quantities RMeff = " << << 207 G4cout << "effective quantities RMeff = "<<Rmeff/g *cm2<<" g/cm2" << G4endl; 248 Rmeff = Rmeff / Rhoeff; << 208 Rmeff = Rmeff/ Rhoeff; 249 G4cout << "effective quantities Eceff = " << << 209 G4cout << "effective quantities Eceff = "<<Eceff/MeV<< " MeV"<< G4endl; 250 G4cout << "effective quantities Fs = " << Fs << 210 G4cout << "effective quantities Fs = "<<Fs<<G4endl; 251 G4cout << "effective quantities ehat = " << << 211 G4cout << "effective quantities ehat = "<<ehat<<G4endl; 252 G4cout << "/******************************** << 212 G4cout << "/********************************************/ " <<G4endl; 253 } 213 } 254 214 255 // ------------------------------------------- 215 // ------------------------------------------------------------ 256 216 257 void GFlashSamplingShowerParameterisation::Gen << 217 void GFlashSamplingShowerParameterisation:: >> 218 GenerateLongitudinalProfile(G4double Energy) 258 { 219 { 259 if ((material1 == 0) || (material2 == 0)) { << 220 if ((material1==0) || (material2 ==0)) >> 221 { 260 G4Exception("GFlashSamplingShowerParameter 222 G4Exception("GFlashSamplingShowerParameterisation::GenerateLongitudinalProfile()", 261 "InvalidSetup", FatalException 223 "InvalidSetup", FatalException, "No material initialized!"); 262 } << 224 } 263 G4double y = Energy / Eceff; << 225 G4double y = Energy/Eceff; 264 ComputeLongitudinalParameters(y); << 226 ComputeLongitudinalParameters(y); 265 GenerateEnergyProfile(y); 227 GenerateEnergyProfile(y); 266 GenerateNSpotProfile(y); 228 GenerateNSpotProfile(y); 267 } 229 } 268 230 269 // ------------------------------------------- 231 // ------------------------------------------------------------ 270 232 271 void GFlashSamplingShowerParameterisation::Com << 233 void >> 234 GFlashSamplingShowerParameterisation::ComputeLongitudinalParameters(G4double y) 272 { 235 { 273 AveLogTmaxh = std::log(std::max(ParAveT1 + s << 236 AveLogTmaxh = std::log(std::max(ParAveT1 +std::log(y),0.1)); //ok 274 AveLogAlphah = << 237 AveLogAlphah = std::log(std::max(ParAveA1 + (ParAveA2+ParAveA3/Zeff)*std::log(y),.1)); //ok 275 std::log(std::max(ParAveA1 + (ParAveA2 + P << 238 //hom 276 // hom << 239 SigmaLogTmaxh = std::min(0.5,1.00/( ParSigLogT1 + ParSigLogT2*std::log(y)) ); //ok 277 SigmaLogTmaxh = std::min(0.5, 1.00 / (ParSig << 240 SigmaLogAlphah = std::min(0.5,1.00/( ParSigLogA1 + ParSigLogA2*std::log(y))); //ok 278 SigmaLogAlphah = std::min(0.5, 1.00 / (ParSi << 241 Rhoh = ParRho1+ParRho2*std::log(y);//ok 279 Rhoh = ParRho1 + ParRho2 * std::log(y); // << 242 // if sampling 280 // if sampling << 243 AveLogTmax = std::max(0.1,std::log(std::exp(AveLogTmaxh) 281 AveLogTmax = << 244 + ParsAveT1/Fs + ParsAveT2*(1-ehat))); //ok 282 std::max(0.1, std::log(std::exp(AveLogTmax << 245 AveLogAlpha = std::max(0.1,std::log(std::exp(AveLogAlphah) 283 AveLogAlpha = std::max(0.1, std::log(std::ex << 246 + (ParsAveA1/Fs))); //ok 284 // 247 // 285 SigmaLogTmax = std::min(0.5, 1.00 / (ParsSig << 248 SigmaLogTmax = std::min(0.5,1.00/( ParsSigLogT1 286 SigmaLogAlpha = std::min(0.5, 1.00 / (ParsSi << 249 + ParsSigLogT2*std::log(y)) ); //ok 287 Rho = ParsRho1 + ParsRho2 * std::log(y); // << 250 SigmaLogAlpha = std::min(0.5,1.00/( ParsSigLogA1 288 << 251 + ParsSigLogA2*std::log(y))); //ok 289 if (0) { << 252 Rho = ParsRho1+ParsRho2*std::log(y); //ok 290 G4cout << " y = " << y << G4end << 291 G4cout << " std::log(std::exp(AveLogTmaxh) << 292 << " std::log(" << std::exp(AveLogT << 293 << ParsAveT2 * (1 - ehat) << ") = " << 294 << " std::log(" << std::exp(AveLogT << 295 << ParsAveT2 << "*" << (1 - ehat) < << 296 << " std::log(" << std::exp(AveLogT << 297 << G4endl; << 298 G4cout << " AveLogTmaxh " << AveLogTmax << 299 G4cout << " AveLogAlphah " << AveLogAlph << 300 G4cout << " SigmaLogTmaxh " << SigmaLogTm << 301 G4cout << " 1.00/( ParSigLogT1 + ParSigLog << 302 << (ParSigLogT1 + ParSigLogT2 * std << 303 << ParSigLogT1 << " + " << ParSigLo << 304 << ParSigLogT1 << " + " << ParSigLo << 305 G4cout << " SigmaLogAlphah " << SigmaLogAl << 306 G4cout << " Rhoh " << Rhoh << G4 << 307 G4cout << " AveLogTmax " << AveLogTmax << 308 G4cout << " AveLogAlpha " << AveLogAlph << 309 G4cout << " SigmaLogTmax " << SigmaLogTm << 310 G4cout << " SigmaLogAlpha " << SigmaLogAl << 311 G4cout << " Rho " << Rho << G4e << 312 } << 313 } 253 } 314 254 315 // ------------------------------------------- 255 // ------------------------------------------------------------ 316 256 317 void GFlashSamplingShowerParameterisation::Gen 257 void GFlashSamplingShowerParameterisation::GenerateEnergyProfile(G4double /* y */) 318 { << 258 { 319 G4double Correlation1 = std::sqrt((1 + Rho) << 259 G4double Correlation1 = std::sqrt((1+Rho)/2); 320 G4double Correlation2 = std::sqrt((1 - Rho) << 260 G4double Correlation2 = std::sqrt((1-Rho)/2); 321 G4double Correlation1h = std::sqrt((1 + Rhoh << 261 G4double Correlation1h = std::sqrt((1+Rhoh)/2); 322 G4double Correlation2h = std::sqrt((1 - Rhoh << 262 G4double Correlation2h = std::sqrt((1-Rhoh)/2); 323 G4double Random1 = G4RandGauss::shoot(); 263 G4double Random1 = G4RandGauss::shoot(); 324 G4double Random2 = G4RandGauss::shoot(); 264 G4double Random2 = G4RandGauss::shoot(); 325 265 326 Tmax = std::max( << 266 Tmax = std::max(1.,std::exp( AveLogTmax + SigmaLogTmax * 327 1., std::exp(AveLogTmax + SigmaLogTmax * ( << 267 (Correlation1*Random1 + Correlation2*Random2) )); 328 Alpha = std::max( << 268 Alpha = std::max(1.1,std::exp( AveLogAlpha + SigmaLogAlpha * 329 1.1, std::exp(AveLogAlpha + SigmaLogAlpha << 269 (Correlation1*Random1 - Correlation2*Random2) )); 330 Beta = (Alpha - 1.00) / Tmax; << 270 Beta = (Alpha-1.00)/Tmax; 331 // Parameters for Enenrgy Profile including << 271 //Parameters for Enenrgy Profile including correaltion and sigmas 332 Tmaxh = << 272 Tmaxh = std::exp( AveLogTmaxh + SigmaLogTmaxh * 333 std::exp(AveLogTmaxh + SigmaLogTmaxh * (Co << 273 (Correlation1h*Random1 + Correlation2h*Random2) ); 334 Alphah = << 274 Alphah = std::exp( AveLogAlphah + SigmaLogAlphah * 335 std::exp(AveLogAlphah + SigmaLogAlphah * ( << 275 (Correlation1h*Random1 - Correlation2h*Random2) ); 336 Betah = (Alphah - 1.00) / Tmaxh; << 276 Betah = (Alphah-1.00)/Tmaxh; 337 } 277 } 338 278 339 // ------------------------------------------- 279 // ------------------------------------------------------------ 340 280 341 void GFlashSamplingShowerParameterisation::Gen 281 void GFlashSamplingShowerParameterisation::GenerateNSpotProfile(const G4double y) 342 { 282 { 343 TNSpot = Tmaxh * (ParsSpotT1 + ParsSpotT2 * << 283 TNSpot = Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff); //ok. 344 TNSpot = std::max(0.5, Tmaxh * (ParsSpotT1 + << 284 TNSpot = std::max(0.5,Tmaxh * (ParsSpotT1+ParsSpotT2*Zeff)); 345 AlphaNSpot = Alphah * (ParsSpotA1 + ParsSpot << 285 AlphaNSpot = Alphah * (ParsSpotA1+ParsSpotA2*Zeff); 346 BetaNSpot = (AlphaNSpot - 1.00) / TNSpot; / << 286 BetaNSpot = (AlphaNSpot-1.00)/TNSpot; // ok 347 NSpot = ParsSpotN1 / SamplingResolution * st << 287 NSpot = ParsSpotN1 /SamplingResolution * std::pow(y*Eceff/GeV,ParsSpotN2 ); 348 } 288 } 349 289 350 // ------------------------------------------- 290 // ------------------------------------------------------------ 351 291 352 G4double GFlashSamplingShowerParameterisation: << 292 G4double >> 293 GFlashSamplingShowerParameterisation:: >> 294 ApplySampling(const G4double DEne, const G4double ) 353 { 295 { 354 G4double DEneFluctuated = DEne; 296 G4double DEneFluctuated = DEne; 355 G4double Resolution = std::pow(SamplingResol << 297 G4double Resolution = std::pow(SamplingResolution,2); 356 298 357 // +pow(NoiseResolution,2)/ //@@@@@@@ << 299 // +pow(NoiseResolution,2)/ //@@@@@@@@ FIXME 358 // Energy*(1.*MeV)+ 300 // Energy*(1.*MeV)+ 359 // pow(ConstantResol 301 // pow(ConstantResolution,2)* 360 // Energy/(1.*MeV); 302 // Energy/(1.*MeV); 361 303 362 if (Resolution > 0.0 && DEne > 0.00) { << 304 if(Resolution >0.0 && DEne > 0.00) 363 // G4float x1 = DEne / Resolution; << 305 { 364 // G4float x2 = G4RandGamma::shoot(x1, 1.0 << 306 G4float x1=DEne/Resolution; 365 // DEneFluctuated = x2; << 307 G4float x2 = CLHEP::RandGamma::shoot(x1, 1.0)*Resolution; 366 G4double x1 = DEne / Resolution; << 308 DEneFluctuated=x2; 367 G4double x2 = 1.0 / Resolution; << 368 DEneFluctuated = G4RandGamma::shoot(x1, x2 << 369 } 309 } 370 return DEneFluctuated; 310 return DEneFluctuated; 371 } 311 } 372 312 373 // ------------------------------------------- 313 // ------------------------------------------------------------ 374 314 375 G4double GFlashSamplingShowerParameterisation: 315 G4double GFlashSamplingShowerParameterisation:: 376 IntegrateEneLongitudinal(G4double Longitudinal 316 IntegrateEneLongitudinal(G4double LongitudinalStep) 377 { 317 { 378 G4double LongitudinalStepInX0 = Longitudinal 318 G4double LongitudinalStepInX0 = LongitudinalStep / X0eff; 379 G4float x1= Betah*LongitudinalStepInX0; 319 G4float x1= Betah*LongitudinalStepInX0; 380 G4float x2= Alphah; 320 G4float x2= Alphah; 381 float x3 = gam(x1,x2); 321 float x3 = gam(x1,x2); 382 G4double DEne=x3; 322 G4double DEne=x3; 383 return DEne; 323 return DEne; 384 } 324 } 385 325 386 // ------------------------------------------- 326 // ------------------------------------------------------------ 387 327 388 G4double GFlashSamplingShowerParameterisation: << 328 G4double GFlashSamplingShowerParameterisation:: >> 329 IntegrateNspLongitudinal(G4double LongitudinalStep) 389 { 330 { 390 G4double LongitudinalStepInX0 = Longitudinal << 331 G4double LongitudinalStepInX0 = LongitudinalStep / X0eff; 391 G4float x1 = BetaNSpot * LongitudinalStepInX << 332 G4float x1 = BetaNSpot*LongitudinalStepInX0; 392 G4float x2 = AlphaNSpot; 333 G4float x2 = AlphaNSpot; 393 G4float x3 = gam(x1, x2); << 334 G4float x3 = gam(x1,x2); 394 G4double DNsp = x3; 335 G4double DNsp = x3; 395 return DNsp; 336 return DNsp; 396 } 337 } 397 338 398 // ------------------------------------------- 339 // ------------------------------------------------------------ 399 340 400 G4double GFlashSamplingShowerParameterisation: << 341 G4double GFlashSamplingShowerParameterisation:: 401 << 342 GenerateRadius(G4int ispot, G4double Energy, G4double LongitudinalPosition) 402 { 343 { 403 if (ispot < 1) { << 344 if(ispot < 1) >> 345 { 404 // Determine lateral parameters in the mid 346 // Determine lateral parameters in the middle of the step. 405 // They depend on energy & position along 347 // They depend on energy & position along step 406 // 348 // 407 G4double Tau = ComputeTau(LongitudinalPosi 349 G4double Tau = ComputeTau(LongitudinalPosition); 408 ComputeRadialParameters(Energy, Tau); << 350 ComputeRadialParameters(Energy,Tau); 409 } 351 } 410 << 352 411 G4double Radius; 353 G4double Radius; 412 G4double Random1 = G4UniformRand(); 354 G4double Random1 = G4UniformRand(); 413 G4double Random2 = G4UniformRand(); << 355 G4double Random2 = G4UniformRand(); 414 if (Random1 < WeightCore) // WeightCore = p << 356 if(Random1 <WeightCore) //WeightCore = p < w_i 415 { 357 { 416 Radius = Rmeff * RadiusCore * std::sqrt(Ra << 358 Radius = Rmeff * RadiusCore * std::sqrt( Random2/(1. - Random2) ); 417 } << 418 else { << 419 Radius = Rmeff * RadiusTail * std::sqrt(Ra << 420 } 359 } 421 Radius = std::min(Radius, DBL_MAX); << 360 else >> 361 { >> 362 Radius = Rmeff * RadiusTail * std::sqrt( Random2/(1. - Random2) ); >> 363 } >> 364 Radius = std::min(Radius,DBL_MAX); 422 return Radius; 365 return Radius; 423 } 366 } 424 367 425 // ------------------------------------------- 368 // ------------------------------------------------------------ 426 369 427 G4double GFlashSamplingShowerParameterisation: << 370 G4double 428 { << 371 GFlashSamplingShowerParameterisation:: 429 G4double tau = LongitudinalPosition / Tmax / << 372 ComputeTau(G4double LongitudinalPosition) 430 * (Alpha - 1.00) / Alpha * st << 373 { 431 / (std::exp(AveLogAlpha) - 1. << 374 G4double tau = LongitudinalPosition / Tmax/ X0eff //<t> = T* a /(a - 1) >> 375 * (Alpha-1.00) /Alpha >> 376 * std::exp(AveLogAlpha)/(std::exp(AveLogAlpha)-1.); //ok 432 return tau; 377 return tau; 433 } 378 } 434 379 435 // ------------------------------------------- 380 // ------------------------------------------------------------ 436 381 437 void GFlashSamplingShowerParameterisation::Com << 382 void GFlashSamplingShowerParameterisation:: >> 383 ComputeRadialParameters(G4double Energy, G4double Tau) 438 { 384 { 439 G4double z1 = ParRC1 + ParRC2 * std::log(Ene << 385 G4double z1 = ParRC1 + ParRC2* std::log(Energy/GeV); //ok 440 G4double z2 = ParRC3 + ParRC4 * Zeff; // ok << 386 G4double z2 = ParRC3+ParRC4*Zeff; //ok 441 RadiusCore = z1 + z2 * Tau; // ok << 387 RadiusCore = z1 + z2 * Tau; //ok 442 G4double p1 = ParWC1 + ParWC2 * Zeff; // ok << 388 G4double p1 = ParWC1+ParWC2*Zeff; //ok 443 G4double p2 = ParWC3 + ParWC4 * Zeff; // ok << 389 G4double p2 = ParWC3+ParWC4*Zeff; //ok 444 G4double p3 = ParWC5 + ParWC6 * std::log(Ene << 390 G4double p3 = ParWC5+ParWC6*std::log(Energy/GeV); //ok 445 WeightCore = p1 * std::exp((p2 - Tau) / p3 - << 391 WeightCore = p1 * std::exp( (p2-Tau)/p3- std::exp( (p2-Tau) /p3) ); //ok 446 << 392 447 G4double k1 = ParRT1 + ParRT2 * Zeff; // ok << 393 G4double k1 = ParRT1+ParRT2*Zeff; // ok 448 G4double k2 = ParRT3; // ok << 394 G4double k2 = ParRT3; // ok 449 G4double k3 = ParRT4; // ok << 395 G4double k3 = ParRT4; // ok 450 G4double k4 = ParRT5 + ParRT6 * std::log(Ene << 396 G4double k4 = ParRT5+ParRT6* std::log(Energy/GeV); // ok 451 << 397 452 RadiusTail = k1 * (std::exp(k3 * (Tau - k2)) << 398 RadiusTail = k1*(std::exp(k3*(Tau-k2)) 453 << 399 + std::exp(k4*(Tau-k2)) ); //ok 454 // sampling calorimeter << 400 455 << 401 // sampling calorimeter 456 RadiusCore = RadiusCore + ParsRC1 * (1 - eha << 402 457 WeightCore = << 403 RadiusCore = RadiusCore + ParsRC1*(1-ehat) + ParsRC2/Fs*std::exp(-Tau); //ok 458 WeightCore + (1 - ehat) * (ParsWC1 + ParsW << 404 WeightCore = WeightCore + (1-ehat) 459 RadiusTail = RadiusTail + (1 - ehat) * ParsR << 405 * (ParsWC1+ParsWC2/Fs * std::exp(-std::pow((Tau-1.),2))); //ok >> 406 RadiusTail = RadiusTail + (1-ehat)* ParsRT1+ ParsRT2/Fs *std::exp(-Tau); //ok 460 } 407 } 461 408 462 // ------------------------------------------- 409 // ------------------------------------------------------------ 463 410 464 G4double GFlashSamplingShowerParameterisation: 411 G4double GFlashSamplingShowerParameterisation:: 465 GenerateExponential(const G4double /* Energy * 412 GenerateExponential(const G4double /* Energy */ ) 466 { 413 { 467 G4double ParExp1 = 9./7.*X0eff; 414 G4double ParExp1 = 9./7.*X0eff; 468 G4double random = -ParExp1*G4RandExponentia << 415 G4double random = -ParExp1*CLHEP::RandExponential::shoot() ; 469 return random; 416 return random; 470 } 417 } 471 418