CMSSW Reference Manual

ZdcSD::ZdcSD	(	G4String	name,
		const DDCompactView &	cpv,
		SensitiveDetectorCatalog &	clg,
		edm::ParameterSet const &	p,
		const SimTrackManager *	manager
	)

Definition at line 25 of file ZdcSD.cc.

References edm::ParameterSet::getParameter(), setNumberingScheme(), and verbosity.

00027                                                                        : 
00028   CaloSD(name, cpv, clg, p, manager), numberingScheme(0) {
00029   edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
00030   verbosity  = m_ZdcSD.getParameter<int>("Verbosity");
00031   int verbn  = verbosity/10;
00032   verbosity %= 10;
00033   setNumberingScheme(new ZdcNumberingScheme(verbn));
00034 
00035   edm::LogInfo("ForwardSim")
00036     << "***************************************************\n"
00037     << "*                                                 *\n"
00038     << "* Constructing a ZdcSD  with name " << name << "\n"
00039     << "*                                                 *\n"
00040     << "***************************************************";
00041 }

ZdcSD::~ZdcSD ( ) [virtual]

Definition at line 43 of file ZdcSD.cc.

00043 {}

double ZdcSD::getEnergyDeposit	(	G4Step *	aStep,
		edm::ParameterSet const &	p
	)			`[virtual]`

Definition at line 45 of file ZdcSD.cc.

References a, DeDxTools::beta(), d, eta, Exception, edm::ParameterSet::getParameter(), edm::Service< T >::isAvailable(), funct::log(), LogDebug, max, min, NULL, phi, pi, CaloSD::preStepPoint, r, funct::sqrt(), StDecayID::status, funct::tan(), theta, CaloSD::theTrack, and thFibDir.

