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