CMSSW Reference Manual

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

Definition at line 26 of file ECalSD.cc.

References birk1, birk2, birk3, birkCut, birkSlope, g, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), initMap(), CaloSD::kmaxIon, CaloSD::kmaxNeutron, CaloSD::kmaxProton, LogDebug, setNumberingScheme(), slopeLY, CaloSD::suppressHeavy, useBirk, useBirkL3, and useWeight.

00028                                                                           : 
00029   CaloSD(name, cpv, clg, p, manager), numberingScheme(0) {
00030   
00031   //   static SimpleConfigurable<bool>   on1(false,  "ECalSD:UseBirkLaw");
00032   //   static SimpleConfigurable<double> bk1(0.00463,"ECalSD:BirkC1");
00033   //   static SimpleConfigurable<double> bk2(-0.03,  "ECalSD:BirkC2");
00034   //   static SimpleConfigurable<double> bk3(1.0,    "ECalSD:BirkC3");
00035   // Values from NIM A484 (2002) 239-244: as implemented in Geant3
00036   //   useBirk          = on1.value();
00037   //   birk1            = bk1.value()*(g/(MeV*cm2));
00038   //   birk2            = bk2.value()*(g/(MeV*cm2))*(g/(MeV*cm2));
00039 
00040   edm::ParameterSet m_EC = p.getParameter<edm::ParameterSet>("ECalSD");
00041   useBirk    = m_EC.getParameter<bool>("UseBirkLaw");
00042   useBirkL3  = m_EC.getParameter<bool>("BirkL3Parametrization");
00043   birk1      = m_EC.getParameter<double>("BirkC1")*(g/(MeV*cm2));
00044   birk2      = m_EC.getParameter<double>("BirkC2");
00045   birk3      = m_EC.getParameter<double>("BirkC3");
00046   birkSlope  = m_EC.getParameter<double>("BirkSlope");
00047   birkCut    = m_EC.getParameter<double>("BirkCut");
00048   slopeLY    = m_EC.getParameter<double>("SlopeLightYield");
00049   bool isItTB= m_EC.getUntrackedParameter<bool>("TestBeam", false);
00050   useWeight= true;
00051 
00052   EcalNumberingScheme* scheme=0;
00053   if      (name == "EcalHitsEB") scheme = dynamic_cast<EcalNumberingScheme*>(new EcalBarrelNumberingScheme());
00054   else if (name == "EcalHitsEE") scheme = dynamic_cast<EcalNumberingScheme*>(new EcalEndcapNumberingScheme());
00055   else if (name == "EcalHitsES") {
00056     if (isItTB) scheme = dynamic_cast<EcalNumberingScheme*>(new ESTBNumberingScheme());
00057     else        scheme = dynamic_cast<EcalNumberingScheme*>(new EcalPreshowerNumberingScheme());
00058     useWeight = false;
00059   } else {edm::LogWarning("EcalSim") << "ECalSD: ReadoutName not supported\n";}
00060 
00061   if (scheme)  setNumberingScheme(scheme);
00062   LogDebug("EcalSim") 
00063     << "Constructing a ECalSD  with name " << GetName() << "\n";
00064   edm::LogInfo("EcalSim")  << "ECalSD:: Use of Birks law is set to      " 
00065                            << useBirk << "        with three constants kB = "
00066                            << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3
00067                            << "\n         Use of L3 parametrization " 
00068                            << useBirkL3 << " with slope " << birkSlope
00069                            << " and cut off " << birkCut << "\n"
00070                            << "         Slope for Light yield is set to "
00071                            << slopeLY;
00072 
00073   edm::LogInfo("EcalSim") << "ECalSD:: Suppression Flag " << suppressHeavy
00074                           << " protons below " << kmaxProton << " MeV,"
00075                           << " neutrons below " << kmaxNeutron << " MeV and"
00076                           << " ions below " << kmaxIon << " MeV";
00077 
00078   if (useWeight) initMap(name,cpv);
00079 
00080 }

ECalSD::~ECalSD ( ) [virtual]