00045                                                                          {
00046 
00047   float NCherPhot = 0.;
00048 
00049   if (aStep == NULL) {
00050     LogDebug("ForwardSim") << "ZdcSD::  getEnergyDeposit: aStep is NULL!";
00051     return 0;
00052   }
00053   else {
00054     // preStepPoint information
00055     G4SteppingControl  stepControlFlag = aStep->GetControlFlag();
00056     G4StepPoint*       preStepPoint = aStep->GetPreStepPoint();
00057     G4VPhysicalVolume* currentPV    = preStepPoint->GetPhysicalVolume();
00058     G4String           nameVolume   = currentPV->GetName();
00059 
00060     G4ThreeVector      hitPoint = preStepPoint->GetPosition();  
00061     G4ThreeVector      hit_mom = preStepPoint->GetMomentumDirection();
00062     G4double           stepL = aStep->GetStepLength()/cm;
00063     G4double           beta     = preStepPoint->GetBeta();
00064     G4double           charge   = preStepPoint->GetCharge();
00065 
00066     // G4VProcess*        curprocess   = preStepPoint->GetProcessDefinedStep();
00067     // G4String           namePr   = preStepPoint->GetProcessDefinedStep()->GetProcessName();
00068     // G4LogicalVolume*   lv    = currentPV->GetLogicalVolume();
00069     // G4Material*        mat   = lv->GetMaterial();
00070     // G4double           rad   = mat->GetRadlen();
00071 
00072     // postStepPoint information
00073     G4StepPoint* postStepPoint = aStep->GetPostStepPoint();   
00074     G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
00075     G4String postnameVolume = postPV->GetName();
00076 
00077     // theTrack information
00078     G4Track* theTrack = aStep->GetTrack();   
00079     G4String particleType = theTrack->GetDefinition()->GetParticleName();
00080     G4int primaryID = theTrack->GetTrackID();
00081     G4double entot = theTrack->GetTotalEnergy();
00082     G4ThreeVector vert_mom = theTrack->GetVertexMomentumDirection();
00083     G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
00084 
00085     // calculations
00086     float costheta = vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+
00087                                        vert_mom.y()*vert_mom.y()+
00088                                        vert_mom.z()*vert_mom.z());
00089     float theta = acos(std::min(std::max(costheta,float(-1.)),float(1.)));
00090     float eta = -log(tan(theta/2));
00091     float phi = -100.;
00092     if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x()); 
00093     if (phi < 0.) phi += twopi;
00094 
00095     // Get the total energy deposit
00096     double stepE   = aStep->GetTotalEnergyDeposit();
00097     LogDebug("ForwardSim") 
00098       << "ZdcSD::  getEnergyDeposit: "
00099       <<"*****************HHHHHHHHHHHHHHHHHHHHHHHHHHLLLLLLLLLlllllllllll&&&&&&&&&&\n"
00100       << "  preStepPoint: " << nameVolume << "," << stepL << "," << stepE 
00101       << "," << beta << "," << charge << "\n"
00102       << "  postStepPoint: " << postnameVolume << "," << costheta << "," 
00103       << theta << "," << eta << "," << phi << "," << particleType << "," 
00104       << primaryID;
00105 
00106     float bThreshold = 0.67;
00107     int status = 0;
00108     if ((beta > bThreshold) && (charge != 0) && (nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber")) {
00109       status = 1;
00110       LogDebug("ForwardSim") << "ZdcSD::  getEnergyDeposit:  pass "; 
00111 
00112       float nMedium = 1.4925;
00113       // float photEnSpectrDL = 10714.285714;
00114       //       photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
00115 
00116       float photEnSpectrDE = 1.24;
00117       // E = 2pi*(1./137.)*(eV*cm/370.)/lambda = 12.389184*(eV*cm)/lambda
00118       // Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm  = 3.01 eV
00119       // Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm  = 1.77 eV
00120       // delE = Emax - Emin = 1.24 eV
00121 
00122       float effPMTandTransport = 0.15;
00123 
00124       // Check these values
00125       float thFullRefl = 23.;
00126       float thFullReflRad = thFullRefl*pi/180.;
00127 
00128       edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
00129       thFibDir  = m_ZdcSD.getParameter<double>("FiberDirection");
00130       //float thFibDir = 90.;
00131       float thFibDirRad = thFibDir*pi/180.;
00132 
00133       // at which theta the point is located:
00134       //   float th1 = hitPoint.theta();
00135 
00136       // theta of charged particle in LabRF(hit momentum direction):
00137       float costh = hit_mom.z()/sqrt(hit_mom.x()*hit_mom.x()+
00138                                      hit_mom.y()*hit_mom.y()+
00139                                      hit_mom.z()*hit_mom.z());
00140       float th = acos(std::min(std::max(costh,float(-1.)),float(1.)));
00141       // just in case (can do both standard ranges of phi):
00142       if (th < 0.) th += twopi;
00143 
00144       // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
00145       float costhcher =1./(nMedium*beta);
00146       float thcher = acos(std::min(std::max(costhcher,float(-1.)),float(1.)));
00147 
00148       // diff thetas of charged part. and quartz direction in LabRF:
00149       float DelFibPart = fabs(th - thFibDirRad);
00150 
00151       // define real distances:
00152       float d = fabs(tan(th)-tan(thFibDirRad));   
00153 
00154       // float a = fabs(tan(thFibDirRad)-tan(thFibDirRad+thFullReflRad));   
00155       // float r = fabs(tan(th)-tan(th+thcher));   
00156       float a = tan(thFibDirRad)+tan(fabs(thFibDirRad-thFullReflRad));   
00157       float r = tan(th)+tan(fabs(th-thcher));   
00158       
00159       // std::cout.testOut << "  d=|tan(" << th << ")-tan(" << thFibDirRad << ")| "
00160       //              << "=|" << tan(th) << "-" << tan(thFibDirRad) << "| = " << d;
00161       // std::cout.testOut << "  a=tan(" << thFibDirRad << ")=" << tan(thFibDirRad) 
00162       //              << " + tan(|" << thFibDirRad << " - " << thFullReflRad << "|)="
00163       //              << tan(fabs(thFibDirRad-thFullReflRad)) << " = " << a;
00164       // std::cout.testOut << "  r=tan(" << th << ")=" << tan(th) << " + tan(|" << th 
00165       //              << " - " << thcher << "|)=" << tan(fabs(th-thcher)) << " = " << r;
00166 
00167       // define losses d_qz in cone of full reflection inside quartz direction
00168       float d_qz = -1;
00169       float variant = -1;
00170 
00171       // if (d > (r+a))
00172       if (DelFibPart > (thFullReflRad + thcher) ) {
00173         variant = 0.; d_qz = 0.;
00174       }
00175       else {
00176         // if ((DelFibPart + thcher) < thFullReflRad )  [(d+r) < a]
00177         if ((th + thcher) < (thFibDirRad+thFullReflRad) && (th - thcher) > (thFibDirRad-thFullReflRad) ) {
00178           variant = 1.; d_qz = 1.;
00179         }
00180         else {
00181           // if ((thcher - DelFibPart ) > thFullReflRad )  [(r-d) > a]
00182           if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher) ) {
00183             variant = 2.; d_qz = 0.;
00184           }
00185           else {
00186             // if ((thcher + DelFibPart ) > thFullReflRad && thcher < (DelFibPart+thFullReflRad) ) {  [(r+d) > a && (r-d) < a)]
00187             variant = 3.; // d_qz is calculated below
00188 
00189             // use crossed length of circles(cone projection) - dC1/dC2 : 
00190             float arg_arcos = 0.;
00191             float tan_arcos = 2.*a*d;
00192             if (tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos; 
00193             // std::cout.testOut << "  d_qz: " << r << "," << a << "," << d << " " << tan_arcos << " " << arg_arcos;
00194             arg_arcos = fabs(arg_arcos);
00195             // std::cout.testOut << "," << arg_arcos;
00196             float th_arcos = acos(std::min(std::max(arg_arcos,float(-1.)),float(1.)));
00197             // std::cout.testOut << " " << th_arcos;
00198             d_qz = th_arcos/pi/2.;
00199             // std::cout.testOut << " " << d_qz;
00200             d_qz = fabs(d_qz);
00201             // std::cout.testOut << "," << d_qz;
00202           }
00203         }
00204       }
00205 
00206       //  std::cout<< std::endl;
00207 
00208       double meanNCherPhot = 0.;
00209       G4int poissNCherPhot = 0;
00210       if (d_qz > 0) {
00211         meanNCherPhot = 370.*charge*charge*( 1. - 1./(nMedium*nMedium*beta*beta) ) * photEnSpectrDE * stepL;
00212 
00213         // dLamdX:  meanNCherPhot = (2.*pi/137.)*charge*charge* 
00214         //                          ( 1. - 1./(nMedium*nMedium*beta*beta) ) * photEnSpectrDL * stepL;
00215         edm::Service<edm::RandomNumberGenerator> rng;
00216         if ( ! rng.isAvailable()) {
00217           throw cms::Exception("Configuration")
00218             << "ZdcSD requires the RandomNumberGeneratorService\n"
00219             << "which is not present in the configuration file.  "
00220             << "You must add the service\n" << "in the configuration file "
00221             << "or remove the modules that require it.";
00222         }
00223         CLHEP::RandPoissonQ randPoisson(rng->getEngine());
00224         poissNCherPhot = (G4int) randPoisson.fire(meanNCherPhot);
00225         // poissNCherPhot = (G4int) G4Poisson(meanNCherPhot);
00226 
00227         if (poissNCherPhot < 0) poissNCherPhot = 0; 
00228 
00229         // NCherPhot = meanNCherPhot;
00230         NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
00231       }
00232 
00233       LogDebug("ForwardSim") 
00234         << "ZdcSD::  getEnergyDeposit:  gED: "
00235         << stepE
00236         << "," << costh
00237         << "," << th
00238         << "," << costhcher
00239         << "," << thcher
00240         << "," << DelFibPart
00241         << "," << d
00242         << "," << a
00243         << "," << r
00244         << "," << hitPoint
00245         << "," << hit_mom
00246         << "," << stepControlFlag
00247         << "," << entot
00248         << "," << vert_mom
00249         << "," << localPoint
00250         << "," << charge
00251         << "," << beta
00252         << "," << stepL
00253         << "," << d_qz
00254         << "," << variant
00255         << "," << meanNCherPhot
00256         << "," << poissNCherPhot
00257         << "," << NCherPhot;
00258       // --constants-----------------
00259       // << "," << photEnSpectrDE
00260       // << "," << nMedium
00261       // << "," << bThreshold
00262       // << "," << thFibDirRad
00263       // << "," << thFullReflRad
00264       // << "," << effPMTandTransport
00265       // --other variables-----------
00266       // << "," << curprocess
00267       // << "," << nameProcess
00268       // << "," << name
00269       // << "," << rad
00270       // << "," << mat
00271 
00272     }
00273     else {
00274       // determine failure mode: beta, charge, and/or nameVolume
00275       status = 0;
00276       if (beta <= bThreshold)
00277         LogDebug("ForwardSim") 
00278           << "ZdcSD::  getEnergyDeposit: fail beta=" << beta;
00279       if (charge == 0)
00280         LogDebug("ForwardSim") 
00281           << "ZdcSD::  getEnergyDeposit: fail charge=0";
00282       if ( !(nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber") )
00283         LogDebug("ForwardSim") 
00284           << "ZdcSD::  getEnergyDeposit: fail nv=" << nameVolume;
00285     }
00286 
00287 
00288     return NCherPhot;
00289   } 
00290 }

