#include <HFShowerParam.h>

Classes
struct	Hit

Public Member Functions
std::vector< Hit >	getHits (G4Step *aStep, double weight)

	HFShowerParam (std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p)

void	initRun (G4ParticleTable *)

virtual	~HFShowerParam ()

Private Member Functions
std::vector< double >	getDDDArray (const std::string &, const DDsvalues_type &)

Private Attributes
double	aperture

bool	applyFidCut

double	attLMeanInv

double	edMin

TH2F *	em_2d_1

TH2F *	em_2d_2

TH1F *	em_lateral_1

TH1F *	em_lateral_2

TH1F *	em_long_1

TH1F *	em_long_1_tuned

TH1F *	em_long_2

TH1F *	em_long_gflash

TH1F *	em_long_sl

G4int	emPDG

G4int	epPDG

HFFibre *	fibre

bool	fillHisto

G4int	gammaPDG

HFGflash *	gflash

std::vector< double >	gpar

TH1F *	hzvem

TH1F *	hzvhad

bool	onlyLong

bool	parametrizeLast

double	pePerGeV

double	ref_index

HFShowerLibrary *	showerLibrary

bool	trackEM

Detailed Description

Definition at line 26 of file HFShowerParam.h.

Constructor & Destructor Documentation

HFShowerParam::HFShowerParam	(	std::string &	name,
		const DDCompactView &	cpv,
		edm::ParameterSet const &	p
	)

Definition at line 30 of file HFShowerParam.cc.

