CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/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 "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 //#define DebugLog
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   //static SimpleConfigurable<bool>   on1(false, "HCalSD:UseBirkLaw");
00042   //static SimpleConfigurable<double> bk1(0.013, "HCalSD:BirkC1");
00043   //static SimpleConfigurable<double> bk2(0.0568,"HCalSD:BirkC2");
00044   //static SimpleConfigurable<double> bk3(1.75,  "HCalSD:BirkC3");
00045   // Values from NIM 80 (1970) 239-244: as implemented in Geant3
00046   //static SimpleConfigurable<bool> on2(true,"HCalSD:UseShowerLibrary");
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     // HF volume names
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     // HF Fibre volume names
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     // HF PMT volume names
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     // HF Fibre bundles
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   //Material list for HB/HE/HO sensitive detectors
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   //Layer0 Weight
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   //  TimeMe t1( theHitTimer, false);
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())) { // Only secondary particles
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     //get the ID's as eta, phi, depth, ... indices
00537     HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, pos, 
00538                                                                 depth, lay);
00539     //get the ID
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   int    primaryID = 0;
00674   if (etrack >= energyCut) {
00675     primaryID    = theTrack->GetTrackID();
00676   } else {
00677     primaryID    = theTrack->GetParentID();
00678     if (primaryID == 0) primaryID = theTrack->GetTrackID();
00679   }
00680   */
00681 
00682   // Reset entry point for new primary
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     // check if it is in the same unit and timeslice as the previous one
00711     if (currentID == previousID) {
00712       updateHit(currentHit);
00713     } else {
00714       if (!checkHit()) currentHit = createNewHit();
00715     }
00716   }
00717 
00718   //Now kill the current track
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       // check if it is in the same unit and timeslice as the previous one
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       // check if it is in the same unit and timeslice as the previous one
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     // Reset entry point for new primary
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     // check if it is in the same unit and timeslice as the previous one
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     // Reset entry point for new primary
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     // check if it is in the same unit and timeslice as the previous one
00940     if (currentID == previousID) {
00941       updateHit(currentHit);
00942     } else {
00943       if (!checkHit()) currentHit = createNewHit();
00944     }
00945   } // non-zero energy deposit
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     //get the ID's as eta, phi, depth, ... indices
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 }