00001
00002
00003
00005
00006 #include "SimG4CMS/Calo/interface/HCalSD.h"
00007 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
00008 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
00009 #include "SimG4Core/Notification/interface/TrackInformation.h"
00010 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00011 #include "DetectorDescription/Core/interface/DDFilter.h"
00012 #include "DetectorDescription/Core/interface/DDFilteredView.h"
00013 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
00014 #include "DetectorDescription/Core/interface/DDMaterial.h"
00015 #include "DetectorDescription/Core/interface/DDValue.h"
00016 #include "FWCore/Utilities/interface/Exception.h"
00017
00018 #include "FWCore/ServiceRegistry/interface/Service.h"
00019 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00020
00021 #include "G4LogicalVolumeStore.hh"
00022 #include "G4LogicalVolume.hh"
00023 #include "G4Step.hh"
00024 #include "G4Track.hh"
00025 #include "G4ParticleTable.hh"
00026 #include "G4VProcess.hh"
00027
00028 #include <iostream>
00029 #include <fstream>
00030 #include <iomanip>
00031
00032
00033
00034
00035 HCalSD::HCalSD(G4String name, const DDCompactView & cpv,
00036 SensitiveDetectorCatalog & clg,
00037 edm::ParameterSet const & p, const SimTrackManager* manager) :
00038 CaloSD(name, cpv, clg, p, manager,
00039 p.getParameter<edm::ParameterSet>("HCalSD").getParameter<int>("TimeSliceUnit"),
00040 p.getParameter<edm::ParameterSet>("HCalSD").getParameter<bool>("IgnoreTrackID")),
00041 numberingFromDDD(0), numberingScheme(0), showerLibrary(0), hfshower(0),
00042 showerParam(0), showerPMT(0), showerBundle(0), m_HEDarkening(0),
00043 m_HFDarkening(0) {
00044
00045
00046
00047
00048
00049
00050
00051 edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("HCalSD");
00052 useBirk = m_HC.getParameter<bool>("UseBirkLaw");
00053 birk1 = m_HC.getParameter<double>("BirkC1")*(g/(MeV*cm2));
00054 birk2 = m_HC.getParameter<double>("BirkC2");
00055 birk3 = m_HC.getParameter<double>("BirkC3");
00056 useShowerLibrary = m_HC.getParameter<bool>("UseShowerLibrary");
00057 useParam = m_HC.getParameter<bool>("UseParametrize");
00058 bool testNumber = m_HC.getParameter<bool>("TestNumberingScheme");
00059 usePMTHit = m_HC.getParameter<bool>("UsePMTHits");
00060 betaThr = m_HC.getParameter<double>("BetaThreshold");
00061 eminHitHB = m_HC.getParameter<double>("EminHitHB")*MeV;
00062 eminHitHE = m_HC.getParameter<double>("EminHitHE")*MeV;
00063 eminHitHO = m_HC.getParameter<double>("EminHitHO")*MeV;
00064 eminHitHF = m_HC.getParameter<double>("EminHitHF")*MeV;
00065 useFibreBundle = m_HC.getParameter<bool>("UseFibreBundleHits");
00066 deliveredLumi = m_HC.getParameter<double>("DelivLuminosity");
00067 bool ageingFlagHE= m_HC.getParameter<bool>("HEDarkening");
00068 bool ageingFlagHF= m_HC.getParameter<bool>("HFDarkening");
00069 useHF = m_HC.getUntrackedParameter<bool>("UseHF",true);
00070 bool forTBH2 = m_HC.getUntrackedParameter<bool>("ForTBH2",false);
00071 useLayerWt = m_HC.getUntrackedParameter<bool>("UseLayerWt",false);
00072 std::string file = m_HC.getUntrackedParameter<std::string>("WtFile","None");
00073 edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
00074 applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
00075
00076 #ifdef DebugLog
00077 LogDebug("HcalSim") << "***************************************************"
00078 << "\n"
00079 << "* *"
00080 << "\n"
00081 << "* Constructing a HCalSD with name " << name << "\n"
00082 << "* *"
00083 << "\n"
00084 << "***************************************************";
00085 #endif
00086 edm::LogInfo("HcalSim") << "HCalSD:: Use of HF code is set to " << useHF
00087 << "\nUse of shower parametrization set to "
00088 << useParam << "\nUse of shower library is set to "
00089 << useShowerLibrary << "\nUse PMT Hit is set to "
00090 << usePMTHit << " with beta Threshold "<< betaThr
00091 << "\nUSe of FibreBundle Hit set to "<<useFibreBundle
00092 << "\n Use of Birks law is set to "
00093 << useBirk << " with three constants kB = "
00094 << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3;
00095 edm::LogInfo("HcalSim") << "HCalSD:: Suppression Flag " << suppressHeavy
00096 << " protons below " << kmaxProton << " MeV,"
00097 << " neutrons below " << kmaxNeutron << " MeV and"
00098 << " ions below " << kmaxIon << " MeV\n"
00099 << " Threshold for storing hits in HB: "
00100 << eminHitHB << " HE: " << eminHitHE << " HO: "
00101 << eminHitHO << " HF: " << eminHitHF << "\n"
00102 << "Delivered luminosity for Darkening "
00103 << deliveredLumi << " Flag (HE) " << ageingFlagHE
00104 << " Flag (HF) " << ageingFlagHF << "\n"
00105 << "Application of Fiducial Cut " << applyFidCut;
00106
00107 numberingFromDDD = new HcalNumberingFromDDD(name, cpv);
00108 HcalNumberingScheme* scheme;
00109 if (testNumber || forTBH2)
00110 scheme = dynamic_cast<HcalNumberingScheme*>(new HcalTestNumberingScheme(forTBH2));
00111 else
00112 scheme = new HcalNumberingScheme();
00113 setNumberingScheme(scheme);
00114
00115 const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
00116 std::vector<G4LogicalVolume *>::const_iterator lvcite;
00117 G4LogicalVolume* lv;
00118 std::string attribute, value;
00119 if (useHF) {
00120 if (useParam) {
00121 showerParam = new HFShowerParam(name, cpv, p);
00122 } else {
00123 if (useShowerLibrary) showerLibrary = new HFShowerLibrary(name, cpv, p);
00124 hfshower = new HFShower(name, cpv, p, 0);
00125 }
00126
00127
00128 attribute = "Volume";
00129 value = "HF";
00130 DDSpecificsFilter filter0;
00131 DDValue ddv0(attribute, value, 0);
00132 filter0.setCriteria(ddv0, DDSpecificsFilter::equals);
00133 DDFilteredView fv0(cpv);
00134 fv0.addFilter(filter0);
00135 hfNames = getNames(fv0);
00136 fv0.firstChild();
00137 DDsvalues_type sv0(fv0.mergedSpecifics());
00138 std::vector<double> temp = getDDDArray("Levels",sv0);
00139 edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
00140 << " = " << value << " has " << hfNames.size()
00141 << " elements";
00142 for (unsigned int i=0; i < hfNames.size(); ++i) {
00143 G4String namv = hfNames[i];
00144 lv = 0;
00145 for(lvcite=lvs->begin(); lvcite!=lvs->end(); lvcite++)
00146 if((*lvcite)->GetName()==namv) {
00147 lv = (*lvcite);
00148 break;
00149 }
00150 hfLV.push_back(lv);
00151 int level = static_cast<int>(temp[i]);
00152 hfLevels.push_back(level);
00153 edm::LogInfo("HcalSim") << "HCalSD: HF[" << i << "] = " << hfNames[i]
00154 << " LV " << hfLV[i] << " at level "
00155 << hfLevels[i];
00156 }
00157
00158
00159 value = "HFFibre";
00160 DDSpecificsFilter filter1;
00161 DDValue ddv1(attribute,value,0);
00162 filter1.setCriteria(ddv1, DDSpecificsFilter::equals);
00163 DDFilteredView fv1(cpv);
00164 fv1.addFilter(filter1);
00165 fibreNames = getNames(fv1);
00166 edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
00167 << " = " << value << ":";
00168 for (unsigned int i=0; i<fibreNames.size(); ++i) {
00169 G4String namv = fibreNames[i];
00170 lv = 0;
00171 for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite) {
00172 if ((*lvcite)->GetName() == namv) {
00173 lv = (*lvcite);
00174 break;
00175 }
00176 }
00177 fibreLV.push_back(lv);
00178 edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
00179 << " LV " << fibreLV[i];
00180 }
00181
00182
00183 value = "HFPMT";
00184 DDSpecificsFilter filter3;
00185 DDValue ddv3(attribute,value,0);
00186 filter3.setCriteria(ddv3,DDSpecificsFilter::equals);
00187 DDFilteredView fv3(cpv);
00188 fv3.addFilter(filter3);
00189 std::vector<G4String> pmtNames = getNames(fv3);
00190 edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
00191 << " = " << value << " have " << pmtNames.size()
00192 << " entries";
00193 for (unsigned int i=0; i<pmtNames.size(); ++i) {
00194 G4String namv = pmtNames[i];
00195 lv = 0;
00196 for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite)
00197 if ((*lvcite)->GetName() == namv) {
00198 lv = (*lvcite);
00199 break;
00200 }
00201 pmtLV.push_back(lv);
00202 edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << pmtNames[i]
00203 << " LV " << pmtLV[i];
00204 }
00205 if (pmtNames.size() > 0) showerPMT = new HFShowerPMT (name, cpv, p);
00206
00207
00208 value = "HFFibreBundleStraight";
00209 DDSpecificsFilter filter4;
00210 DDValue ddv4(attribute,value,0);
00211 filter4.setCriteria(ddv4,DDSpecificsFilter::equals);
00212 DDFilteredView fv4(cpv);
00213 fv4.addFilter(filter4);
00214 std::vector<G4String> fibreNames = getNames(fv4);
00215 edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
00216 << " = " << value << " have " << fibreNames.size()
00217 << " entries";
00218 for (unsigned int i=0; i<fibreNames.size(); ++i) {
00219 G4String namv = fibreNames[i];
00220 lv = 0;
00221 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
00222 if ((*lvcite)->GetName() == namv) {
00223 lv = (*lvcite);
00224 break;
00225 }
00226 fibre1LV.push_back(lv);
00227 edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
00228 << " LV " << fibre1LV[i];
00229 }
00230
00231
00232 value = "HFFibreBundleConical";
00233 DDSpecificsFilter filter5;
00234 DDValue ddv5(attribute,value,0);
00235 filter5.setCriteria(ddv5,DDSpecificsFilter::equals);
00236 DDFilteredView fv5(cpv);
00237 fv5.addFilter(filter5);
00238 fibreNames = getNames(fv5);
00239 edm::LogInfo("HcalSim") << "HCalSD: Names to be tested for " << attribute
00240 << " = " << value << " have " << fibreNames.size()
00241 << " entries";
00242 for (unsigned int i=0; i<fibreNames.size(); ++i) {
00243 G4String namv = fibreNames[i];
00244 lv = 0;
00245 for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite)
00246 if ((*lvcite)->GetName() == namv) {
00247 lv = (*lvcite);
00248 break;
00249 }
00250 fibre2LV.push_back(lv);
00251 edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << fibreNames[i]
00252 << " LV " << fibre2LV[i];
00253 }
00254 if (fibre1LV.size() > 0 || fibre2LV.size() > 0)
00255 showerBundle = new HFShowerFibreBundle (name, cpv, p);
00256
00257 attribute = "ReadOutName";
00258 value = name;
00259 DDSpecificsFilter filter6;
00260 DDValue ddv6(attribute,value,0);
00261 filter6.setCriteria(ddv6,DDSpecificsFilter::equals);
00262 DDFilteredView fv6(cpv);
00263 fv6.addFilter(filter6);
00264 if (fv6.firstChild()) {
00265 DDsvalues_type sv(fv6.mergedSpecifics());
00266
00267 gpar = getDDDArray("gparHF",sv);
00268 edm::LogInfo("HFShower") << "HFShowerParam: " << gpar.size()
00269 << " gpar (cm)";
00270 for (unsigned int ig=0; ig<gpar.size(); ig++)
00271 edm::LogInfo("HFShower") << "HFShowerParam: gpar[" << ig << "] = "
00272 << gpar[ig]/cm << " cm";
00273 } else {
00274 edm::LogWarning("HFShower") << "HFShowerParam: cannot get filtered "
00275 << " view for " << attribute << " matching "
00276 << name;
00277 }
00278 }
00279
00280
00281 attribute = "ReadOutName";
00282 DDSpecificsFilter filter2;
00283 DDValue ddv2(attribute,name,0);
00284 filter2.setCriteria(ddv2,DDSpecificsFilter::equals);
00285 DDFilteredView fv2(cpv);
00286 fv2.addFilter(filter2);
00287 bool dodet = fv2.firstChild();
00288
00289 DDsvalues_type sv(fv2.mergedSpecifics());
00290
00291 layer0wt = getDDDArray("Layer0Wt",sv);
00292 edm::LogInfo("HcalSim") << "HCalSD: " << layer0wt.size() << " Layer0Wt";
00293 for (unsigned int it=0; it<layer0wt.size(); ++it)
00294 edm::LogInfo("HcalSim") << "HCalSD: [" << it << "] " << layer0wt[it];
00295
00296 const G4MaterialTable * matTab = G4Material::GetMaterialTable();
00297 std::vector<G4Material*>::const_iterator matite;
00298 while (dodet) {
00299 const DDLogicalPart & log = fv2.logicalPart();
00300 G4String namx = log.name().name();
00301 if (!isItHF(namx) && !isItFibre(namx)) {
00302 bool notIn = true;
00303 for (unsigned int i=0; i<matNames.size(); ++i) {
00304 if (!strcmp(matNames[i].c_str(),log.material().name().name().c_str())){
00305 notIn = false;
00306 break;
00307 }
00308 }
00309 if (notIn) {
00310 namx = log.material().name().name();
00311 matNames.push_back(namx);
00312 G4Material* mat;
00313 for (matite = matTab->begin(); matite != matTab->end(); ++matite) {
00314 if ((*matite)->GetName() == namx) {
00315 mat = (*matite);
00316 break;
00317 }
00318 }
00319 materials.push_back(mat);
00320 }
00321 }
00322 dodet = fv2.next();
00323 }
00324
00325 edm::LogInfo("HcalSim") << "HCalSD: Material names for " << attribute
00326 << " = " << name << ":";
00327 for (unsigned int i=0; i<matNames.size(); ++i)
00328 edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << matNames[i]
00329 << " pointer " << materials[i];
00330
00331 mumPDG = mupPDG = 0;
00332
00333 if (useLayerWt) readWeightFromFile(file);
00334
00335 for (int i=0; i<9; ++i) hit_[i] = time_[i]= dist_[i] = 0;
00336 hzvem = hzvhad = 0;
00337
00338 if (ageingFlagHE) m_HEDarkening = new HEDarkening();
00339 if (ageingFlagHF) m_HFDarkening = new HFDarkening();
00340 #ifdef plotDebug
00341 edm::Service<TFileService> tfile;
00342
00343 if ( tfile.isAvailable() ) {
00344 static std::string labels[9] = {"HB", "HE", "HO", "HF Absorber", "HF PMT",
00345 "HF Absorber Long", "HF Absorber Short",
00346 "HF PMT Long", "HF PMT Short"};
00347 TFileDirectory hcDir = tfile->mkdir("ProfileFromHCalSD");
00348 char name[20], title[60];
00349 for (int i=0; i<9; ++i) {
00350 sprintf (title, "Hit energy in %s", labels[i].c_str());
00351 sprintf (name, "HCalSDHit%d", i);
00352 hit_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
00353 sprintf (title, "Energy (MeV)");
00354 hit_[i]->GetXaxis()->SetTitle(title);
00355 hit_[i]->GetYaxis()->SetTitle("Hits");
00356 sprintf (title, "Time of the hit in %s", labels[i].c_str());
00357 sprintf (name, "HCalSDTime%d", i);
00358 time_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
00359 sprintf (title, "Time (ns)");
00360 time_[i]->GetXaxis()->SetTitle(title);
00361 time_[i]->GetYaxis()->SetTitle("Hits");
00362 sprintf (title, "Longitudinal profile in %s", labels[i].c_str());
00363 sprintf (name, "HCalSDDist%d", i);
00364 dist_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
00365 sprintf (title, "Distance (mm)");
00366 dist_[i]->GetXaxis()->SetTitle(title);
00367 dist_[i]->GetYaxis()->SetTitle("Hits");
00368 }
00369 if (useHF && (!useParam)) {
00370 hzvem = hcDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part)",330,0.0,1650.0);
00371 hzvem->GetXaxis()->SetTitle("Longitudinal Profile (EM Part)");
00372 hzvhad = hcDir.make<TH1F>("hzvhad","Longitudinal Profile (Had Part)",330,0.0,1650.0);
00373 hzvhad->GetXaxis()->SetTitle("Longitudinal Profile (Hadronic Part)");
00374 }
00375 }
00376 #endif
00377 }
00378
00379 HCalSD::~HCalSD() {
00380
00381 if (numberingFromDDD) delete numberingFromDDD;
00382 if (numberingScheme) delete numberingScheme;
00383 if (showerLibrary) delete showerLibrary;
00384 if (hfshower) delete hfshower;
00385 if (showerParam) delete showerParam;
00386 if (showerPMT) delete showerPMT;
00387 if (showerBundle) delete showerBundle;
00388 if (m_HEDarkening) delete m_HEDarkening;
00389 if (m_HFDarkening) delete m_HFDarkening;
00390 }
00391
00392 bool HCalSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
00393
00394 NaNTrap( aStep ) ;
00395
00396 if (aStep == NULL) {
00397 return true;
00398 } else {
00399 G4LogicalVolume* lv =
00400 aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
00401 G4String nameVolume = lv->GetName();
00402 if (isItHF(aStep)) {
00403 G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
00404 double weight(1.0);
00405 if (m_HFDarkening != 0) {
00406 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
00407 double r = hitPoint.perp()/CLHEP::cm;
00408 double z = std::abs(hitPoint.z())/CLHEP::cm;
00409 float dose_acquired = 0.;
00410 if (z>=1100 && z <= 1300) {
00411 int hfZLayer = (int)((z - 1100)/20);
00412 if (hfZLayer > 9) hfZLayer = 9;
00413 float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
00414 for (int i = hfZLayer; i <= 9; ++i) {
00415 dose_acquired = m_HFDarkening->dose(i,r);
00416 weight *= m_HFDarkening->degradation(normalized_lumi*dose_acquired);
00417 }
00418 }
00419 #ifdef DebugLog
00420 LogDebug("HcalSim") << "HCalSD: HFLumiDarkening at r = " << r
00421 << ", z = " << z << " Dose " << dose_acquired
00422 << " weight " << weight;
00423 #endif
00424 }
00425 if (useParam) {
00426 #ifdef DebugLog
00427 LogDebug("HcalSim") << "HCalSD: " << getNumberOfHits()
00428 << " hits from parametrization in " << nameVolume
00429 << " for Track " << aStep->GetTrack()->GetTrackID()
00430 <<" (" << aStep->GetTrack()->GetDefinition()->GetParticleName()
00431 <<")";
00432 #endif
00433 getFromParam(aStep, weight);
00434 #ifdef DebugLog
00435 LogDebug("HcalSim") << "HCalSD: " << getNumberOfHits()
00436 << " hits afterParamS*";
00437 #endif
00438 } else {
00439 bool notaMuon = true;
00440 if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
00441 if (useShowerLibrary && notaMuon) {
00442 #ifdef DebugLog
00443 LogDebug("HcalSim") << "HCalSD: Starts shower library from "
00444 << nameVolume << " for Track "
00445 << aStep->GetTrack()->GetTrackID() << " ("
00446 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
00447 #endif
00448 getFromLibrary(aStep, weight);
00449 } else if (isItFibre(lv)) {
00450 #ifdef DebugLog
00451 LogDebug("HcalSim") << "HCalSD: Hit at Fibre in " << nameVolume
00452 << " for Track "
00453 << aStep->GetTrack()->GetTrackID() <<" ("
00454 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
00455 #endif
00456 hitForFibre(aStep, weight);
00457 }
00458 }
00459 } else if (isItPMT(lv)) {
00460 #ifdef DebugLog
00461 LogDebug("HcalSim") << "HCalSD: Hit from PMT parametrization from "
00462 << nameVolume << " for Track "
00463 << aStep->GetTrack()->GetTrackID() << " ("
00464 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
00465 #endif
00466 if (usePMTHit && showerPMT) getHitPMT(aStep);
00467 } else if (isItStraightBundle(lv) || isItConicalBundle(lv)) {
00468 #ifdef DebugLog
00469 LogDebug("HcalSim") << "HCalSD: Hit from FibreBundle from "
00470 << nameVolume << " for Track "
00471 << aStep->GetTrack()->GetTrackID() << " ("
00472 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
00473 #endif
00474 if (useFibreBundle && showerBundle)
00475 getHitFibreBundle(aStep, isItConicalBundle(lv));
00476 } else {
00477 #ifdef DebugLog
00478 LogDebug("HcalSim") << "HCalSD: Hit from standard path from "
00479 << nameVolume << " for Track "
00480 << aStep->GetTrack()->GetTrackID() << " ("
00481 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
00482 #endif
00483 if (getStepInfo(aStep)) {
00484 #ifdef plotDebug
00485 if (edepositEM+edepositHAD > 0)
00486 plotProfile(aStep, aStep->GetPreStepPoint()->GetPosition(),
00487 edepositEM+edepositHAD,aStep->GetPostStepPoint()->GetGlobalTime(),0);
00488 #endif
00489 if (hitExists() == false && edepositEM+edepositHAD>0.) currentHit = createNewHit();
00490 }
00491 }
00492 return true;
00493 }
00494 }
00495
00496 double HCalSD::getEnergyDeposit(G4Step* aStep) {
00497 double destep = aStep->GetTotalEnergyDeposit();
00498 double weight = 1;
00499 G4Track* theTrack = aStep->GetTrack();
00500
00501 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
00502 int depth = (touch->GetReplicaNumber(0))%10;
00503 int det = ((touch->GetReplicaNumber(1))/1000)-3;
00504 if (depth==0 && (det==0 || det==1)) weight = layer0wt[det];
00505 if (useLayerWt) {
00506 int lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
00507 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
00508 weight = layerWeight(det+3, hitPoint, depth, lay);
00509 }
00510
00511 if (m_HEDarkening !=0 && det == 1) {
00512 int lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
00513 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
00514 float r = (hitPoint.perp())/CLHEP::cm;
00515 #ifdef DebugLog
00516 LogDebug("HcalSim") << "HCalSD:HE_Darkening >>> position: "<< hitPoint
00517 << " lay: " << lay << " R: " << r << " cm ";
00518 #endif
00519 float normalized_lumi = m_HEDarkening->int_lumi(deliveredLumi);
00520 float dose_acquired = m_HEDarkening->dose(lay-2,r);
00521 weight *= m_HEDarkening->degradation(normalized_lumi * dose_acquired);
00522 #ifdef DebugLog
00523 LogDebug("HcalSim") << "HCalSD: >>> norm_Lumi: " << normalized_lumi
00524 << " dose: " << dose_acquired
00525 << " coefficient = " << weight;
00526 #endif
00527 }
00528
00529 if (suppressHeavy) {
00530 TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
00531 if (trkInfo) {
00532 int pdg = theTrack->GetDefinition()->GetPDGEncoding();
00533 if (!(trkInfo->isPrimary())) {
00534 double ke = theTrack->GetKineticEnergy()/MeV;
00535 if ( pdg/1000000000 == 1 && (pdg/10000)%100 > 0 &&
00536 (pdg/10)%100 > 0 && ke <kmaxIon ) weight = 0;
00537 if ((pdg == 2212) && (ke < kmaxProton)) weight = 0;
00538 if ((pdg == 2112) && (ke < kmaxNeutron)) weight = 0;
00539 #ifdef DebugLog
00540 if (weight == 0)
00541 edm::LogInfo("HcalSim") << "HCalSD:Ignore Track "
00542 << theTrack->GetTrackID() << " Type "
00543 << theTrack->GetDefinition()->GetParticleName()
00544 << " Kinetic Energy " << ke << " MeV";
00545 #endif
00546 }
00547 }
00548 }
00549 #ifdef DebugLog
00550 double weight0 = weight;
00551 #endif
00552 if (useBirk) {
00553 G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
00554 if (isItScintillator(mat)) weight *= getAttenuation(aStep, birk1, birk2, birk3);
00555 }
00556 double wt1 = getResponseWt(theTrack);
00557 double wt2 = theTrack->GetWeight();
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572 #ifdef DebugLog
00573 edm::LogInfo("HcalSim") << "HCalSD: Detector " << det+3 << " Depth " << depth
00574 << " weight " << weight0 << " " << weight << " " << wt1
00575 << " " << wt2;
00576 #endif
00577 double edep = weight*wt1*destep;
00578 if(wt2 > 0.0) { edep *= wt2; }
00579 return edep;
00580 }
00581
00582 uint32_t HCalSD::setDetUnitId(G4Step * aStep) {
00583
00584 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
00585 const G4VTouchable* touch = preStepPoint->GetTouchable();
00586 G4ThreeVector hitPoint = preStepPoint->GetPosition();
00587
00588 int depth = (touch->GetReplicaNumber(0))%10 + 1;
00589 int lay = (touch->GetReplicaNumber(0)/10)%100 + 1;
00590 int det = (touch->GetReplicaNumber(1))/1000;
00591
00592 return setDetUnitId (det, hitPoint, depth, lay);
00593 }
00594
00595 void HCalSD::setNumberingScheme(HcalNumberingScheme * scheme) {
00596 if (scheme != 0) {
00597 edm::LogInfo("HcalSim") << "HCalSD: updates numbering scheme for " << GetName();
00598 if (numberingScheme) delete numberingScheme;
00599 numberingScheme = scheme;
00600 }
00601 }
00602
00603 void HCalSD::initRun() {
00604 G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
00605 G4String particleName;
00606 mumPDG = theParticleTable->FindParticle(particleName="mu-")->GetPDGEncoding();
00607 mupPDG = theParticleTable->FindParticle(particleName="mu+")->GetPDGEncoding();
00608 #ifdef DebugLog
00609 LogDebug("HcalSim") << "HCalSD: Particle code for mu- = " << mumPDG
00610 << " for mu+ = " << mupPDG;
00611 #endif
00612 if (showerLibrary) showerLibrary->initRun(theParticleTable);
00613 if (showerParam) showerParam->initRun(theParticleTable);
00614 if (hfshower) hfshower->initRun(theParticleTable);
00615 }
00616
00617 bool HCalSD::filterHit(CaloG4Hit* aHit, double time) {
00618 double threshold=0;
00619 DetId theId(aHit->getUnitID());
00620 switch (theId.subdetId()) {
00621 case HcalBarrel:
00622 threshold = eminHitHB; break;
00623 case HcalEndcap:
00624 threshold = eminHitHE; break;
00625 case HcalOuter:
00626 threshold = eminHitHO; break;
00627 case HcalForward:
00628 threshold = eminHitHF; break;
00629 default:
00630 break;
00631 }
00632 return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > threshold));
00633 }
00634
00635
00636 uint32_t HCalSD::setDetUnitId (int det, G4ThreeVector pos, int depth, int lay=1) {
00637 uint32_t id = 0;
00638 if (numberingFromDDD) {
00639
00640 HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, pos, depth, lay);
00641
00642 if (numberingScheme) id = numberingScheme->getUnitID(tmp);
00643 }
00644 return id;
00645 }
00646
00647 std::vector<double> HCalSD::getDDDArray(const std::string & str,
00648 const DDsvalues_type & sv) {
00649 #ifdef DebugLog
00650 LogDebug("HcalSim") << "HCalSD:getDDDArray called for " << str;
00651 #endif
00652 DDValue value(str);
00653 if (DDfetch(&sv,value)) {
00654 #ifdef DebugLog
00655 LogDebug("HcalSim") << value;
00656 #endif
00657 const std::vector<double> & fvec = value.doubles();
00658 int nval = fvec.size();
00659 if (nval < 1) {
00660 edm::LogError("HcalSim") << "HCalSD : # of " << str << " bins " << nval
00661 << " < 2 ==> illegal";
00662 throw cms::Exception("Unknown", "HCalSD") << "nval < 2 for array " << str << "\n";
00663 }
00664
00665 return fvec;
00666 } else {
00667 edm::LogError("HcalSim") << "HCalSD : cannot get array " << str;
00668 throw cms::Exception("Unknown", "HCalSD") << "cannot get array " << str << "\n";
00669 }
00670 }
00671
00672 std::vector<G4String> HCalSD::getNames(DDFilteredView& fv) {
00673
00674 std::vector<G4String> tmp;
00675 bool dodet = fv.firstChild();
00676 while (dodet) {
00677 const DDLogicalPart & log = fv.logicalPart();
00678 bool ok = true;
00679
00680 for (unsigned int i=0; i<tmp.size(); ++i) {
00681 if (!strcmp(tmp[i].c_str(), log.name().name().c_str())) {
00682 ok = false;
00683 break;
00684 }
00685 }
00686 if (ok) tmp.push_back(log.name().name());
00687 dodet = fv.next();
00688 }
00689 return tmp;
00690 }
00691
00692 bool HCalSD::isItHF(G4Step * aStep) {
00693 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
00694 int levels = (touch->GetHistoryDepth()) + 1;
00695 for (unsigned int it=0; it < hfNames.size(); ++it) {
00696 if (levels >= hfLevels[it]) {
00697 G4LogicalVolume* lv = touch->GetVolume(levels-hfLevels[it])->GetLogicalVolume();
00698 if (lv == hfLV[it]) return true;
00699 }
00700 }
00701 return false;
00702 }
00703
00704 bool HCalSD::isItHF (G4String name) {
00705 std::vector<G4String>::const_iterator it = hfNames.begin();
00706 for (; it != hfNames.end(); ++it) if (name == *it) return true;
00707 return false;
00708 }
00709
00710 bool HCalSD::isItFibre (G4LogicalVolume* lv) {
00711 std::vector<G4LogicalVolume*>::const_iterator ite = fibreLV.begin();
00712 for (; ite != fibreLV.end(); ++ite) if (lv == *ite) return true;
00713 return false;
00714 }
00715
00716 bool HCalSD::isItFibre (G4String name) {
00717 std::vector<G4String>::const_iterator it = fibreNames.begin();
00718 for (; it != fibreNames.end(); ++it) if (name == *it) return true;
00719 return false;
00720 }
00721
00722 bool HCalSD::isItPMT (G4LogicalVolume* lv) {
00723 std::vector<G4LogicalVolume*>::const_iterator ite = pmtLV.begin();
00724 for (; ite != pmtLV.end(); ++ite) if (lv == *ite) return true;
00725 return false;
00726 }
00727
00728 bool HCalSD::isItStraightBundle (G4LogicalVolume* lv) {
00729 std::vector<G4LogicalVolume*>::const_iterator ite = fibre1LV.begin();
00730 for (; ite != fibre1LV.end(); ++ite) if (lv == *ite) return true;
00731 return false;
00732 }
00733
00734 bool HCalSD::isItConicalBundle (G4LogicalVolume* lv) {
00735 std::vector<G4LogicalVolume*>::const_iterator ite = fibre2LV.begin();
00736 for (; ite != fibre2LV.end(); ++ite) if (lv == *ite) return true;
00737 return false;
00738 }
00739
00740 bool HCalSD::isItScintillator (G4Material* mat) {
00741 std::vector<G4Material*>::const_iterator ite = materials.begin();
00742 for (; ite != materials.end(); ++ite) if (mat == *ite) return true;
00743 return false;
00744 }
00745
00746 bool HCalSD::isItinFidVolume (G4ThreeVector& hitPoint) {
00747 bool flag = true;
00748 if (applyFidCut) {
00749 int npmt = HFFibreFiducial:: PMTNumber(hitPoint);
00750 #ifdef DebugLog
00751 edm::LogInfo("HcalSim") << "HCalSD::isItinFidVolume:#PMT= " << npmt
00752 << " for hit point " << hitPoint;
00753 #endif
00754 if (npmt <= 0) flag = false;
00755 }
00756 #ifdef DebugLog
00757 edm::LogInfo("HcalSim") << "HCalSD::isItinFidVolume: point " << hitPoint
00758 << " return flag " << flag;
00759 #endif
00760 return flag;
00761 }
00762
00763 void HCalSD::getFromLibrary (G4Step* aStep, double weight) {
00764 preStepPoint = aStep->GetPreStepPoint();
00765 theTrack = aStep->GetTrack();
00766 int det = 5;
00767 bool ok;
00768
00769 std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, ok, weight, false);
00770
00771 double etrack = preStepPoint->GetKineticEnergy();
00772 int primaryID = setTrackID(aStep);
00773
00774
00775 posGlobal = preStepPoint->GetPosition();
00776 resetForNewPrimary(posGlobal, etrack);
00777
00778 G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
00779 if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
00780 edepositEM = 1.*GeV;
00781 edepositHAD = 0.;
00782 } else {
00783 edepositEM = 0.;
00784 edepositHAD = 1.*GeV;
00785 }
00786 #ifdef DebugLog
00787 edm::LogInfo("HcalSim") << "HCalSD::getFromLibrary " <<hits.size()
00788 << " hits for " << GetName() << " of " << primaryID
00789 << " with " << theTrack->GetDefinition()->GetParticleName()
00790 << " of " << preStepPoint->GetKineticEnergy()/GeV << " GeV";
00791 #endif
00792 for (unsigned int i=0; i<hits.size(); ++i) {
00793 G4ThreeVector hitPoint = hits[i].position;
00794 if (isItinFidVolume (hitPoint)) {
00795 int depth = hits[i].depth;
00796 double time = hits[i].time;
00797 unsigned int unitID = setDetUnitId(det, hitPoint, depth);
00798 currentID.setID(unitID, time, primaryID, 0);
00799 #ifdef plotDebug
00800 plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
00801 bool emType = false;
00802 if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
00803 emType = true;
00804 plotHF(hitPoint,emType);
00805 #endif
00806
00807
00808 if (currentID == previousID) {
00809 updateHit(currentHit);
00810 } else {
00811 if (!checkHit()) currentHit = createNewHit();
00812 }
00813 }
00814 }
00815
00816
00817 if (ok) {
00818 theTrack->SetTrackStatus(fStopAndKill);
00819 G4TrackVector tv = *(aStep->GetSecondary());
00820 for (unsigned int kk=0; kk<tv.size(); ++kk)
00821 if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume())
00822 tv[kk]->SetTrackStatus(fStopAndKill);
00823 }
00824 }
00825
00826 void HCalSD::hitForFibre (G4Step* aStep, double weight) {
00827
00828 preStepPoint = aStep->GetPreStepPoint();
00829 theTrack = aStep->GetTrack();
00830 int primaryID = setTrackID(aStep);
00831
00832 int det = 5;
00833 std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight);
00834
00835 G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
00836 if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
00837 edepositEM = 1.*GeV;
00838 edepositHAD = 0.;
00839 } else {
00840 edepositEM = 0.;
00841 edepositHAD = 1.*GeV;
00842 }
00843
00844 #ifdef DebugLog
00845 edm::LogInfo("HcalSim") << "HCalSD::hitForFibre " << hits.size()
00846 << " hits for " << GetName() << " of " << primaryID
00847 << " with " << theTrack->GetDefinition()->GetParticleName()
00848 << " of " << preStepPoint->GetKineticEnergy()/GeV
00849 << " GeV in detector type " << det;
00850 #endif
00851 if (hits.size() > 0) {
00852 for (unsigned int i=0; i<hits.size(); ++i) {
00853 G4ThreeVector hitPoint = hits[i].position;
00854 if (isItinFidVolume (hitPoint)) {
00855 int depth = hits[i].depth;
00856 double time = hits[i].time;
00857 unsigned int unitID = setDetUnitId(det, hitPoint, depth);
00858 currentID.setID(unitID, time, primaryID, 0);
00859 #ifdef plotDebug
00860 plotProfile(aStep, hitPoint, edepositEM, time, depth);
00861 bool emType = false;
00862 if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
00863 emType = true;
00864 plotHF(hitPoint,emType);
00865 #endif
00866
00867 if (currentID == previousID) {
00868 updateHit(currentHit);
00869 } else {
00870 posGlobal = preStepPoint->GetPosition();
00871 if (!checkHit()) currentHit = createNewHit();
00872 }
00873 }
00874 }
00875 }
00876 }
00877
00878 void HCalSD::getFromParam (G4Step* aStep, double weight) {
00879 std::vector<HFShowerParam::Hit> hits = showerParam->getHits(aStep, weight);
00880 int nHit = static_cast<int>(hits.size());
00881
00882 if (nHit > 0) {
00883 preStepPoint = aStep->GetPreStepPoint();
00884 int primaryID = setTrackID(aStep);
00885
00886 int det = 5;
00887 #ifdef DebugLog
00888 edm::LogInfo("HcalSim") << "HCalSD::getFromParam " << nHit << " hits for "
00889 << GetName() << " of " << primaryID << " with "
00890 << aStep->GetTrack()->GetDefinition()->GetParticleName()
00891 << " of " << preStepPoint->GetKineticEnergy()/GeV
00892 << " GeV in detector type " << det;
00893 #endif
00894 for (int i = 0; i<nHit; ++i) {
00895 G4ThreeVector hitPoint = hits[i].position;
00896 int depth = hits[i].depth;
00897 double time = hits[i].time;
00898 unsigned int unitID = setDetUnitId(det, hitPoint, depth);
00899 currentID.setID(unitID, time, primaryID, 0);
00900 edepositEM = hits[i].edep*GeV;
00901 edepositHAD = 0.;
00902 #ifdef plotDebug
00903 plotProfile(aStep, hitPoint, edepositEM, time, depth);
00904 #endif
00905
00906
00907 if (currentID == previousID) {
00908 updateHit(currentHit);
00909 } else {
00910 posGlobal = preStepPoint->GetPosition();
00911 if (!checkHit()) currentHit = createNewHit();
00912 }
00913 }
00914 }
00915 }
00916
00917 void HCalSD::getHitPMT (G4Step * aStep) {
00918
00919 preStepPoint = aStep->GetPreStepPoint();
00920 theTrack = aStep->GetTrack();
00921 double edep = showerPMT->getHits(aStep);
00922
00923 if (edep >= 0) {
00924 double etrack = preStepPoint->GetKineticEnergy();
00925 int primaryID = 0;
00926 if (etrack >= energyCut) {
00927 primaryID = theTrack->GetTrackID();
00928 } else {
00929 primaryID = theTrack->GetParentID();
00930 if (primaryID == 0) primaryID = theTrack->GetTrackID();
00931 }
00932
00933 posGlobal = preStepPoint->GetPosition();
00934 resetForNewPrimary(posGlobal, etrack);
00935
00936
00937 int det = static_cast<int>(HcalForward);
00938 G4ThreeVector hitPoint = preStepPoint->GetPosition();
00939 double rr = (hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
00940 double phi = (rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x()));
00941 double etaR = showerPMT->getRadius();
00942 int depth = 1;
00943 if (etaR < 0) {
00944 depth = 2;
00945 etaR =-etaR;
00946 }
00947 if (hitPoint.z() < 0) etaR =-etaR;
00948 #ifdef DebugLog
00949 edm::LogInfo("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR "
00950 << etaR << " phi " << phi/deg << " depth " <<depth;
00951 #endif
00952 double time = (aStep->GetPostStepPoint()->GetGlobalTime());
00953 uint32_t unitID = 0;
00954 if (numberingFromDDD) {
00955 HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det,etaR,phi,
00956 depth,1);
00957 if (numberingScheme) unitID = numberingScheme->getUnitID(tmp);
00958 }
00959 currentID.setID(unitID, time, primaryID, 1);
00960
00961 edepositHAD = aStep->GetTotalEnergyDeposit();
00962 edepositEM =-edepositHAD + (edep*GeV);
00963 #ifdef plotDebug
00964 plotProfile(aStep, hitPoint, edep*GeV, time, depth);
00965 #endif
00966 #ifdef DebugLog
00967 double beta = preStepPoint->GetBeta();
00968 LogDebug("HcalSim") << "HCalSD::getHitPMT 1 hit for " << GetName()
00969 << " of " << primaryID << " with "
00970 << theTrack->GetDefinition()->GetParticleName()
00971 << " of " << preStepPoint->GetKineticEnergy()/GeV
00972 << " GeV with velocity " << beta << " UnitID "
00973 << std::hex << unitID << std::dec;
00974 #endif
00975
00976 if (currentID == previousID) {
00977 updateHit(currentHit);
00978 } else {
00979 if (!checkHit()) currentHit = createNewHit();
00980 }
00981 }
00982 }
00983
00984 void HCalSD::getHitFibreBundle (G4Step* aStep, bool type) {
00985 preStepPoint = aStep->GetPreStepPoint();
00986 theTrack = aStep->GetTrack();
00987 double edep = showerBundle->getHits(aStep, type);
00988
00989 if (edep >= 0) {
00990 double etrack = preStepPoint->GetKineticEnergy();
00991 int primaryID = 0;
00992 if (etrack >= energyCut) {
00993 primaryID = theTrack->GetTrackID();
00994 } else {
00995 primaryID = theTrack->GetParentID();
00996 if (primaryID == 0) primaryID = theTrack->GetTrackID();
00997 }
00998
00999 posGlobal = preStepPoint->GetPosition();
01000 resetForNewPrimary(posGlobal, etrack);
01001
01002
01003 int det = static_cast<int>(HcalForward);
01004 G4ThreeVector hitPoint = preStepPoint->GetPosition();
01005 double rr = (hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
01006 double phi = (rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x()));
01007 double etaR = showerBundle->getRadius();
01008 int depth = 1;
01009 if (etaR < 0) {
01010 depth = 2;
01011 etaR =-etaR;
01012 }
01013 if (hitPoint.z() < 0) etaR =-etaR;
01014 #ifdef DebugLog
01015 LogDebug("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR "
01016 << etaR << " phi " << phi/deg << " depth " <<depth;
01017 #endif
01018 double time = (aStep->GetPostStepPoint()->GetGlobalTime());
01019 uint32_t unitID = 0;
01020 if (numberingFromDDD) {
01021 HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det,etaR,phi,depth,1);
01022 if (numberingScheme) unitID = numberingScheme->getUnitID(tmp);
01023 }
01024 if (type) currentID.setID(unitID, time, primaryID, 3);
01025 else currentID.setID(unitID, time, primaryID, 2);
01026
01027 edepositHAD = aStep->GetTotalEnergyDeposit();
01028 edepositEM =-edepositHAD + (edep*GeV);
01029 #ifdef plotDebug
01030 plotProfile(aStep, hitPoint, edep*GeV, time, depth);
01031 #endif
01032 #ifdef DebugLog
01033 double beta = preStepPoint->GetBeta();
01034 LogDebug("HcalSim") << "HCalSD::getHitFibreBundle 1 hit for " << GetName()
01035 << " of " << primaryID << " with "
01036 << theTrack->GetDefinition()->GetParticleName()
01037 << " of " << preStepPoint->GetKineticEnergy()/GeV
01038 << " GeV with velocity " << beta << " UnitID "
01039 << std::hex << unitID << std::dec;
01040 #endif
01041
01042 if (currentID == previousID) updateHit(currentHit);
01043 else if (!checkHit()) currentHit = createNewHit();
01044 }
01045 }
01046
01047 int HCalSD::setTrackID (G4Step* aStep) {
01048 theTrack = aStep->GetTrack();
01049
01050 double etrack = preStepPoint->GetKineticEnergy();
01051 TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
01052 int primaryID = trkInfo->getIDonCaloSurface();
01053 if (primaryID == 0) {
01054 #ifdef DebugLog
01055 edm::LogInfo("HcalSim") << "HCalSD: Problem with primaryID **** set by "
01056 << "force to TkID **** " <<theTrack->GetTrackID();
01057 #endif
01058 primaryID = theTrack->GetTrackID();
01059 }
01060
01061 if (primaryID != previousID.trackID())
01062 resetForNewPrimary(preStepPoint->GetPosition(), etrack);
01063
01064 return primaryID;
01065 }
01066
01067 void HCalSD::readWeightFromFile(std::string fName) {
01068
01069 std::ifstream infile;
01070 int entry=0;
01071 infile.open(fName.c_str(), std::ios::in);
01072 if (infile) {
01073 int det, zside, etaR, phi, lay;
01074 double wt;
01075 while (infile >> det >> zside >> etaR >> phi >> lay >> wt) {
01076 uint32_t id = HcalTestNumbering::packHcalIndex(det,zside,1,etaR,phi,lay);
01077 layerWeights.insert(std::pair<uint32_t,double>(id,wt));
01078 ++entry;
01079 #ifdef DebugLog
01080 edm::LogInfo("HcalSim") << "HCalSD::readWeightFromFile:Entry " << entry
01081 << " ID " << std::hex << id << std::dec << " ("
01082 << det << "/" << zside << "/1/" << etaR << "/"
01083 << phi << "/" << lay << ") Weight " << wt;
01084 #endif
01085 }
01086 infile.close();
01087 }
01088 edm::LogInfo("HcalSim") << "HCalSD::readWeightFromFile: reads " << entry
01089 << " weights from " << fName;
01090 if (entry <= 0) useLayerWt = false;
01091 }
01092
01093 double HCalSD::layerWeight(int det, G4ThreeVector pos, int depth, int lay) {
01094
01095 double wt = 1.;
01096 if (numberingFromDDD) {
01097
01098 HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, pos,
01099 depth, lay);
01100 uint32_t id = HcalTestNumbering::packHcalIndex(tmp.subdet, tmp.zside, 1,
01101 tmp.etaR, tmp.phis,tmp.lay);
01102 std::map<uint32_t,double>::const_iterator ite = layerWeights.find(id);
01103 if (ite != layerWeights.end()) wt = ite->second;
01104 #ifdef DebugLog
01105 edm::LogInfo("HcalSim") << "HCalSD::layerWeight: ID " << std::hex << id
01106 << std::dec << " (" << tmp.subdet << "/"
01107 << tmp.zside << "/1/" << tmp.etaR << "/"
01108 << tmp.phis << "/" << tmp.lay << ") Weight " <<wt;
01109 #endif
01110 }
01111 return wt;
01112 }
01113
01114 void HCalSD::plotProfile(G4Step* aStep, G4ThreeVector global, double edep,
01115 double time, int id) {
01116
01117 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
01118 static G4String modName[8] = {"HEModule", "HVQF" , "HBModule", "MBAT",
01119 "MBBT" , "MBBTC", "MBBT_R1P", "MBBT_R1M"};
01120 G4ThreeVector local;
01121 bool found=false;
01122 double depth=-2000;
01123 int idx = 4;
01124 for (int n=0; n<touch->GetHistoryDepth(); ++n) {
01125 G4String name = touch->GetVolume(n)->GetName();
01126 #ifdef DebugLog
01127 LogDebug("HcalSim") << "plotProfile Depth " << n << " Name " << name;
01128 #endif
01129 for (unsigned int ii=0; ii<8; ++ii) {
01130 if (name == modName[ii]) {
01131 found = true;
01132 int dn = touch->GetHistoryDepth() - n;
01133 local = touch->GetHistory()->GetTransform(dn).TransformPoint(global);
01134 if (ii == 0) {depth = local.z() - 4006.5; idx = 1;}
01135 else if (ii == 1) {depth = local.z() + 825.0; idx = 3;}
01136 else if (ii == 2) {depth = local.x() - 1775.; idx = 0;}
01137 else {depth = local.y() + 15.; idx = 2;}
01138 break;
01139 }
01140 }
01141 if (found) break;
01142 }
01143 if (!found) depth = std::abs(global.z()) - 11500;
01144 #ifdef DebugLog
01145 LogDebug("HcalSim") << "plotProfile Found " << found << " Global " << global
01146 << " Local " << local << " depth " << depth << " ID "
01147 << id << " EDEP " << edep << " Time " << time;
01148 #endif
01149 if (hit_[idx] != 0) hit_[idx]->Fill(edep);
01150 if (time_[idx] != 0) time_[idx]->Fill(time,edep);
01151 if (dist_[idx] != 0) dist_[idx]->Fill(depth,edep);
01152 int jd = 2*idx + id - 7;
01153 if (jd >= 0 && jd < 4) {
01154 jd += 5;
01155 if (hit_[jd] != 0) hit_[jd]->Fill(edep);
01156 if (time_[jd] != 0) time_[jd]->Fill(time,edep);
01157 if (dist_[jd] != 0) dist_[jd]->Fill(depth,edep);
01158 }
01159 }
01160
01161 void HCalSD::plotHF(G4ThreeVector& hitPoint, bool emType) {
01162 double zv = std::abs(hitPoint.z()) - gpar[4];
01163 if (emType) {
01164 if (hzvem != 0) hzvem->Fill(zv);
01165 } else {
01166 if (hzvhad != 0) hzvhad->Fill(zv);
01167 }
01168 }