Definition at line 82 of file ECalSD.cc.

References numberingScheme.

00082                 {
00083   if (numberingScheme) delete numberingScheme;
00084 }

double ECalSD::crystalLength ( G4LogicalVolume * lv ) [private]

Definition at line 247 of file ECalSD.cc.

References xtalLMap.

Referenced by curve_LY().

00247                                                 {
00248 
00249   double length= 230.;
00250   std::map<G4LogicalVolume*,double>::const_iterator ite = xtalLMap.find(lv);
00251   if (ite != xtalLMap.end()) length = ite->second;
00252   return length;
00253 }

double ECalSD::curve_LY ( G4Step * aStep ) [private]

Definition at line 218 of file ECalSD.cc.

References crystalLength(), LogDebug, lv, CaloSD::setToLocal(), slopeLY, and weight.

Referenced by getEnergyDeposit().

00218                                      {
00219 
00220   G4StepPoint*     stepPoint = aStep->GetPreStepPoint();
00221   G4LogicalVolume* lv        = stepPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();
00222   G4String         nameVolume= lv->GetName();
00223 
00224   double weight = 1.;
00225   G4ThreeVector  localPoint = setToLocal(stepPoint->GetPosition(),
00226                                          stepPoint->GetTouchable());
00227   double crlength = crystalLength(lv);
00228   double dapd = 0.5 * crlength - localPoint.z();
00229   if (dapd >= -0.1 || dapd <= crlength+0.1) {
00230     if (dapd <= 100.)
00231       weight = 1.0 + slopeLY - dapd * 0.01 * slopeLY;
00232   } else {
00233     edm::LogWarning("EcalSim") << "ECalSD: light coll curve : wrong distance "
00234                                << "to APD " << dapd << " crlength = " 
00235                                << crlength << " crystal name = " << nameVolume 
00236                                << " z of localPoint = " << localPoint.z() 
00237                                << " take weight = " << weight;
00238   }
00239   LogDebug("EcalSim") << "ECalSD, light coll curve : " << dapd 
00240                       << " crlength = " << crlength
00241                       << " crystal name = " << nameVolume 
00242                       << " z of localPoint = " << localPoint.z() 
00243                       << " take weight = " << weight;
00244   return weight;
00245 }

void ECalSD::getBaseNumber ( const G4Step * aStep ) [private]

Definition at line 255 of file ECalSD.cc.

References EcalBaseNumber::addLevel(), EcalBaseNumber::getCapacity(), LogDebug, EcalBaseNumber::reset(), EcalBaseNumber::setSize(), and theBaseNumber.

Referenced by setDetUnitId().

00255                                               {
00256 
00257   theBaseNumber.reset();
00258   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
00259   int theSize = touch->GetHistoryDepth()+1;
00260   if ( theBaseNumber.getCapacity() < theSize ) theBaseNumber.setSize(theSize);
00261   //Get name and copy numbers
00262   if ( theSize > 1 ) {
00263     for (int ii = 0; ii < theSize ; ii++) {
00264       theBaseNumber.addLevel(touch->GetVolume(ii)->GetName(),touch->GetReplicaNumber(ii));
00265       LogDebug("EcalSim") << "ECalSD::getBaseNumber(): Adding level " << ii
00266                           << ": " << touch->GetVolume(ii)->GetName() << "["
00267                           << touch->GetReplicaNumber(ii) << "]";
00268     }
00269   }
00270 
00271 }

double ECalSD::getBirkL3 ( G4Step * aStep ) [private]

Definition at line 273 of file ECalSD.cc.

References birk1, birkCut, birkSlope, funct::log(), LogDebug, and weight.

Referenced by getEnergyDeposit().