References DDFilteredView::addFilter(), aperture, applyFidCut, HFFibre::attLength(), attLMeanInv, funct::cos(), edMin, em_2d_1, em_2d_2, em_lateral_1, em_lateral_2, em_long_1, em_long_1_tuned, em_long_2, em_long_gflash, em_long_sl, DDSpecificsFilter::equals, edm::hlt::Exception, fibre, fillHisto, alcazmumu_cfi::filter, DDFilteredView::firstChild(), getDDDArray(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), gflash, gpar, hzvem, hzvhad, edm::Service< T >::isAvailable(), LogDebug, TFileDirectory::make(), DDFilteredView::mergedSpecifics(), TFileService::mkdir(), mergeVDriftHistosByStation::name, onlyLong, parametrizeLast, pePerGeV, ref_index, DDSpecificsFilter::setCriteria(), showerLibrary, trackEM, and relativeConstraints::value.

                                                         : showerLibrary(0), 
                                                             fibre(0), gflash(0),
                                                             fillHisto(false) { 
   edm::ParameterSet m_HF  = p.getParameter<edm::ParameterSet>("HFShower");
   pePerGeV                = m_HF.getParameter<double>("PEPerGeV");
   trackEM                 = m_HF.getParameter<bool>("TrackEM");
   bool useShowerLibrary   = m_HF.getParameter<bool>("UseShowerLibrary");
   bool useGflash          = m_HF.getParameter<bool>("UseHFGflash");
   edMin                   = m_HF.getParameter<double>("EminLibrary");
   onlyLong                = m_HF.getParameter<bool>("OnlyLong");
   ref_index               = m_HF.getParameter<double>("RefIndex");
   double lambdaMean       = m_HF.getParameter<double>("LambdaMean");
   aperture                = cos(asin(m_HF.getParameter<double>("Aperture")));
   applyFidCut             = m_HF.getParameter<bool>("ApplyFiducialCut");
   parametrizeLast         = m_HF.getUntrackedParameter<bool>("ParametrizeLast",false);
   edm::LogInfo("HFShower") << "HFShowerParam::Use of shower library is set to "
                            << useShowerLibrary << " Use of Gflash is set to "
                            << useGflash << " P.E. per GeV " << pePerGeV
                            << ", ref. index of fibre " << ref_index 
                            << ", Track EM Flag " << trackEM << ", edMin "
                            << edMin << " GeV, use of Short fibre info in"
                            << " shower library set to " << !(onlyLong)
                            << ", use of parametrization for last part set to "
                            << parametrizeLast << ", Mean lambda " << lambdaMean
                << ", aperture (cutoff) " << aperture
                            << ", Application of Fiducial Cut " << applyFidCut;
 
 #ifdef plotDebug
   edm::Service<TFileService> tfile;
   if (tfile.isAvailable()) {
     fillHisto = true;
     LogDebug("HFShower") << "HFShowerParam::Save histos in directory "
              << "ProfileFromParam";
     TFileDirectory showerDir = tfile->mkdir("ProfileFromParam");
     hzvem           = showerDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part);Number of PE",330,0.0,1650.0);
     hzvhad          = showerDir.make<TH1F>("hzvhad","Longitudinal Profile (Had Part);Number of PE",330,0.0,1650.0);
     em_2d_1         = showerDir.make<TH2F>("em_2d_1","Lateral Profile vs. Shower Depth;cm;Events",800,800.0,1600.0,100,50.0,150.0);
     em_long_1       = showerDir.make<TH1F>("em_long_1","Longitudinal Profile;Radiation Length;Number of Spots",800,800.0,1600.0);
     em_long_1_tuned = showerDir.make<TH1F>("em_long_1_tuned","Longitudinal Profile;Radiation Length;Number of Spots",800,800.0,1600.0);
     em_lateral_1    = showerDir.make<TH1F>("em_lateral_1","Lateral Profile;cm;Events",100,50.0,150.0);
     em_2d_2         = showerDir.make<TH2F>("em_2d_2","Lateral Profile vs. Shower Depth;cm;Events",800,800.0,1600.0,100,50.0,150.0);
     em_long_2       = showerDir.make<TH1F>("em_long_2","Longitudinal Profile;Radiation Length;Number of Spots",800,800.0,1600.0);
     em_lateral_2    = showerDir.make<TH1F>("em_lateral_2","Lateral Profile;cm;Events",100,50.0,150.0);
     em_long_gflash  = showerDir.make<TH1F>("em_long_gflash","Longitudinal Profile From GFlash;cm;Number of Spots",800,800.0,1600.0);
     em_long_sl      = showerDir.make<TH1F>("em_long_sl","Longitudinal Profile From Shower Library;cm;Number of Spots",800,800.0,1600.0);
   } else {
     fillHisto = false;
     edm::LogInfo("HFShower") << "HFShowerParam::No file is available for "
                  << "saving histos so the flag is set to false";
   }
 #endif
   
   G4String attribute = "ReadOutName";
   G4String value     = name;
   DDSpecificsFilter filter;
   DDValue           ddv(attribute,value,0);
   filter.setCriteria(ddv,DDSpecificsFilter::equals);
   DDFilteredView fv(cpv);
   fv.addFilter(filter);
   bool dodet = fv.firstChild();
   if (dodet) {
     DDsvalues_type sv(fv.mergedSpecifics());
     //Special Geometry parameters
     gpar      = getDDDArray("gparHF",sv);
     edm::LogInfo("HFShower") << "HFShowerParam: " <<gpar.size() <<" gpar (cm)";
     for (unsigned int ig=0; ig<gpar.size(); ig++)
       edm::LogInfo("HFShower") << "HFShowerParam: gpar[" << ig << "] = "
                                << gpar[ig]/cm << " cm";
   } else {
     edm::LogError("HFShower") << "HFShowerParam: cannot get filtered "
                               << " view for " << attribute << " matching " << name;
     throw cms::Exception("Unknown", "HFShowerParam") << "cannot match " << attribute
                                                      << " to " << name <<"\n";
   }
   
   if (useShowerLibrary) showerLibrary = new HFShowerLibrary(name, cpv, p);
   if (useGflash)        gflash        = new HFGflash(p);
   fibre = new HFFibre(name, cpv, p);
   attLMeanInv = fibre->attLength(lambdaMean);
   edm::LogInfo("HFShower") << "att. length used for (lambda=" << lambdaMean
                            << ") = " << 1/(attLMeanInv*cm) << " cm";
 }

HFShowerParam::~HFShowerParam ( )

