CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/SimG4CMS/Calo/src/HCalSD.cc

Go to the documentation of this file.
00001 
00002 // File: HCalSD.cc
00003 // Description: Sensitive Detector class for calorimeters
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 //#define DebugLog
00033 //#define plotDebug
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   //static SimpleConfigurable<bool>   on1(false, "HCalSD:UseBirkLaw");
00046   //static SimpleConfigurable<double> bk1(0.013, "HCalSD:BirkC1");
00047   //static SimpleConfigurable<double> bk2(0.0568,"HCalSD:BirkC2");
00048   //static SimpleConfigurable<double> bk3(1.75,  "HCalSD:BirkC3");
00049   // Values from NIM 80 (1970) 239-244: as implemented in Geant3
00050   //static SimpleConfigurable<bool> on2(true,"HCalSD:UseShowerLibrary");
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     // HF volume names
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     // HF Fibre volume names
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     // HF PMT volume names
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     // HF Fibre bundles
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     // Geometry parameters for HF
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       //Special Geometry parameters
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   //Material list for HB/HE/HO sensitive detectors
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   //Layer0 Weight
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);//NB:diff. layer count
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())) {         // Only secondary particles
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   if(wt2 != 1.0) { 
00560     std::cout << "HCalSD: Detector " << det+3 << " Depth " << depth
00561               << " weight= " << weight << " wt1= " << wt1 
00562               << " wt2= " << wt2 << std::endl;
00563     std::cout << theTrack->GetDefinition()->GetParticleName()
00564               << " " << theTrack->GetKineticEnergy()
00565               << " Id=" << theTrack->GetTrackID()
00566               << " IdP=" << theTrack->GetParentID();
00567     const G4VProcess* pr = theTrack->GetCreatorProcess();
00568     if(pr) std::cout << " from  " << pr->GetProcessName();
00569     std::cout << std::endl;
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     //get the ID's as eta, phi, depth, ... indices
00640     HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, pos, depth, lay);
00641     //get the ID
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   // Reset entry point for new primary
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       // check if it is in the same unit and timeslice as the previous one
00808       if (currentID == previousID) {
00809         updateHit(currentHit);
00810       } else {
00811         if (!checkHit()) currentHit = createNewHit();
00812       }
00813     }
00814   }
00815 
00816   //Now kill the current track
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) { // if not ParamShower
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         // check if it is in the same unit and timeslice as the previous one
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       // check if it is in the same unit and timeslice as the previous one
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     // Reset entry point for new primary
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     // check if it is in the same unit and timeslice as the previous one
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     // Reset entry point for new primary
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     // check if it is in the same unit and timeslice as the previous one
01042     if (currentID == previousID) updateHit(currentHit);
01043     else if (!checkHit()) currentHit = createNewHit();
01044   } // non-zero energy deposit
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     //get the ID's as eta, phi, depth, ... indices
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 }