uint32_t ZdcSD::setDetUnitId ( G4Step * step ) [virtual]

Implements CaloSD.

Definition at line 292 of file ZdcSD.cc.

References ZdcNumberingScheme::getUnitID(), and numberingScheme.

00292                                           {
00293   uint32_t returnNumber = 0;
00294   if(numberingScheme != 0)returnNumber = numberingScheme->getUnitID(aStep);
00295   // edm: return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(aStep));
00296   return returnNumber;
00297 }

void ZdcSD::setNumberingScheme ( ZdcNumberingScheme * scheme )

Definition at line 299 of file ZdcSD.cc.

References numberingScheme.

Referenced by ZdcSD().

00299                                                          {
00300   if (scheme != 0) {
00301     edm::LogInfo("ForwardSim") << "ZdcSD: updates numbering scheme for " 
00302                                << GetName();
00303     if (numberingScheme) delete numberingScheme;
00304     numberingScheme = scheme;
00305   }
00306 }

ZdcNumberingScheme* ZdcSD::numberingScheme [private]

Definition at line 26 of file ZdcSD.h.

Referenced by setDetUnitId(), and setNumberingScheme().

double ZdcSD::thFibDir [private]

Definition at line 25 of file ZdcSD.h.

Referenced by getEnergyDeposit().

int ZdcSD::verbosity [private]

Definition at line 24 of file ZdcSD.h.

Referenced by ZdcSD().


Public Member Functions
virtual double	getEnergyDeposit (G4Step *, edm::ParameterSet const &)
virtual uint32_t	setDetUnitId (G4Step *step)
void	setNumberingScheme (ZdcNumberingScheme *scheme)
	ZdcSD (G4String, const DDCompactView &, SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *)
virtual	~ZdcSD ()
Private Attributes
ZdcNumberingScheme *	numberingScheme
double	thFibDir
int	verbosity

ZdcSD Class Reference

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation