Reference Manual

Go to the documentation of this file.
00001 #include "SimG4CMS/Calo/interface/HFGflash.h"
00002 
00003 #include "G4VPhysicalVolume.hh"
00004 #include "G4Step.hh"
00005 #include "G4Track.hh"
00006 #include "G4Navigator.hh"
00007 #include "G4NavigationHistory.hh"
00008 #include "CLHEP/Units/PhysicalConstants.h"
00009 #include "CLHEP/Units/SystemOfUnits.h"
00010 
00011 #include "Randomize.hh"
00012 #include "G4TransportationManager.hh"
00013 #include "G4VPhysicalVolume.hh"
00014 #include "G4LogicalVolume.hh"
00015 #include "G4VSensitiveDetector.hh"
00016 #include "G4EventManager.hh"
00017 #include "G4SteppingManager.hh"
00018 #include "G4FastTrack.hh"
00019 #include "G4ParticleTable.hh"
00020 
00021 #include "CLHEP/GenericFunctions/IncompleteGamma.hh"
00022 
00023 #include "SimG4Core/Application/interface/SteppingAction.h"
00024 #include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
00025 #include "SimGeneral/GFlash/interface/GflashTrajectoryPoint.h"
00026 
00027 #include "FWCore/ServiceRegistry/interface/Service.h"
00028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00029 
00030 #include <math.h>
00031 
00032 //#define DebugLog
00033 
00034 HFGflash::HFGflash(edm::ParameterSet const & p) {
00035 
00036   edm::ParameterSet m_HF  = p.getParameter<edm::ParameterSet>("HFGflash");
00037   theBField               = m_HF.getUntrackedParameter<double>("BField", 3.8);
00038   theWatcherOn            = m_HF.getUntrackedParameter<bool>("WatcherOn",true);
00039   theFillHisto            = m_HF.getUntrackedParameter<bool>("FillHisto",true);
00040   edm::LogInfo("HFShower") << "HFGFlash:: Set B-Field to " << theBField
00041                            << ", WatcherOn to " << theWatcherOn
00042                            << " and FillHisto to " << theFillHisto;
00043 
00044   theHelix = new GflashTrajectory;
00045   theGflashStep = new G4Step();
00046   theGflashNavigator = new G4Navigator();
00047   //  theGflashNavigator = 0;
00048   theGflashTouchableHandle = new G4TouchableHistory();
00049 
00050 #ifdef DebugLog
00051   if (theFillHisto) {
00052     edm::Service<TFileService> tfile;
00053     if ( tfile.isAvailable() ) {
00054       TFileDirectory showerDir = tfile->mkdir("GflashEMShowerProfile");
00055 
00056       em_incE = showerDir.make<TH1F>("em_incE","Incoming energy (GeV)",500,0,500.);
00057       em_ssp_rho    = showerDir.make<TH1F>("em_ssp_rho","Shower starting position;#rho (cm);Number of Events",100,100.0,200.0);
00058       em_ssp_z      = showerDir.make<TH1F>("em_ssp_z","Shower starting position;z (cm);Number of Events",2000,0.0,2000.0);
00059       em_long       = showerDir.make<TH1F>("em_long","Longitudinal Profile;Radiation Length;Number of Spots",800,800.0,1600.0);
00060       em_lateral    = showerDir.make<TH1F>("em_lateral","Lateral Profile;Radiation Length;Moliere Radius",100,0.0,5.0);
00061       em_2d         = showerDir.make<TH2F>("em_2d","Lateral Profile vs. Shower Depth;Radiation Length;Moliere Radius",800,800.0,1600.0,100,0.0,5.0);
00062       em_long_sd    = showerDir.make<TH1F>("em_long_sd","Longitudinal Profile in Sensitive Detector;Radiation Length;Number of Spots",800,800.0,1600.0);
00063       em_lateral_sd = showerDir.make<TH1F>("em_lateral_sd","Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",100,0.0,5.0);
00064       em_2d_sd      = showerDir.make<TH2F>("em_2d_sd","Lateral Profile vs. Shower Depth in Sensitive Detector;Radiation Length;Moliere Radius",800,800.0,1600.0,100,0.0,5.0);
00065       em_ze         = showerDir.make<TH2F>("em_ze","Profile vs. Energy;Radiation Length;Moliere Radius",800,800.0,1600.0,1000,0.0,1.0);
00066       em_ratio      = showerDir.make<TH2F>("em_ratio","Profile vs. Energy;Radiation Length;Moliere Radius",800,800.0,1600.0,1000,0.0,100.0);
00067       em_ratio_selected      = showerDir.make<TH2F>("em_ratio_selected","Profile vs. Energy;Radiation Length;Moliere Radius",800,800.0,1600.0,1000,0.0,100.0);
00068       em_nSpots_sd  = showerDir.make<TH1F>("em_nSpots_sd","Number of Gflash Spots in Sensitive Detector;Number of Spots;Number of Events",100,0.0,100);
00069       em_ze_ratio   = showerDir.make<TH1F>("em_ze_ratio","Ratio of Energy and Z Position",1000,0.0,0.001);
00070     } else {
00071       theFillHisto = false;
00072       edm::LogInfo("HFShower") << "HFGFlash::No file is available for saving"
00073                                << " histos so the flag is set to false";
00074     }
00075   }
00076 #endif
00077   jCalorimeter = Gflash::kHF;
00078 
00079 }
00080 
00081 HFGflash::~HFGflash() {
00082   if (theHelix)      delete theHelix;
00083   if (theGflashStep) delete theGflashStep;
00084   if (theGflashNavigator) delete theGflashNavigator;
00085 }
00086 
00087 
00088 std::vector<HFGflash::Hit> HFGflash::gfParameterization(G4Step * aStep,bool & ok,bool onlyLong) {
00089   double tempZCalo = 26;  // Gflash::Z[jCalorimeter]
00090   double hfcriticalEnergy = 0.021;  // Gflash::criticalEnergy
00091 
00092   std::vector<HFGflash::Hit> hit;
00093   HFGflash::Hit oneHit;
00094 
00095   G4StepPoint * preStepPoint  = aStep->GetPreStepPoint(); 
00096   //G4StepPoint * postStepPoint = aStep->GetPostStepPoint(); 
00097   G4Track *     track    = aStep->GetTrack();
00098   // Get Z-direction 
00099   const G4DynamicParticle *aParticle = track->GetDynamicParticle();
00100   G4ThreeVector momDir = aParticle->GetMomentumDirection();
00101 
00102   G4ThreeVector hitPoint = preStepPoint->GetPosition();   
00103   G4String      partType = track->GetDefinition()->GetParticleName();
00104   //  int           parCode  = track->GetDefinition()->GetPDGEncoding();
00105 
00106   // This part of code is copied from the original GFlash Fortran code.
00107   // reference : hep-ex/0001020v1
00108 
00109   const G4double energyCutoff     = 1; 
00110   const G4int    maxNumberOfSpots = 10000000;
00111 
00112   G4ThreeVector showerStartingPosition = track->GetPosition()/cm;
00113   G4ThreeVector showerMomentum = preStepPoint->GetMomentum()/GeV;
00114   jCalorimeter = Gflash::kHF;
00115 
00116   G4double logEinc = std::log((preStepPoint->GetTotalEnergy())/GeV);
00117 
00118   G4double y = ((preStepPoint->GetTotalEnergy())/GeV) / hfcriticalEnergy; // y = E/Ec, criticalEnergy is in GeV
00119   G4double logY = std::log(y);
00120 
00121 
00122   G4double nSpots = 93.0 * std::log(tempZCalo) * ((preStepPoint->GetTotalEnergy())/GeV); // total number of spot due linearization
00123   if(preStepPoint->GetTotalEnergy()/GeV < 1.6)  nSpots = 140.4 * std::log(tempZCalo) * std::pow(((preStepPoint->GetTotalEnergy())/GeV),0.876); 
00124 
00125 
00126   //   // implementing magnetic field effects
00127   double charge = track->GetStep()->GetPreStepPoint()->GetCharge();
00128   theHelix->initializeTrajectory(showerMomentum,showerStartingPosition,charge,theBField);
00129 
00130   G4double pathLength0 = theHelix->getPathLengthAtZ(showerStartingPosition.getZ());
00131   G4double pathLength = pathLength0; // this will grow along the shower development
00132 
00133   //--- 2.2  Fix intrinsic properties of em. showers.
00134 
00135   G4double fluctuatedTmax = std::log(logY - 0.7157);
00136   G4double fluctuatedAlpha= std::log(0.7996 +(0.4581 + 1.8628/tempZCalo)*logY);
00137 
00138   G4double sigmaTmax = 1.0/( -1.4  + 1.26 * logY);
00139   G4double sigmaAlpha = 1.0/( -0.58 + 0.86 * logY);
00140   G4double rho = 0.705  - 0.023 * logY;
00141   G4double sqrtPL = std::sqrt((1.0+rho)/2.0);
00142   G4double sqrtLE = std::sqrt((1.0-rho)/2.0);
00143 
00144   G4double norm1 = G4RandGauss::shoot();
00145   G4double norm2 = G4RandGauss::shoot();
00146   G4double tempTmax = fluctuatedTmax + sigmaTmax*(sqrtPL*norm1 + sqrtLE*norm2);
00147   G4double tempAlpha = fluctuatedAlpha + sigmaAlpha*(sqrtPL*norm1 - sqrtLE*norm2);
00148 
00149   // tmax, alpha, beta : parameters of gamma distribution
00150   G4double tmax = std::exp(tempTmax);
00151   G4double alpha = std::exp(tempAlpha);
00152   G4double beta = (alpha - 1.0)/tmax;
00153 
00154   if (!alpha)          return hit; 
00155   if (!beta)           return hit;
00156   if (alpha < 0.00001) return hit;
00157   if (beta < 0.00001)  return hit;
00158  
00159   // spot fluctuations are added to tmax, alpha, beta
00160   G4double averageTmax = logY-0.858;
00161   G4double averageAlpha = 0.21+(0.492+2.38/tempZCalo)*logY;
00162   G4double spotTmax  = averageTmax * (0.698 + .00212*tempZCalo);
00163   G4double spotAlpha= averageAlpha * (0.639 + .00334*tempZCalo);
00164   G4double spotBeta = (spotAlpha-1.0)/spotTmax;
00165 
00166   if (!spotAlpha)          return hit;
00167   if (!spotBeta)           return hit;
00168   if (spotAlpha < 0.00001) return hit;
00169   if (spotBeta < 0.00001)  return hit;
00170 
00171 #ifdef DebugLog  
00172   LogDebug("HFShower") << "Incoming energy = " << ((preStepPoint->GetTotalEnergy())/GeV) << " Position (rho,z) = (" << showerStartingPosition.rho() << ", " << showerStartingPosition.z() << ")";
00173 
00174   if(theFillHisto) {
00175     em_incE->Fill(((preStepPoint->GetTotalEnergy())/GeV));
00176     em_ssp_rho->Fill(showerStartingPosition.rho());
00177     em_ssp_z->Fill(std::abs(showerStartingPosition.z()));
00178   }
00179 #endif
00180   //  parameters for lateral distribution and fluctuation
00181   G4double z1=0.0251+0.00319*logEinc;
00182   G4double z2=0.1162-0.000381*tempZCalo;
00183 
00184   G4double k1=0.659 - 0.00309 * tempZCalo;
00185   G4double k2=0.645;
00186   G4double k3=-2.59;
00187   G4double k4=0.3585+ 0.0421*logEinc;
00188 
00189   G4double p1=2.623 -0.00094*tempZCalo;
00190   G4double p2=0.401 +0.00187*tempZCalo;
00191   G4double p3=1.313 -0.0686*logEinc;
00192 
00193   //   // @@@ dwjang, intial tuning by comparing 20-150GeV TB data
00194   //   // the width of energy response is not yet tuned.
00195   G4double e25Scale = 1.03551;
00196   z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790e+00*std::log(((preStepPoint->GetTotalEnergy())/GeV)) - 3.66237e+00);
00197   p1 *= 0.96;
00198 
00199   G4double stepLengthLeft = 10000;
00200   G4int    nSpots_sd = 0; // count total number of spots in SD
00201   G4double zInX0 = 0.0; // shower depth in X0 unit
00202   G4double deltaZInX0 = 0.0; // segment of depth in X0 unit
00203   G4double deltaZ = 0.0; // segment of depth in cm
00204   G4double stepLengthLeftInX0 = 0.0; // step length left in X0 unit
00205 
00206   const G4double divisionStepInX0 = 0.1; //step size in X0 unit
00207   G4double energy = ((preStepPoint->GetTotalEnergy())/GeV); // energy left in GeV
00208 
00209   Genfun::IncompleteGamma gammaDist;
00210 
00211   G4double energyInGamma = 0.0; // energy in a specific depth(z) according to Gamma distribution
00212   G4double preEnergyInGamma = 0.0; // energy calculated in a previous depth
00213   G4double sigmaInGamma  = 0.; // sigma of energy in a specific depth(z) according to Gamma distribution
00214   G4double preSigmaInGamma = 0.0; // sigma of energy in a previous depth
00215 
00216   //energy segment in Gamma distribution of shower in each step  
00217   G4double deltaEnergy =0.0 ; // energy in deltaZ
00218   G4int spotCounter = 0; // keep track of number of spots generated
00219 
00220   //step increment along the shower direction
00221   G4double deltaStep = 0.0;
00222 
00223   // Uniqueness of G4Step is important otherwise hits won't be created.
00224   G4double timeGlobal = track->GetStep()->GetPreStepPoint()->GetGlobalTime();
00225 
00226   // this needs to be deleted manually at the end of this loop.
00227   //  theGflashNavigator = new G4Navigator();
00228   theGflashNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
00229 
00230   //   // loop for longitudinal integration
00231 
00232 #ifdef DebugLog  
00233   LogDebug("HFShower") << " Energy = " << energy << " Step Length Left = "  << stepLengthLeft;
00234 #endif
00235   while(energy > 0.0 && stepLengthLeft > 0.0) { 
00236     stepLengthLeftInX0 = stepLengthLeft / Gflash::radLength[jCalorimeter];
00237 
00238     if ( stepLengthLeftInX0 < divisionStepInX0 ) {
00239       deltaZInX0 = stepLengthLeftInX0;
00240       deltaZ     = deltaZInX0 * Gflash::radLength[jCalorimeter];
00241       stepLengthLeft = 0.0;
00242     } else {
00243       deltaZInX0 = divisionStepInX0;
00244       deltaZ     = deltaZInX0 * Gflash::radLength[jCalorimeter];
00245       stepLengthLeft -= deltaZ;
00246     }
00247 
00248     zInX0 += deltaZInX0;
00249     
00250 #ifdef DebugLog  
00251     LogDebug("HFShower") << " zInX0 = " << zInX0 << " spotBeta*zInX0 = " << spotBeta*zInX0;
00252 #endif
00253     if ((!zInX0) || (!spotBeta*zInX0) || (zInX0 < 0.01) || 
00254         (spotBeta*zInX0 < 0.00001) || (!zInX0*beta) || (zInX0*beta < 0.00001)) 
00255       return hit;
00256 
00257     G4int nSpotsInStep = 0;
00258 
00259 #ifdef DebugLog  
00260     LogDebug("HFShower") << " Energy - Energy Cut off = " << energy - energyCutoff;
00261 #endif
00262 
00263     if ( energy > energyCutoff  ) {
00264       preEnergyInGamma  = energyInGamma;
00265       gammaDist.a().setValue(alpha);  //alpha
00266       
00267       energyInGamma = gammaDist(beta*zInX0); //beta
00268       G4double energyInDeltaZ  = energyInGamma - preEnergyInGamma;
00269       deltaEnergy   = std::min(energy,((preStepPoint->GetTotalEnergy())/GeV)*energyInDeltaZ);
00270  
00271       preSigmaInGamma  = sigmaInGamma;
00272       gammaDist.a().setValue(spotAlpha);  //alpha spot
00273       sigmaInGamma = gammaDist(spotBeta*zInX0); //beta spot
00274       nSpotsInStep = std::max(1,int(nSpots * (sigmaInGamma - preSigmaInGamma)));
00275     } else {
00276       deltaEnergy = energy;
00277       preSigmaInGamma  = sigmaInGamma;
00278       nSpotsInStep = std::max(1,int(nSpots * (1.0 - preSigmaInGamma)));
00279     }
00280 
00281     if ( deltaEnergy > energy || (energy-deltaEnergy) < energyCutoff ) deltaEnergy = energy;
00282 
00283     energy  -= deltaEnergy;
00284 
00285     if ( spotCounter+nSpotsInStep > maxNumberOfSpots ) {
00286       nSpotsInStep = maxNumberOfSpots - spotCounter;
00287       if (nSpotsInStep < 1) nSpotsInStep = 1;
00288     }
00289 
00290 
00291     //     // It begins with 0.5 of deltaZ and then icreases by 1 deltaZ
00292     deltaStep  += 0.5*deltaZ;
00293     pathLength += deltaStep;
00294     deltaStep   =  0.5*deltaZ;
00295 
00296 
00297     //lateral shape and fluctuations for  homogenous calo.
00298     G4double tScale = tmax *alpha/(alpha-1.0) * (std::exp(fluctuatedAlpha)-1.0)/std::exp(fluctuatedAlpha);
00299     G4double tau = std::min(10.0,(zInX0 - 0.5*deltaZInX0)/tScale);
00300     G4double rCore = z1 + z2 * tau; 
00301     G4double rTail = k1 *( std::exp(k3*(tau-k2)) + std::exp(k4*(tau-k2))); // @@ check RT3 sign
00302     G4double p23 = (p2 - tau)/p3;
00303     G4double probabilityWeight = p1 *  std::exp( p23 - std::exp(p23) );
00304 
00305 
00306     // Deposition of spots according to lateral distr.
00307     // Apply absolute energy scale
00308     // Convert into MeV unit
00309     G4double emSpotEnergy = deltaEnergy / nSpotsInStep * e25Scale * GeV;
00310 
00311 
00312 #ifdef DebugLog  
00313     LogDebug("HFShower") << " nSpotsInStep = " << nSpotsInStep;
00314 #endif
00315 
00316 
00317 
00318     for (G4int ispot = 0 ;  ispot < nSpotsInStep ; ispot++) {
00319       spotCounter++;
00320       G4double u1 = G4UniformRand();
00321       G4double u2 = G4UniformRand();
00322       G4double rInRM = 0.0;
00323       
00324       if (u1 < probabilityWeight) {
00325         rInRM = rCore * std::sqrt( u2/(1.0-u2) );
00326       } else {
00327         rInRM = rTail * std::sqrt( u2/(1.0-u2) );
00328       }
00329   
00330       G4double rShower =  rInRM * Gflash::rMoliere[jCalorimeter];
00331       
00332       //Uniform & random rotation of spot along the azimuthal angle
00333       G4double azimuthalAngle = twopi*G4UniformRand();
00334       
00335       //Compute global position of generated spots with taking into account magnetic field
00336       //Divide deltaZ into nSpotsInStep and give a spot a global position
00337       G4double incrementPath = (deltaZ/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
00338 
00339       // trajectoryPoint give a spot an imaginary point along the shower development
00340       GflashTrajectoryPoint trajectoryPoint;
00341       theHelix->getGflashTrajectoryPoint(trajectoryPoint,pathLength+incrementPath);
00342 
00343 
00344       G4ThreeVector SpotPosition0 = trajectoryPoint.getPosition() + rShower*std::cos(azimuthalAngle)*trajectoryPoint.getOrthogonalUnitVector() + rShower*std::sin(azimuthalAngle)*trajectoryPoint.getCrossUnitVector();
00345 
00346       // Convert into mm unit
00347       SpotPosition0 *= cm;
00348 
00349 
00350       //---------------------------------------------------
00351       // fill a fake step to send it to hit maker
00352       //---------------------------------------------------
00353 
00354       // to make a different time for each fake step. (0.03 nsec is corresponding to 1cm step size)
00355       timeGlobal += 0.0001*nanosecond;
00356 
00357       //fill equivalent changes to a (fake) step associated with a spot 
00358       
00359       G4double zInX0_spot = std::abs(pathLength+incrementPath - pathLength0)/Gflash::radLength[jCalorimeter];
00360 
00361 #ifdef DebugLog  
00362       LogDebug("HFShower") <<  "zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot << " , " << emSpotEnergy/GeV <<  "emSpotEnergy/GeV =" << emSpotEnergy/GeV;
00363 #endif
00364 
00365       if ((!zInX0_spot) || (zInX0_spot < 0)) continue;
00366       if ((!emSpotEnergy/GeV) ||  (emSpotEnergy < 0)) continue;
00367       if ((!rShower/Gflash::rMoliere[jCalorimeter]) || (rShower/Gflash::rMoliere[jCalorimeter] < 0)) continue;
00368 
00369       oneHit.depth    = 1;
00370 
00371 #ifdef DebugLog
00372       if (theFillHisto) {
00373         em_long->Fill(SpotPosition0.z()/cm,emSpotEnergy/GeV);
00374         em_lateral->Fill(rShower/Gflash::rMoliere[jCalorimeter],emSpotEnergy/GeV);
00375         em_2d->Fill(SpotPosition0.z()/cm,rShower/Gflash::rMoliere[jCalorimeter],emSpotEnergy/GeV);
00376       }
00377 #endif
00378 
00379       if(SpotPosition0 == 0) continue;
00380 
00381       double energyratio = emSpotEnergy/(preStepPoint->GetTotalEnergy()/(nSpots*e25Scale));
00382 
00383       if (emSpotEnergy/GeV < 0.0001) continue;
00384       if (energyratio > 80) continue;
00385 
00386       double zshift =0;
00387       if(SpotPosition0.z() > 0) zshift = 18;
00388       if(SpotPosition0.z() < 0) zshift = -18;
00389 
00390 
00391       G4ThreeVector gfshift(0,0,zshift*(pow(100,0.1)/pow(preStepPoint->GetTotalEnergy()/GeV,0.1)));
00392 
00393       G4ThreeVector SpotPosition = gfshift + SpotPosition0;
00394 
00395       double LengthWeight = std::fabs(std::pow(SpotPosition0.z()/11370,1));
00396       if (G4UniformRand()>  0.0021 * energyratio * LengthWeight) continue;
00397 
00398       oneHit.position = SpotPosition;
00399       oneHit.time     = timeGlobal;
00400       oneHit.edep     = emSpotEnergy/GeV;
00401       hit.push_back(oneHit);
00402       nSpots_sd++;
00403       
00404     } // end of for spot iteration
00405 
00406   } // end of while for longitudinal integration
00407 #ifdef DebugLog
00408   if (theFillHisto) {
00409     em_nSpots_sd->Fill(nSpots_sd);
00410   }
00411 #endif
00412   //  delete theGflashNavigator;
00413   return hit;
00414 }
/data/refman/pasoursint/CMSSW_5_3_3/src/SimG4CMS/Calo/src/HFGflash.cc
