CMSSW Reference Manual

00001 
00002 // File: HFShowerLibrary.cc
00003 // Description: Shower library for Very forward hadron calorimeter
00005 
00006 #include "SimG4CMS/Calo/interface/HFShowerLibrary.h"
00007 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
00008 #include "DetectorDescription/Core/interface/DDFilter.h"
00009 #include "DetectorDescription/Core/interface/DDFilteredView.h"
00010 #include "DetectorDescription/Core/interface/DDValue.h"
00011 
00012 #include "FWCore/Utilities/interface/Exception.h"
00013 
00014 #include "G4VPhysicalVolume.hh"
00015 #include "G4Step.hh"
00016 #include "G4Track.hh"
00017 #include "Randomize.hh"
00018 #include "CLHEP/Units/SystemOfUnits.h"
00019 
00020 HFShowerLibrary::HFShowerLibrary(std::string & name, const DDCompactView & cpv,
00021                                  edm::ParameterSet const & p) : fibre(0),hf(0),
00022                                                                 emBranch(0),
00023                                                                 hadBranch(0),
00024                                                                 nHit(0), 
00025                                                                 npe(0) {
00026   
00027 
00028   edm::ParameterSet m_HF  = p.getParameter<edm::ParameterSet>("HFShower");
00029   probMax                 = m_HF.getParameter<double>("ProbMax");
00030 
00031   edm::ParameterSet m_HS= p.getParameter<edm::ParameterSet>("HFShowerLibrary");
00032   edm::FileInPath fp       = m_HS.getParameter<edm::FileInPath>("FileName");
00033   std::string pTreeName    = fp.fullPath();
00034   backProb                 = m_HS.getParameter<double>("BackProbability");  
00035   std::string emName       = m_HS.getParameter<std::string>("TreeEMID");
00036   std::string hadName      = m_HS.getParameter<std::string>("TreeHadID");
00037   std::string branchEvInfo = m_HS.getUntrackedParameter<std::string>("BranchEvt","HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
00038   std::string branchPre    = m_HS.getUntrackedParameter<std::string>("BranchPre","HFShowerPhotons_hfshowerlib_");
00039   std::string branchPost   = m_HS.getUntrackedParameter<std::string>("BranchPost","_R.obj");
00040   verbose                  = m_HS.getUntrackedParameter<bool>("Verbosity",false);
00041 
00042   if (pTreeName.find(".") == 0) pTreeName.erase(0,2);
00043   const char* nTree = pTreeName.c_str();
00044   hf                = TFile::Open(nTree);
00045 
00046   if (!hf->IsOpen()) { 
00047     edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree 
00048                               << " failed";
00049     throw cms::Exception("Unknown", "HFShowerLibrary") 
00050       << "Opening of " << pTreeName << " fails\n";
00051   } else {
00052     edm::LogInfo("HFShower") << "HFShowerLibrary: opening " << nTree 
00053                              << " successfully"; 
00054   }
00055 
00056   TTree* event = (TTree *) hf ->Get("Events");
00057   if (event) {
00058     std::string info = branchEvInfo + branchPost;
00059     TBranch *evtInfo = event->GetBranch(info.c_str());
00060     if (evtInfo) {
00061       loadEventInfo(evtInfo);
00062     } else {
00063       edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
00064                                 << " Branch does not exist in Event";
00065       throw cms::Exception("Unknown", "HFShowerLibrary")
00066         << "Event information absent\n";
00067     }
00068   } else {
00069     edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
00070                               << "exist";
00071     throw cms::Exception("Unknown", "HFShowerLibrary")
00072       << "Events tree absent\n";
00073   }
00074   
00075   edm::LogInfo("HFShower") << "HFShowerLibrary: Library " << libVers 
00076                            << " ListVersion "   << listVersion 
00077                            << " Events Total " << totEvents << " and "
00078                            << evtPerBin << " per bin";
00079   edm::LogInfo("HFShower") << "HFShowerLibrary: Energies (GeV) with " 
00080                            << nMomBin   << " bins";
00081   for (int i=0; i<nMomBin; i++)
00082     edm::LogInfo("HFShower") << "HFShowerLibrary: pmom[" << i << "] = "
00083                              << pmom[i]/GeV << " GeV";
00084 
00085   std::string nameBr = branchPre + emName + branchPost;
00086   emBranch         = event->GetBranch(nameBr.c_str());
00087   if (verbose) emBranch->Print();
00088   nameBr           = branchPre + hadName + branchPost;
00089   hadBranch        = event->GetBranch(nameBr.c_str());
00090   if (verbose) hadBranch->Print();
00091   edm::LogInfo("HFShower") << "HFShowerLibrary:Branch " << emName 
00092                            << " has " << emBranch->GetEntries() 
00093                            << " entries and Branch " << hadName 
00094                            << " has " << hadBranch->GetEntries() 
00095                            << " entries";
00096   edm::LogInfo("HFShower") << "HFShowerLibrary::No packing information -"
00097                            << " Assume x, y, z are not in packed form";
00098 
00099   edm::LogInfo("HFShower") << "HFShowerLibrary: Maximum probability cut off " 
00100                            << probMax << "  Back propagation of light prob. "
00101                            << backProb ;
00102   
00103   G4String attribute = "ReadOutName";
00104   G4String value     = name;
00105   DDSpecificsFilter filter;
00106   DDValue           ddv(attribute,value,0);
00107   filter.setCriteria(ddv,DDSpecificsFilter::equals);
00108   DDFilteredView fv(cpv);
00109   fv.addFilter(filter);
00110   bool dodet = fv.firstChild();
00111   if (dodet) {
00112     DDsvalues_type sv(fv.mergedSpecifics());
00113 
00114     //Radius (minimum and maximum)
00115     int nR     = -1;
00116     std::vector<double> rTable = getDDDArray("rTable",sv,nR);
00117     rMin = rTable[0];
00118     rMax = rTable[nR-1];
00119     edm::LogInfo("HFShower") << "HFShowerLibrary: rMIN " << rMin/cm 
00120                              << " cm and rMax " << rMax/cm;
00121 
00122     //Delta phi
00123     int nEta   = -1;
00124     std::vector<double> etaTable = getDDDArray("etaTable",sv,nEta);
00125     int nPhi   = nEta + nR - 2;
00126     std::vector<double> phibin   = getDDDArray("phibin",sv,nPhi);
00127     dphi       = phibin[nEta-1];
00128     edm::LogInfo("HFShower") << "HFShowerLibrary: (Half) Phi Width of wedge " 
00129                              << dphi/deg;
00130 
00131     //Special Geometry parameters
00132     int ngpar = 7;
00133     gpar      = getDDDArray("gparHF",sv,ngpar);
00134     edm::LogInfo("HFShower") << "HFShowerLibrary: " << ngpar << " gpar (cm)";
00135     for (int ig=0; ig<ngpar; ig++)
00136       edm::LogInfo("HFShower") << "HFShowerLibrary: gpar[" << ig << "] = "
00137                                << gpar[ig]/cm << " cm";
00138   } else {
00139     edm::LogError("HFShower") << "HFShowerLibrary: cannot get filtered "
00140                               << " view for " << attribute << " matching "
00141                               << name;
00142     throw cms::Exception("Unknown", "HFShowerLibrary")
00143       << "cannot match " << attribute << " to " << name <<"\n";
00144   }
00145   
00146   fibre = new HFFibre(name, cpv, p);
00147   emPDG = epPDG = gammaPDG = 0;
00148   pi0PDG = etaPDG = nuePDG = numuPDG = nutauPDG= 0;
00149   anuePDG= anumuPDG = anutauPDG = geantinoPDG = 0;
00150 }
00151 
00152 HFShowerLibrary::~HFShowerLibrary() {
00153   if (hf)     hf->Close();
00154   if (fibre)  delete   fibre;  fibre  = 0;
00155 }
00156 
00157 void HFShowerLibrary::initRun(G4ParticleTable * theParticleTable) {
00158 
00159   G4String parName;
00160   emPDG = theParticleTable->FindParticle(parName="e-")->GetPDGEncoding();
00161   epPDG = theParticleTable->FindParticle(parName="e+")->GetPDGEncoding();
00162   gammaPDG = theParticleTable->FindParticle(parName="gamma")->GetPDGEncoding();
00163   pi0PDG = theParticleTable->FindParticle(parName="pi0")->GetPDGEncoding();
00164   etaPDG = theParticleTable->FindParticle(parName="eta")->GetPDGEncoding();
00165   nuePDG = theParticleTable->FindParticle(parName="nu_e")->GetPDGEncoding();
00166   numuPDG = theParticleTable->FindParticle(parName="nu_mu")->GetPDGEncoding();
00167   nutauPDG= theParticleTable->FindParticle(parName="nu_tau")->GetPDGEncoding();
00168   anuePDG = theParticleTable->FindParticle(parName="anti_nu_e")->GetPDGEncoding();
00169   anumuPDG= theParticleTable->FindParticle(parName="anti_nu_mu")->GetPDGEncoding();
00170   anutauPDG= theParticleTable->FindParticle(parName="anti_nu_tau")->GetPDGEncoding();
00171   geantinoPDG= theParticleTable->FindParticle(parName="geantino")->GetPDGEncoding();
00172   LogDebug("HFShower") << "HFShowerLibrary: Particle codes for e- = " << emPDG
00173                        << ", e+ = " << epPDG << ", gamma = " << gammaPDG 
00174                        << ", pi0 = " << pi0PDG << ", eta = " << etaPDG
00175                        << ", geantino = " << geantinoPDG << "\n        nu_e = "
00176                        << nuePDG << ", nu_mu = " << numuPDG << ", nu_tau = "
00177                        << nutauPDG << ", anti_nu_e = " << anuePDG
00178                        << ", anti_nu_mu = " << anumuPDG << ", anti_nu_tau = "
00179                        << anutauPDG;
00180 }
00181 
00182 int HFShowerLibrary::getHits(G4Step * aStep) {
00183 
00184 
00185   G4StepPoint * preStepPoint  = aStep->GetPreStepPoint(); 
00186   G4StepPoint * postStepPoint = aStep->GetPostStepPoint(); 
00187   G4Track *     track    = aStep->GetTrack();
00188   // Get Z-direction 
00189   const G4DynamicParticle *aParticle = track->GetDynamicParticle();
00190   G4ThreeVector momDir = aParticle->GetMomentumDirection();
00191   //  double mom = aParticle->GetTotalMomentum();
00192 
00193   G4ThreeVector hitPoint = preStepPoint->GetPosition();   
00194   G4String      partType = track->GetDefinition()->GetParticleName();
00195   int           parCode  = track->GetDefinition()->GetPDGEncoding();
00196 
00197   double tSlice = (postStepPoint->GetGlobalTime())/nanosecond;
00198   double pin    = preStepPoint->GetTotalEnergy();
00199   
00200   double px   = momDir.x(); 
00201   double py   = momDir.y(); 
00202   double pz   = momDir.z(); 
00203 
00204   double xint = hitPoint.x(); 
00205   double yint = hitPoint.y(); 
00206   double zint = hitPoint.z(); 
00207 
00208   // if particle moves from interaction point or "backwards (halo)
00209   int backward = 0;
00210   if(pz * zint < 0.) backward = 1;
00211   
00212   double sphi   = sin(momDir.phi());
00213   double cphi   = cos(momDir.phi());
00214   double ctheta = cos(momDir.theta());
00215   double stheta = sin(momDir.theta());
00216  
00217   /*  // Obsolete piece, only valid if particle comes from IP...
00218   double sphi   = sin(hitPoint.phi());
00219   double cphi   = cos(hitPoint.phi());
00220   double ctheta = cos(hitPoint.theta());
00221   double stheta = sin(hitPoint.theta());
00222   */
00223 
00224   LogDebug("HFShower") << "HFShowerLibrary: getHits " << partType
00225                        << " of energy " << pin/GeV << " GeV"
00226                        << "  dir.orts " << px << ", " << py << ", " << pz
00227                        << "  Pos x,y,z = " << xint << "," << yint << "," 
00228                        << zint << "   sphi,cphi,stheta,ctheta  =" 
00229                        << sphi << "," << cphi << ","   
00230                        << stheta << "," << ctheta ; 
00231     
00232                        
00233   if (parCode == pi0PDG || parCode == etaPDG || parCode == nuePDG ||
00234       parCode == numuPDG || parCode == nutauPDG || parCode == anuePDG ||
00235       parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG) {
00236     return -1;
00237   } else if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
00238     if (pin<pmom[nMomBin-1]) {
00239       interpolate(0, pin);
00240     } else {
00241       extrapolate(0, pin);
00242     }
00243   } else {
00244     if (pin<pmom[nMomBin-1]) {
00245       interpolate(1, pin);
00246     } else {
00247       extrapolate(1, pin);
00248     }
00249   }
00250     
00251   nHit = 0;
00252   if (npe > 0) {
00253     hit.clear(); hit.resize(npe);
00254   }
00255   for (int i = 0; i < npe; i++) {
00256     LogDebug("HFShower") << "HFShowerLibrary: Hit " << i << " " << pe[i];
00257     double zv = std::abs(pe[i].z()); // abs local z  
00258     if (zv <= gpar[1] && pe[i].lambda() > 0 &&
00259         (pe[i].z() >= 0 || pe[i].z() <= -gpar[0])) {
00260       int depth = 1;
00261       if(backward == 0) {            // fully valid only for "front" particles
00262         if (pe[i].z() < 0) depth = 2;// with "front"-simulated shower lib.
00263       }
00264       else {                         // for "backward" particles - almost equal
00265         double r = G4UniformRand();  // share between L and S fibers
00266         if(r > 0.5) depth = 2;
00267       } 
00268       
00269 
00270       // Updated coordinate transformation from local
00271       //  back to global using two Euler angles: phi and theta
00272       double pex = pe[i].x();
00273       double pey = pe[i].y();
00274 
00275       double xx = pex*ctheta*cphi - pey*sphi + zv*stheta*cphi; 
00276       double yy = pex*ctheta*sphi + pey*cphi + zv*stheta*sphi;
00277       double zz = -pex*stheta + zv*ctheta;
00278 
00279       G4ThreeVector pos = hitPoint + G4ThreeVector(xx,yy,zz);
00280 
00281       if(backward == 0) zv = gpar[1] - zv;      // remaining distance to PMT !
00282 
00283       double r  = pos.perp();
00284       double p  = fibre->attLength(pe[i].lambda());
00285       double fi = pos.phi();
00286       if (fi < 0) fi += twopi;
00287       int    isect = int(fi/dphi) + 1;
00288       isect        = (isect + 1) / 2;
00289       double dfi   = ((isect*2-1)*dphi - fi);
00290       if (dfi < 0) dfi = -dfi;
00291       double dfir  = r * sin(dfi);
00292       LogDebug("HFShower") << "HFShowerLibrary: Position shift " << xx 
00293                            << ", " << yy 
00294                            << ", "  << zz << ": " << pos << " R " << r 
00295                            << " Phi " << fi << " Section " << isect 
00296                            << " R*Dfi " << dfir;
00297       zz           = ((pos.z()) >= 0 ? (pos.z()) : -(pos.z()));
00298       double r1    = G4UniformRand();
00299       double r2    = G4UniformRand();
00300       double r3    = -9999.;
00301       if(backward) r3 = G4UniformRand();
00302 
00303       LogDebug("HFShower") << "  rLimits " << rInside(r)
00304                            << " attenuation " << r1 <<":" << exp(-p*zv) 
00305                            << " r2 " << r2 << " r3 " << r3 << " rDfi "  
00306                            << gpar[5] << " zz " 
00307                            << zz << " zLim " << gpar[4] << ":" 
00308                            << gpar[4]+gpar[1];
00309       LogDebug("HFShower") << "  rInside(r) :" << rInside(r) 
00310                            << "  r1 <= exp(-p*zv) :" <<  (r1 <= exp(-p*zv))
00311                            << "  r2 <= probMax :"    <<  (r2 <= probMax)
00312                            << "  r3 <= backProb :"   <<  (r3 <= backProb) 
00313                            << "  dfir > gpar[5] :"   <<  (dfir > gpar[5])
00314                            << "  zz >= gpar[4] :"    <<  (zz >= gpar[4])
00315                            << "  zz <= gpar[4]+gpar[1] :" 
00316                            << (zz <= gpar[4]+gpar[1]);   
00317 
00318       if (rInside(r) && r1 <= exp(-p*zv) && r2 <= probMax && dfir > gpar[5] &&
00319           zz >= gpar[4] && zz <= gpar[4]+gpar[1] && r3 <= backProb ){
00320         hit[nHit].position = pos;
00321         hit[nHit].depth    = depth;
00322         hit[nHit].time     = (tSlice+(pe[i].t())+(fibre->tShift(pos,depth,true)));
00323         LogDebug("HFShower") << "HFShowerLibrary: Final Hit " << nHit 
00324                              <<" position " << (hit[nHit].position) <<" Depth "
00325                              <<(hit[nHit].depth) <<" Time " <<(hit[nHit].time);
00326         nHit++;
00327       }
00328       else  LogDebug("HFShower") << " REJECTED !!!";
00329     }
00330   }
00331 
00332   LogDebug("HFShower") << "HFShowerLibrary: Total Hits " << nHit
00333                        << " out of " << npe << " PE";
00334   if (nHit > npe)
00335     edm::LogWarning("HFShower") << "HFShowerLibrary: Hit buffer " << npe 
00336                                 << " smaller than " << nHit << " Hits";
00337   return nHit;
00338 
00339 }
00340 
00341 G4ThreeVector HFShowerLibrary::getPosHit(int i) {
00342 
00343   G4ThreeVector pos;
00344   if (i < nHit) pos = (hit[i].position);
00345   LogDebug("HFShower") << " HFShowerLibrary: getPosHit (" << i << "/" << nHit 
00346                        << ") " << pos;
00347   return pos;
00348 }
00349 
00350 int HFShowerLibrary::getDepth(int i) {
00351 
00352   int depth = 0;
00353   if (i < nHit) depth = (hit[i].depth);
00354   LogDebug("HFShower") << " HFShowerLibrary: getDepth (" << i << "/" << nHit 
00355                        << ") "  << depth;
00356   return depth;
00357 }
00358 
00359 double HFShowerLibrary::getTSlice(int i) {
00360   
00361   double tim = 0.;
00362   if (i < nHit) tim = (hit[i].time);
00363   LogDebug("HFShower") << " HFShowerLibrary: Time (" << i << "/" << nHit 
00364                        << ") "  << tim;
00365   return tim;
00366 }
00367 
00368 bool HFShowerLibrary::rInside(double r) {
00369 
00370   if (r >= rMin && r <= rMax) return true;
00371   else                        return false;
00372 }
00373 
00374 void HFShowerLibrary::getRecord(int type, int record) {
00375 
00376   int nrc     = record-1;
00377   int nPhoton = 0;
00378   photon.clear();
00379   if (type > 0) {
00380     hadBranch->SetAddress(&photon);
00381     hadBranch->GetEntry(nrc);
00382   } else {
00383     emBranch->SetAddress(&photon);
00384     emBranch->GetEntry(nrc);
00385   }
00386   nPhoton = photon.size();
00387   LogDebug("HFShower") << "HFShowerLibrary::getRecord: Record " << record
00388                        << " of type " << type << " with " << nPhoton 
00389                        << " photons";
00390   for (int j = 0; j < nPhoton; j++) 
00391     LogDebug("HFShower") << "Photon " << j << photon[j];
00392 }
00393 
00394 void HFShowerLibrary::loadEventInfo(TBranch* branch) {
00395 
00396   std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
00397   branch->SetAddress(&eventInfoCollection);
00398   branch->GetEntry(0);
00399   edm::LogInfo("HFShower") << "HFShowerLibrary::loadEventInfo loads "
00400                            << " EventInfo Collection of size "
00401                            << eventInfoCollection.size() << " records";
00402   totEvents   = eventInfoCollection[0].totalEvents();
00403   nMomBin     = eventInfoCollection[0].numberOfBins();
00404   evtPerBin   = eventInfoCollection[0].eventsPerBin();
00405   libVers     = eventInfoCollection[0].showerLibraryVersion();
00406   listVersion = eventInfoCollection[0].physListVersion();
00407   pmom        = eventInfoCollection[0].energyBins();
00408   for (unsigned int i=0; i<pmom.size(); i++) 
00409     pmom[i] *= GeV;
00410 }
00411 
00412 void HFShowerLibrary::interpolate(int type, double pin) {
00413 
00414   LogDebug("HFShower") << "HFShowerLibrary:: Interpolate for Energy " <<pin/GeV
00415                        << " GeV with " << nMomBin << " momentum bins and " 
00416                        << evtPerBin << " entries/bin -- total " << totEvents;
00417   int irc[2];
00418   double w = 0.;
00419   double r = G4UniformRand();
00420 
00421   if (pin<pmom[0]) {
00422     w = pin/pmom[0];
00423     irc[1] = int(evtPerBin*r) + 1;
00424     irc[0] = 0;
00425   } else {
00426     for (int j=0; j<nMomBin-1; j++) {
00427       if (pin >= pmom[j] && pin < pmom[j+1]) {
00428         w = (pin-pmom[j])/(pmom[j+1]-pmom[j]);
00429         if (j == nMomBin-2) { 
00430           irc[1] = int(evtPerBin*0.5*r);
00431         } else {
00432           irc[1] = int(evtPerBin*r);
00433         }
00434         irc[1] += (j+1)*evtPerBin + 1;
00435         r = G4UniformRand();
00436         irc[0] = int(evtPerBin*r) + 1 + j*evtPerBin;
00437         if (irc[0]<0) {
00438           edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
00439                                       << irc[0] << " now set to 0";
00440           irc[0] = 0;
00441         } else if (irc[0] > totEvents) {
00442           edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = "
00443                                       << irc[0] << " now set to "<< totEvents;
00444           irc[0] = totEvents;
00445         }
00446       }
00447     }
00448   }
00449   if (irc[1]<1) {
00450     edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " 
00451                                 << irc[1] << " now set to 1";
00452     irc[1] = 1;
00453   } else if (irc[1] > totEvents) {
00454     edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " 
00455                                 << irc[1] << " now set to "<< totEvents;
00456     irc[1] = totEvents;
00457   }
00458 
00459   LogDebug("HFShower") << "HFShowerLibrary:: Select records " << irc[0] 
00460                        << " and " << irc[1] << " with weights " << 1-w 
00461                        << " and " << w;
00462 
00463   pe.clear(); 
00464   npe       = 0;
00465   int npold = 0;
00466   for (int ir=0; ir < 2; ir++) {
00467     if (irc[ir]>0) {
00468       getRecord (type, irc[ir]);
00469       int nPhoton = photon.size();
00470       npold      += nPhoton;
00471       for (int j=0; j<nPhoton; j++) {
00472         r = G4UniformRand();
00473         if ((ir==0 && r > w) || (ir > 0 && r < w)) {
00474           storePhoton (j);
00475         }
00476       }
00477     }
00478   }
00479 
00480   if (npe > npold || (npold == 0 && irc[0] > 0)) 
00481     edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
00482                                 << " records " << irc[0] << " and " << irc[1]
00483                                 << " gives a buffer of " << npold 
00484                                 << " photons and fills " << npe << " *****";
00485   else
00486     LogDebug("HFShower") << "HFShowerLibrary: Interpolation == records " 
00487                          << irc[0] << " and " << irc[1] << " gives a "
00488                          << "buffer of " << npold << " photons and fills "
00489                          << npe << " PE";
00490   for (int j=0; j<npe; j++) {
00491     LogDebug("HFShower") << "Photon " << j << " " << pe[j];
00492   }
00493 }
00494 
00495 void HFShowerLibrary::extrapolate(int type, double pin) {
00496 
00497   int nrec   = int(pin/pmom[nMomBin-1]);
00498   double w   = (pin - pmom[nMomBin-1]*nrec)/pmom[nMomBin-1];
00499   nrec++;
00500   LogDebug("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin 
00501                        << " GeV with " << nMomBin << " momentum bins and " 
00502                        << evtPerBin << " entries/bin -- total " << totEvents 
00503                        << " using " << nrec << " records";
00504   std::vector<int> irc(nrec);
00505 
00506   for (int ir=0; ir<nrec; ir++) {
00507     double r = G4UniformRand();
00508     irc[ir] = int(evtPerBin*0.5*r) +(nMomBin-1)*evtPerBin + 1;
00509     if (irc[ir]<1) {
00510       edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir 
00511                                   << "] = " << irc[ir] << " now set to 1";
00512       irc[ir] = 1;
00513     } else if (irc[ir] > totEvents) {
00514       edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir 
00515                                   << "] = " << irc[ir] << " now set to "
00516                                   << totEvents;
00517       irc[ir] = totEvents;
00518     } else {
00519       LogDebug("HFShower") << "HFShowerLibrary::Extrapolation use irc[" 
00520                            << ir  << "] = " << irc[ir];
00521     }
00522   }
00523 
00524   pe.clear(); 
00525   npe       = 0;
00526   int npold = 0;
00527   for (int ir=0; ir<nrec; ir++) {
00528     if (irc[ir]>0) {
00529       getRecord (type, irc[ir]);
00530       int nPhoton = photon.size();
00531       npold      += nPhoton;
00532       for (int j=0; j<nPhoton; j++) {
00533         double r = G4UniformRand();
00534         if (ir != nrec-1 || r < w) {
00535           storePhoton (j);
00536         }
00537       }
00538       LogDebug("HFShower") << "Record [" << ir << "] = " << irc[ir] 
00539                            << " npold = " << npold;
00540     }
00541   }
00542   LogDebug("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
00543 
00544   if (npe > npold || npold == 0)
00545     edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == "
00546                                 << nrec << " records " << irc[0] << ", " 
00547                                 << irc[1] << ", ... gives a buffer of " <<npold
00548                                 << " photons and fills " << npe 
00549                                 << " *****";
00550   else
00551     LogDebug("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec
00552                          << " records " << irc[0] << ", " << irc[1] 
00553                          << ", ... gives a buffer of " << npold 
00554                          << " photons and fills " << npe << " PE";
00555   for (int j=0; j<npe; j++) {
00556     LogDebug("HFShower") << "Photon " << j << " " << pe[j];
00557   }
00558 }
00559 
00560 void HFShowerLibrary::storePhoton(int j) {
00561 
00562   pe.push_back(photon[j]);
00563   LogDebug("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe " 
00564                        << npe << " " << pe[npe];
00565   npe++;
00566 }
00567 
00568 std::vector<double> HFShowerLibrary::getDDDArray(const std::string & str, 
00569                                                  const DDsvalues_type & sv, 
00570                                                  int & nmin) {
00571 
00572   LogDebug("HFShower") << "HFShowerLibrary:getDDDArray called for " << str 
00573                        << " with nMin " << nmin;
00574 
00575   DDValue value(str);
00576   if (DDfetch(&sv,value)) {
00577     LogDebug("HFShower") << value;
00578     const std::vector<double> & fvec = value.doubles();
00579     int nval = fvec.size();
00580     if (nmin > 0) {
00581       if (nval < nmin) {
00582         edm::LogError("HFShower") << "HFShowerLibrary : # of " << str 
00583                                   << " bins " << nval << " < " << nmin 
00584                                   << " ==> illegal";
00585         throw cms::Exception("Unknown", "HFShowerLibrary")
00586           << "nval < nmin for array " << str << "\n";
00587       }
00588     } else {
00589       if (nval < 2) {
00590         edm::LogError("HFShower") << "HFShowerLibrary : # of " << str 
00591                                   << " bins " << nval << " < 2 ==> illegal"
00592                                   << " (nmin=" << nmin << ")";
00593         throw cms::Exception("Unknown", "HFShowerLibrary")
00594           << "nval < 2 for array " << str << "\n";
00595       }
00596     }
00597     nmin = nval;
00598 
00599     return fvec;
00600   } else {
00601     edm::LogError("HFShower") << "HFShowerLibrary : cannot get array " << str;
00602     throw cms::Exception("Unknown", "HFShowerLibrary") 
00603       << "cannot get array " << str << "\n";
00604   }
00605 }