virtual

Definition at line 114 of file HFShowerParam.cc.

References fibre, gflash, and showerLibrary.

                               {
   if (fibre)         delete fibre;
   if (gflash)        delete gflash;
   if (showerLibrary) delete showerLibrary;
 }

Member Function Documentation

std::vector< double > HFShowerParam::getDDDArray	(	const std::string &	str,
		const DDsvalues_type &	sv
	)

private

Definition at line 417 of file HFShowerParam.cc.

References DDfetch(), DDValue::doubles(), edm::hlt::Exception, LogDebug, and relativeConstraints::value.

Referenced by HFShowerParam().

 {
 #ifdef DebugLog
   LogDebug("HFShower") << "HFShowerParam:getDDDArray called for " << str;
 #endif
   DDValue value(str);
   if (DDfetch(&sv,value))
   {
 #ifdef DebugLog
     LogDebug("HFShower") << value;
 #endif
     const std::vector<double> & fvec = value.doubles();
     int nval = fvec.size();
     if (nval < 2) {
       edm::LogError("HFShower") << "HFShowerParam : # of " << str 
                                 << " bins " << nval << " < 2 ==> illegal";
       throw cms::Exception("Unknown", "HFShowerParam") << "nval < 2 for array "
                                                        << str << "\n";
     }
     return fvec;
   } else {
     edm::LogError("HFShower") << "HFShowerParam : cannot get array " << str;
     throw cms::Exception("Unknown", "HFShowerParam")  << "cannot get array "
                                                       << str << "\n";
   }
 }

std::vector< HFShowerParam::Hit > HFShowerParam::getHits	(	G4Step *	aStep,
		double	weight
	)

Definition at line 131 of file HFShowerParam.cc.

References funct::abs(), aperture, applyFidCut, attLMeanInv, gather_cfg::cout, HFShowerParam::Hit::depth, HFShowerParam::Hit::edep, edMin, em_2d_1, em_2d_2, em_lateral_1, em_lateral_2, em_long_1, em_long_2, em_long_gflash, em_long_sl, emPDG, epPDG, create_public_lumi_plots::exp, fibre, fillHisto, gammaPDG, HFShowerLibrary::getHits(), GetVolume(), gflash, HFGflash::gfParameterization(), gpar, hzvem, hzvhad, i, cuy::ii, GetRecoTauVFromDQM_MC_cff::kk, convertSQLiteXML::ok, onlyLong, parametrizeLast, getHLTPrescaleColumns::path, pePerGeV, HFFibreFiducial::PMTNumber(), HFShowerParam::Hit::position, position, funct::pow(), diffTwoXMLs::r1, diffTwoXMLs::r2, ref_index, alignCSCRings::s, showerLibrary, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, HFShowerParam::Hit::time, cond::rpcobgas::time, trackEM, HFFibre::tShift(), histoStyle::weight, detailsBasic3DVector::z, and HFFibre::zShift().