00273                                       {
00274 
00275   double weight = 1.;
00276   double charge = aStep->GetPreStepPoint()->GetCharge();
00277 
00278   if (charge != 0. && aStep->GetStepLength() > 0) {
00279     G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
00280     double density = mat->GetDensity();
00281     double dedx    = aStep->GetTotalEnergyDeposit()/aStep->GetStepLength();
00282     double rkb     = birk1/density;
00283     weight         = 1. - birkSlope*log(rkb*dedx);
00284     if (weight < birkCut) weight = birkCut;
00285     else if (weight > 1.) weight = 1.;
00286     LogDebug("EcalSim") << "ECalSD::getBirkL3 in " << mat->GetName()
00287                         << " Charge " << charge << " dE/dx " << dedx
00288                         << " Birk Const " << rkb << " Weight = " << weight 
00289                         << " dE " << aStep->GetTotalEnergyDeposit();
00290   }
00291   return weight;
00292 
00293 }

double ECalSD::getEnergyDeposit ( G4Step * aStep ) [virtual]

Reimplemented from CaloSD.

Definition at line 86 of file ECalSD.cc.

References birk1, birk2, birk3, curve_LY(), CaloSD::getAttenuation(), getBirkL3(), TrackInformation::isPrimary(), CaloSD::kmaxIon, CaloSD::kmaxNeutron, CaloSD::kmaxProton, LogDebug, NULL, CaloSD::preStepPoint, CaloSD::suppressHeavy, CaloSD::theTrack, useBirk, useBirkL3, useWeight, and weight.

00086                                               {
00087   
00088   if (aStep == NULL) {
00089     return 0;
00090   } else {
00091     preStepPoint        = aStep->GetPreStepPoint();
00092     G4String nameVolume = preStepPoint->GetPhysicalVolume()->GetName();
00093 
00094     // take into account light collection curve for crystals
00095     double weight = 1.;
00096     if (suppressHeavy) {
00097       G4Track* theTrack = aStep->GetTrack();
00098       TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
00099       if (trkInfo) {
00100         int pdg = theTrack->GetDefinition()->GetPDGEncoding();
00101         if (!(trkInfo->isPrimary())) { // Only secondary particles
00102           double ke = theTrack->GetKineticEnergy()/MeV;
00103           if (((pdg/1000000000 == 1 && ((pdg/10000)%100) > 0 &&
00104                 ((pdg/10)%100) > 0)) && (ke<kmaxIon)) weight = 0;
00105           if ((pdg == 2212) && (ke < kmaxProton))     weight = 0;
00106           if ((pdg == 2112) && (ke < kmaxNeutron))    weight = 0;
00107           if (weight == 0) {
00108             LogDebug("EcalSim") << "Ignore Track " << theTrack->GetTrackID()
00109                                 << " Type " << theTrack->GetDefinition()->GetParticleName()
00110                                 << " Kinetic Energy " << ke << " MeV";
00111           }
00112         }
00113       }
00114     }
00115     if (useWeight) {
00116       weight *= curve_LY(aStep);
00117       if (useBirk) {
00118         if (useBirkL3) weight *= getBirkL3(aStep);
00119         else           weight *= getAttenuation(aStep, birk1, birk2, birk3);
00120       }
00121     }
00122     double edep   = aStep->GetTotalEnergyDeposit() * weight;
00123     LogDebug("EcalSim") << "ECalSD:: " << nameVolume
00124                         <<" Light Collection Efficiency " << weight 
00125                         << " Weighted Energy Deposit " << edep/MeV << " MeV";
00126     return edep;
00127   } 
00128 }

int ECalSD::getRadiationLenght ( G4Step * aStep ) [virtual]

Definition at line 130 of file ECalSD.cc.

References int, SensitiveDetector::name, NULL, r, funct::sqrt(), x, y, and z.

00130                                              {
00131   
00132   if (aStep == NULL) {
00133     return 0;
00134   } else {
00135     
00136     G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
00137     G4VPhysicalVolume* currentPV  = aStep->GetPreStepPoint()->GetPhysicalVolume();
00138     G4String name = currentPV->GetName();
00139     std::string crystal;
00140     crystal.assign(name,0,4);
00141     
00142     int thisX0 = 0;
00143     if (crystal == "EFRY"){
00144       float z = hitPoint.z();
00145       float detz = fabs(fabs(z)-3200);
00146       thisX0 = (int)floor( detz/8.9 );   
00147     } 
00148     if(crystal == "EBRY") {
00149       float x = hitPoint.x();
00150       float y = hitPoint.y();
00151       float r = sqrt(x*x +y*y);
00152       float detr = r -1290;
00153       thisX0 = (int)floor( detr/8.9);
00154     }
00155     
00156     return thisX0;
00157   }
00158 }