Referenced by HCalSD::getFromParam().

                                               {
   G4StepPoint * preStepPoint  = aStep->GetPreStepPoint(); 
   G4Track *     track    = aStep->GetTrack();   
   G4ThreeVector hitPoint = preStepPoint->GetPosition();   
   G4int         particleCode = track->GetDefinition()->GetPDGEncoding();
   double        zv = std::abs(hitPoint.z()) - gpar[4] - 0.5*gpar[1];
   G4ThreeVector localPoint = G4ThreeVector(hitPoint.x(),hitPoint.y(),zv);
 
   double pin    = (preStepPoint->GetTotalEnergy())/GeV;
   double zint   = hitPoint.z(); 
   double zz     = std::abs(zint) - gpar[4];
 
 #ifdef DebugLog
   edm::LogInfo("HFShower") << "HFShowerParam: getHits " 
                            << track->GetDefinition()->GetParticleName()
                            << " of energy " << pin << " GeV" 
                            << " Pos x,y,z = " << hitPoint.x() << "," 
                            << hitPoint.y() << "," << zint << " (" << zz << ","
                            << localPoint.z() << ", " 
                            << (localPoint.z()+0.5*gpar[1]) << ") Local " 
                            << localPoint;
 #endif
   std::vector<HFShowerParam::Hit> hits;
   HFShowerParam::Hit hit;
   hit.position = hitPoint;
 
   // look for other charged particles
   bool   other = false;
   double pBeta = track->GetDynamicParticle()->GetTotalMomentum() / track->GetDynamicParticle()->GetTotalEnergy();
   double dirz  = (track->GetDynamicParticle()->GetMomentumDirection()).z();
   if (hitPoint.z() < 0) dirz *= -1.;
 #ifdef DebugLog
   edm::LogInfo("HFShower") << "HFShowerParam: getHits Momentum " 
                <<track->GetDynamicParticle()->GetMomentumDirection()
                << " HitPoint " << hitPoint << " dirz " << dirz;
 #endif  
   if (particleCode != emPDG && particleCode != epPDG && particleCode != gammaPDG ) {
     if (track->GetDefinition()->GetPDGCharge() != 0 && pBeta > (1/ref_index) &&
         aStep->GetTotalEnergyDeposit() > 0) other = true;
   }
 
   // take only e+-/gamma/or special particles
   if (particleCode == emPDG || particleCode == epPDG ||
       particleCode == gammaPDG || other) {
     // Leave out the last part
     double edep = 0.;
     bool   kill = false;
     if ((!trackEM) && ((zz<(gpar[1]-gpar[2])) || parametrizeLast) && (!other)){
       edep = pin;
       kill = true;
     } else if ((track->GetDefinition()->GetPDGCharge() != 0) && 
            (pBeta > (1/ref_index)) && (dirz > aperture)) {
       edep = (aStep->GetTotalEnergyDeposit())/GeV;
     }
     std::string path = "ShowerLibrary";
 #ifdef DebugLog
     edm::LogInfo("HFShower") << "HFShowerParam: getHits edep = " << edep
                  << " weight " << weight << " final " <<edep*weight
                              << ", Kill = " << kill << ", pin = " << pin 
                              << ", edMin = " << edMin << " Other " << other;
 #endif
     edep *= weight;
     if (edep > 0) {
       if ((showerLibrary || gflash) && kill && pin > edMin && (!other)) {
         if (showerLibrary) {
           std::vector<HFShowerLibrary::Hit> hitSL = showerLibrary->getHits(aStep,kill,weight,onlyLong);
           for (unsigned int i=0; i<hitSL.size(); i++) {
             bool ok = true;
 #ifdef DebugLog
             edm::LogInfo("HFShower") << "HFShowerParam: getHits applyFidCut = " << applyFidCut;
 #endif
             if (applyFidCut) { // @@ For showerlibrary no z-cut for Short (no z)
               int npmt = HFFibreFiducial:: PMTNumber(hitSL[i].position);
               if (npmt <= 0) ok = false;
             } 
             if (ok) {
               hit.position = hitSL[i].position;
               hit.depth    = hitSL[i].depth;
               hit.time     = hitSL[i].time;
               hit.edep     = 1;
               hits.push_back(hit);
 #ifdef plotDebug
               if (fillHisto) {
         double zv  = std::abs(hit.position.z()) - gpar[4];
                 hzvem->Fill(zv);
                 em_long_sl->Fill(hit.position.z()/cm);
                 double sq = sqrt(pow(hit.position.x()/cm,2)+pow(hit.position.y()/cm,2));
                 double zp = hit.position.z()/cm;
                 if (hit.depth == 1) {
                   em_2d_1->Fill(zp, sq);  
                   em_lateral_1->Fill(sq);
                   em_long_1->Fill(zp);
                 } else if (hit.depth == 2) {
                   em_2d_2->Fill(zp, sq);
                   em_lateral_2->Fill(sq);
                   em_long_2->Fill(zp);
                 }
               }
 #endif
 #ifdef DebugLog
               edm::LogInfo("HFShower") << "HFShowerParam: Hit at depth " 
                                        << hit.depth << " with edep " << hit.edep
                                        << " Time " << hit.time;
 #endif
             }
           }
         } else { // GFlash clusters with known z
           std::vector<HFGflash::Hit>hitSL=gflash->gfParameterization(aStep,kill, onlyLong);
           for (unsigned int i=0; i<hitSL.size(); ++i) {
             bool ok = true;
             G4ThreeVector pe_effect(hitSL[i].position.x(), hitSL[i].position.y(),
                                     hitSL[i].position.z());
             double zv  = std::abs(pe_effect.z()) - gpar[4];
             //depth
             int depth    = 1;
             int npmt     = 0;
             if (zv < 0. || zv > gpar[1]) {
 #ifdef mkdebug
               std::cout<<"-#Zcut-HFShowerParam::getHits:z="<<zv<<",m="<<gpar[1]<<std::endl;
 #endif
               ok = false;
             }
             if (ok && applyFidCut) {
               npmt = HFFibreFiducial:: PMTNumber(pe_effect);
 #ifdef DebugLog
               edm::LogInfo("HFShower") << "HFShowerParam::getHits:#PMT= "
                                        << npmt << ",z = " << zv;
 #endif
               if (npmt <= 0) {
 #ifdef DebugLog
                 edm::LogInfo("HFShower") << "-#PMT=0 cut-HFShowerParam::"
                      << "getHits: npmt = " << npmt;
 #endif
                 ok = false;
               } else if (npmt > 24) { // a short fibre
                 if    (zv > gpar[0]) {
           depth = 2; 
                 } else {
 #ifdef DebugLog
                   edm::LogInfo("HFShower") << "-SHORT cut-HFShowerParam::"
                        << "getHits:zMin=" << gpar[0];
 #endif
                   ok = false;
                 }
               }
 #ifdef DebugLog
               edm::LogInfo("HFShower") << "HFShowerParam: npmt " << npmt 
                                        << " zv " << std::abs(pe_effect.z()) 
                                        << ":" << gpar[4] << ":" << zv << ":" 
                                        << gpar[0] << " ok " << ok << " depth "
                        << depth;
 #endif
             } else {
               if (G4UniformRand() > 0.5) depth = 2;
               if (depth == 2 && zv < gpar[0]) ok = false;
             }
             //attenuation
             double dist = fibre->zShift(localPoint,depth,0); // distance to PMT
             double r1   = G4UniformRand();
 #ifdef DebugLog
             edm::LogInfo("HFShower") << "HFShowerParam:Distance to PMT (" <<npmt
                                      << ") " << dist << ", exclusion flag " 
                                      << (r1 > exp(-attLMeanInv*zv));
 #endif
             if (r1 > exp(-attLMeanInv*dist)) ok = false;
             if (ok) {
           double r2   = G4UniformRand();
 #ifdef DebugLog
           edm::LogInfo("HFShower") << "HFShowerParam:Extra exclusion "
                        << r2 << ">" << weight << " "
                        << (r2 > weight);
 #endif
           if (r2 < weight) {
         double time = fibre->tShift(localPoint,depth,0);
 
         hit.position = hitSL[i].position;
         hit.depth    = depth;
         hit.time     = time + hitSL[i].time;
         hit.edep     = 1;
         hits.push_back(hit);
 #ifdef plotDebug
         if (fillHisto) {
           em_long_gflash->Fill(pe_effect.z()/cm, hitSL[i].edep);
           hzvem->Fill(zv);
           double sq = sqrt(pow(hit.position.x()/cm,2)+pow(hit.position.y()/cm,2));
           double zp = hit.position.z()/cm;
           if (hit.depth == 1) {
             em_2d_1->Fill(zp, sq);
             em_lateral_1->Fill(s);
             em_long_1->Fill(zp);
           } else if (hit.depth == 2) {
             em_2d_2->Fill(zp, sq);
             em_lateral_2->Fill(sq);
             em_long_2->Fill(zp);
           }
         }
 #endif
 #ifdef DebugLog
         edm::LogInfo("HFShower") << "HFShowerParam: Hit at depth " 
                      << hit.depth << " with edep " 
                      << hit.edep << " Time "  << hit.time;
 #endif
           }
             }
           }
         }
       } else {
         path          = "Rest";
         edep         *= pePerGeV;
         double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
         double time   = fibre->tShift(localPoint,1,0); // remaining part
     bool ok = true;
     if (applyFidCut) { // @@ For showerlibrary no z-cut for Short (no z)
       int npmt = HFFibreFiducial:: PMTNumber(hitPoint);
       if (npmt <= 0) ok = false;
     } 
 #ifdef DebugLog
     edm::LogInfo("HFShower") << "HFShowerParam: getHits hitPoint " << hitPoint << " flag " << ok;
 #endif
     if (ok) {
       hit.depth     = 1;
       hit.time      = tSlice+time;
       hit.edep      = edep;
       hits.push_back(hit);
 #ifdef DebugLog
       edm::LogInfo("HFShower") << "HFShowerParam: Hit at depth 1 with edep "
                    << edep << " Time " << tSlice << ":" << time
                    << ":" << hit.time;
 #endif
 #ifdef plotDebug
       double zv = std::abs(hitPoint.z()) - gpar[4];
       if (fillHisto) {
         hzvhad->Fill(zv);
       }
 #endif
       if (zz >= gpar[0]) {
         time      = fibre->tShift(localPoint,2,0);
         hit.depth = 2;
         hit.time  = tSlice+time;
         hits.push_back(hit);
 #ifdef DebugLog
         edm::LogInfo("HFShower") <<"HFShowerParam: Hit at depth 2 with edep "
                      << edep << " Time " << tSlice << ":" << time
                      << hit.time;
 #endif
 #ifdef plotDebug
         if (fillHisto) {
           hzvhad->Fill(zv);
         }
 #endif
       }
     }
       }
 #ifdef DebugLog
       for (unsigned int ii=0; ii<hits.size(); ++ii) {
         double zv = std::abs(hits[ii].position.z());
         if (zv > 12790) edm::LogInfo("HFShower")<< "HFShowerParam: Abnormal hit along " 
                                                 << path << " in " 
                         << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
                                                 << " at " << hits[ii].position << " zz " 
                                                 << zv << " Edep " << edep << " due to " 
                                                 <<track->GetDefinition()->GetParticleName()
                                                 << " time " << hit.time;
       }
 #endif
       if (kill) {
         track->SetTrackStatus(fStopAndKill);
         G4TrackVector tv = *(aStep->GetSecondary());
         for (unsigned int kk=0; kk<tv.size(); ++kk) {
           if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume())
         tv[kk]->SetTrackStatus(fStopAndKill);
         }
       }
 #ifdef DebugLog
       edm::LogInfo("HFShower") << "HFShowerParam: getHits kill (" << kill
                                << ") track " << track->GetTrackID() 
                                << " at " << hitPoint
                                << " and deposit " << edep << " " << hits.size()
                                << " times" << " ZZ " << zz << " " << gpar[0];
 #endif
     }
   }
   return hits;
 }

void HFShowerParam::initRun ( G4ParticleTable * theParticleTable )

Definition at line 120 of file HFShowerParam.cc.

References emPDG, epPDG, gammaPDG, HFShowerLibrary::initRun(), and showerLibrary.

Referenced by HCalSD::initRun().

                                                               {
   emPDG = theParticleTable->FindParticle("e-")->GetPDGEncoding();
   epPDG = theParticleTable->FindParticle("e+")->GetPDGEncoding();
   gammaPDG = theParticleTable->FindParticle("gamma")->GetPDGEncoding();
 #ifdef DebugLog
   edm::LogInfo("HFShower") << "HFShowerParam: Particle code for e- = " << emPDG
                            << " for e+ = " << epPDG << " for gamma = " << gammaPDG;
 #endif
   if (showerLibrary) showerLibrary->initRun(theParticleTable);
 }