void ECalSD::initMap	(	G4String	sd,
		const DDCompactView &	cpv
	)			`[private]`

Definition at line 174 of file ECalSD.cc.

References DDFilteredView::addFilter(), DDSplit(), ddtrap, DDSpecificsFilter::equals, filter, first, DDFilteredView::firstChild(), i, LogDebug, DDFilteredView::logicalPart(), lv, SensitiveDetector::name, DDBase< N, C >::name(), DDFilteredView::next(), DDSolid::parameters(), DDSpecificsFilter::setCriteria(), DDSolid::shape(), DDLogicalPart::solid(), and xtalLMap.

Referenced by ECalSD().

00174                                                            {
00175 
00176   G4String attribute = "ReadOutName";
00177   DDSpecificsFilter filter;
00178   DDValue           ddv(attribute,sd,0);
00179   filter.setCriteria(ddv,DDSpecificsFilter::equals);
00180   DDFilteredView fv(cpv);
00181   fv.addFilter(filter);
00182   fv.firstChild();
00183 
00184   const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
00185   std::vector<G4LogicalVolume *>::const_iterator lvcite;
00186   bool dodet=true;
00187   while (dodet) {
00188     const DDSolid & sol  = fv.logicalPart().solid();
00189     const std::vector<double> & paras = sol.parameters();
00190     G4String name = DDSplit(sol.name()).first;
00191     G4LogicalVolume* lv=0;
00192     for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) 
00193       if ((*lvcite)->GetName() == name) {
00194         lv = (*lvcite);
00195         break;
00196       }
00197     LogDebug("EcalSim") << "ECalSD::initMap (for " << sd << "): Solid " << name
00198                         << " Shape " << sol.shape() << " Parameter 0 = " 
00199                         << paras[0] << " Logical Volume " << lv;
00200     if (sol.shape() == ddtrap) {
00201       double dz = 2*paras[0];
00202       xtalLMap.insert(std::pair<G4LogicalVolume*,double>(lv,dz));
00203     }
00204     dodet = fv.next();
00205   }
00206   LogDebug("EcalSim") << "ECalSD: Length Table for " << attribute << " = " 
00207                       << sd << ":";   
00208   std::map<G4LogicalVolume*,double>::const_iterator ite = xtalLMap.begin();
00209   int i=0;
00210   for (; ite != xtalLMap.end(); ite++, i++) {
00211     G4String name = "Unknown";
00212     if (ite->first != 0) name = (ite->first)->GetName();
00213     LogDebug("EcalSim") << " " << i << " " << ite->first << " " << name 
00214                         << " L = " << ite->second;
00215   }
00216 }

uint32_t ECalSD::setDetUnitId ( G4Step * aStep ) [virtual]

Implements CaloSD.

Definition at line 160 of file ECalSD.cc.

References getBaseNumber(), EcalNumberingScheme::getUnitID(), numberingScheme, and theBaseNumber.

00160                                             { 
00161   getBaseNumber(aStep);
00162   return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(theBaseNumber));
00163 }

void ECalSD::setNumberingScheme ( EcalNumberingScheme * scheme )

Definition at line 165 of file ECalSD.cc.

References numberingScheme.

Referenced by ECalSD(), and HcalTB04Analysis::update().

00165                                                            {
00166   if (scheme != 0) {
00167     edm::LogInfo("EcalSim") << "EcalSD: updates numbering scheme for " 
00168                             << GetName() << "\n";
00169     if (numberingScheme) delete numberingScheme;
00170     numberingScheme = scheme;
00171   }
00172 }

double ECalSD::birk1 [private]