#include <CosmicGenFilterHelix.h>

Inheritance diagram for CosmicGenFilterHelix:

Public Member Functions
void	beginJob () override

	CosmicGenFilterHelix (const edm::ParameterSet &config)

void	endJob () override

bool	filter (edm::Event &event, const edm::EventSetup &eventSetup) override

	~CosmicGenFilterHelix () override

Public Member Functions inherited from edm::one::EDFilter< edm::one::SharedResources >
	EDFilter ()=default

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

bool	hasAbilityToProduceInLumis () const final

bool	hasAbilityToProduceInRuns () const final

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

Public Member Functions inherited from edm::one::EDFilterBase
	EDFilterBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDFilterBase () override

Public Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

std::vector< edm::ProductResolverIndex > const &	indiciesForPutProducts (BranchType iBranchType) const

	ProducerBase ()

std::vector< edm::ProductResolverIndex > const &	putTokenIndexToProductResolverIndex () const

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

void	resolvePutIndicies (BranchType iBranchType, ModuleToResolverIndicies const &iIndicies, std::string const &moduleLabel)

	~ProducerBase () noexcept(false) override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
bool	charge (int id, int &charge) const
	true if ID selected, return by value its charge More...

void	createHistsEnd (const char *dirName, TObjArray &hists)

void	createHistsStart (const char *dirName, TObjArray &hists)
	for final statistics: events with track reaching target More...

bool	equidistLogBins (double *bins, int nBins, double first, double last) const

const MagneticField *	getMagneticField (const edm::EventSetup &setup) const
	provide magnetic field from Event Setup More...

const Propagator *	getPropagator (const edm::EventSetup &setup) const

void	monitorEnd (const GlobalPoint &endVert, const GlobalVector &endMom, const GlobalPoint &vert, const GlobalVector &mom, double path, TObjArray &hists)

void	monitorStart (const GlobalPoint &vert, const GlobalVector &mom, int charge, TObjArray &hists)

bool	propagateToCutCylinder (const GlobalPoint &vertStart, const GlobalVector &momStart, int charge, const MagneticField field, const Propagator propagator)
	actually propagate to the defined cylinder More...

Private Attributes
const std::vector< int >	theCharges
	requested Ids More...

const bool	theDoMonitor

TObjArray	theHistsAfter
	hists of properties from generator More...

TObjArray	theHistsBefore
	whether or not to fill monitor hists (needs TFileService) More...

const std::vector< int >	theIds

const double	theMinP2

const double	theMinPt2
	minimal momentum^2 after propagation to cylinder More...

unsigned int	theNumPass
	for final statistics: all seen events More...

unsigned int	theNumTotal
	plane closing cylinder at 'positive' side More...

const std::string	thePropagatorName
	charges, parallel to theIds More...

edm::EDGetTokenT< edm::HepMCProduct >	theSrcToken

Cylinder::ConstCylinderPointer	theTargetCylinder
	minimal transverse^2 momentum after propagation to cylinder More...

Plane::ConstPlanePointer	theTargetPlaneMax
	plane closing cylinder at 'negative' side More...

Plane::ConstPlanePointer	theTargetPlaneMin
	target cylinder, around z-axis More...

Additional Inherited Members
Public Types inherited from edm::one::EDFilterBase
typedef EDFilterBase	ModuleType

Public Types inherited from edm::ProducerBase
using	ModuleToResolverIndicies = std::unordered_multimap< std::string, std::tuple< edm::TypeID const , const char , edm::ProductResolverIndex >>

typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from edm::one::EDFilterBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

-*- C++ -*-

Package: GeneratorInterface/GenFilters

Description: Event filter for generated particles reaching a certain cylinder surface (around z-axis).

Implementation: Assumes particles coming from outside of defined cylinder, but might work also otherwise. Uses SteppingHelixPropagator and IdealMagneticFieldRecord.

Original Author: Gero FLUCKE Created: Mon Mar 5 16:32:01 CET 2007

Definition at line 40 of file CosmicGenFilterHelix.h.

Constructor & Destructor Documentation

CosmicGenFilterHelix::CosmicGenFilterHelix ( const edm::ParameterSet & config )

explicit

Definition at line 42 of file CosmicGenFilterHelix.cc.

References Plane::build(), Cylinder::build(), Exception, edm::ParameterSet::getParameter(), TFileService::kSharedResource, pvSelector_cfi::maxZ, HiClusterCompatibility_cfi::minZ, TCMET_cfi::radius, theCharges, theIds, theSrcToken, theTargetCylinder, theTargetPlaneMax, and theTargetPlaneMin.

   : theIds(cfg.getParameter<std::vector<int> >("pdgIds")),
     theCharges(cfg.getParameter<std::vector<int> >("charges")),
     thePropagatorName(cfg.getParameter<std::string>("propagator")),
     theMinP2(cfg.getParameter<double>("minP")*cfg.getParameter<double>("minP")),
     theMinPt2(cfg.getParameter<double>("minPt")*cfg.getParameter<double>("minPt")),
     theDoMonitor(cfg.getUntrackedParameter<bool>("doMonitor"))
 {
   theSrcToken = consumes<edm::HepMCProduct>(cfg.getParameter<edm::InputTag>("src"));
 
   usesResource(TFileService::kSharedResource);
 
   if (theIds.size() != theCharges.size()) {
     throw cms::Exception("BadConfig") << "CosmicGenFilterHelix: "
                       << "'pdgIds' and 'charges' need same length.";
   }
   Surface::Scalar radius = cfg.getParameter<double>("radius");
   Surface::Scalar maxZ   = cfg.getParameter<double>("maxZ");
   Surface::Scalar minZ   = cfg.getParameter<double>("minZ");
   
   if (maxZ < minZ) {
     throw cms::Exception("BadConfig") << "CosmicGenFilterHelix: maxZ (" << maxZ
                       << ") smaller than minZ (" << minZ << ").";
   }
   
   const Surface::RotationType dummyRot;
   theTargetCylinder = Cylinder::build(Surface::PositionType(0.,0.,0.), dummyRot, radius);
   theTargetPlaneMin = Plane::build(Surface::PositionType(0.,0.,minZ), dummyRot);
   theTargetPlaneMax = Plane::build(Surface::PositionType(0.,0.,maxZ), dummyRot);
 }

CosmicGenFilterHelix::~CosmicGenFilterHelix ( )

override

Definition at line 74 of file CosmicGenFilterHelix.cc.

75 {

76 }

Member Function Documentation

void CosmicGenFilterHelix::beginJob ( void )

overridevirtual

Reimplemented from edm::one::EDFilterBase.

Definition at line 165 of file CosmicGenFilterHelix.cc.

References createHistsEnd(), createHistsStart(), theDoMonitor, theHistsAfter, theHistsBefore, theNumPass, and theNumTotal.

 {
   if (theDoMonitor) {
     this->createHistsStart("start", theHistsBefore);
     this->createHistsStart("startAfter", theHistsAfter);
     // must be after the line above: hist indices are static in monitorStart(...)
     this->createHistsEnd("end", theHistsAfter);
   }
 
   theNumTotal = theNumPass = 0;
 }

bool CosmicGenFilterHelix::charge	(	int	id,
		int &	charge
	)		const

private

true if ID selected, return by value its charge

Definition at line 312 of file CosmicGenFilterHelix.cc.

References end, mps_fire::i, theCharges, and theIds.

Referenced by filter().

 {
   std::vector<int>::const_iterator iC = theCharges.begin();
   for (std::vector<int>::const_iterator i = theIds.begin(), end = theIds.end(); i != end;
        ++i,++iC) {
     if (*i == id) {
       charge = *iC;
       return true;
     }
   }
 
   return false;
 }

void CosmicGenFilterHelix::createHistsEnd	(	const char *	dirName,
		TObjArray &	hists
	)

private

Definition at line 224 of file CosmicGenFilterHelix.cc.

References Abs(), equidistLogBins(), pvSelector_cfi::maxZ, HiClusterCompatibility_cfi::minZ, TFileService::mkdir(), Pi, theTargetCylinder, theTargetPlaneMax, and theTargetPlaneMin.

Referenced by beginJob().

 {
   edm::Service<TFileService> fs;
   TFileDirectory fd(fs->mkdir(dirName, dirName));
 
   const int kNumBins = 50;
   double pBinsLog[kNumBins+1] = {0.}; // fully initialised with 0.
   this->equidistLogBins(pBinsLog, kNumBins, 1., 4000.);
 
   // take care: hist names must differ from those in createHistsStart!
   hists.Add(fd.make<TH1F>("pathEnd", "path until cylinder;s [cm]", 100, 0., 2000.));
   hists.Add(fd.make<TH1F>("momEnd", "|p_{end}|;p [GeV]", 100, 0., 1000.));
   hists.Add(fd.make<TH1F>("momEndLog", "|p_{end}|;p [GeV]", kNumBins, pBinsLog));
   hists.Add(fd.make<TH1F>("ptEnd", "p_{t} (end);p_{t} [GeV]", 100, 0., 750.));
   hists.Add(fd.make<TH1F>("ptEndLog", "p_{t} (end);p_{t} [GeV]", kNumBins, pBinsLog));
   hists.Add(fd.make<TH1F>("phiXzEnd", "#phi_{xz} (end);#phi_{xz}", 100,-TMath::Pi(),TMath::Pi()));
   hists.Add(fd.make<TH1F>("thetaYEnd","#theta_{y} (end);#theta_{y}", 100, 0., TMath::Pi()));
   
   hists.Add(fd.make<TH1F>("momStartEnd", "|p_{start}|-|p_{end}|;#Deltap [GeV]",100,0.,15.));
   hists.Add(fd.make<TH1F>("momStartEndRel", "(p_{start}-p_{end})/p_{start};#Deltap_{rel}",
                           100,.0,1.));
   hists.Add(fd.make<TH1F>("phiXzStartEnd", "#phi_{xz,start}-#phi_{xz,end};#Delta#phi_{xz}",
                           100,-1.,1.));
   hists.Add(fd.make<TH1F>("thetaYStartEnd","#theta_{y,start}-#theta_{y,end};#Delta#theta_{y}",
                           100,-1.,1.));
   
   hists.Add(fd.make<TH2F>("phiXzStartVsEnd",
                           "#phi_{xz} start vs end;#phi_{xz}^{end};#phi_{xz}^{start}",
                           50, -TMath::Pi(), TMath::Pi(), 50, -TMath::Pi(), TMath::Pi()));
   hists.Add(fd.make<TH2F>("thetaYStartVsEnd",
                           "#theta_{y} start vs end;#theta_{y}^{end};#theta_{y}^{start}",
                           50, 0., TMath::Pi(), 50, 0., TMath::Pi()/2.));
   
   hists.Add(fd.make<TH2F>("momStartEndRelVsZ",
                           "(p_{start}-p_{end})/p_{start} vs z_{start};z [cm];#Deltap_{rel}",
                           50, -1600., 1600., 50,.0,.8));
   hists.Add(fd.make<TH2F>("phiXzStartEndVsZ", 
                           "#phi_{xz,start}-#phi_{xz,end} vs z_{start};z [cm];#Delta#phi_{xz}",
                           50, -1600., 1600., 50,-1., 1.));
   hists.Add(fd.make<TH2F>("thetaYStartEndVsZ",
                           "#theta_{y,start}-#theta_{y,end} vs z_{start};z [cm];#Delta#theta_{y}",
                           50, -1600., 1600., 50,-.5,.5));
   hists.Add(fd.make<TH2F>("momStartEndRelVsP",
                           "(p_{start}-p_{end})/p_{start} vs p_{start};p [GeV];#Deltap_{rel}",
                           kNumBins, pBinsLog, 50, .0, .8));
   hists.Add(fd.make<TH2F>("phiXzStartEndVsP", 
                           "#phi_{xz,start}-#phi_{xz,end} vs |p|_{start};p [GeV];#Delta#phi_{xz}",
                           kNumBins, pBinsLog, 100,-1.5, 1.5));
   hists.Add(fd.make<TH2F>("thetaYStartEndVsP",
                           "#theta_{y,start}-#theta_{y,end} vs |p|_{start};p [GeV];#Delta#theta_{y}",
                           kNumBins, pBinsLog, 100,-1.,1.));
   
   const double maxR = theTargetCylinder->radius() * 1.1;
   hists.Add(fd.make<TH1F>("radiusEnd", "end radius;r [cm]", 100, 0., maxR));
   double minZ = theTargetPlaneMin->position().z();
   minZ -= TMath::Abs(minZ) * 0.1;
   double maxZ = theTargetPlaneMax->position().z();
   maxZ += TMath::Abs(maxZ) * 0.1;
   hists.Add(fd.make<TH1F>("zEnd", "end z;z [cm]", 100, minZ, maxZ));
   hists.Add(fd.make<TH1F>("zDiff", "z_{start}-z_{end};#Deltaz [cm]", 100, -1000., 1000.));
   hists.Add(fd.make<TH2F>("xyPlaneEnd", "end xy;x [cm];y [cm]", 100, -maxR, maxR, 100,-maxR,maxR));
   
   hists.Add(fd.make<TH2F>("rzPlaneEnd", "end rz (y<0 #Rightarrow r_{#pm}=-r);z [cm];r_{#pm} [cm]",
                           50, minZ, maxZ, 50, -maxR, maxR));
   hists.Add(fd.make<TH2F>("thetaVsZend", "#theta vs z (end);z [cm];#theta",
                           50, minZ, maxZ, 50, 0., TMath::Pi()));
 }

void CosmicGenFilterHelix::createHistsStart	(	const char *	dirName,
		TObjArray &	hists
	)

private

for final statistics: events with track reaching target

Definition at line 178 of file CosmicGenFilterHelix.cc.

References equidistLogBins(), TFileService::mkdir(), and Pi.

Referenced by beginJob().

 {
   edm::Service<TFileService> fs;
   TFileDirectory fd(fs->mkdir(dirName, dirName));
   
   hists.Add(fd.make<TH1F>("momentumP", "|p(#mu^{+})| (start);|p| [GeV]",100, 0., 1000.));
   hists.Add(fd.make<TH1F>("momentumM", "|p(#mu^{-})| (start);|p| [GeV]",100, 0., 1000.));
   hists.Add(fd.make<TH1F>("momentum2", "|p(#mu)| (start);|p| [GeV]",100, 0., 25.));
   const int kNumBins = 50;
   double pBinsLog[kNumBins+1] = {0.}; // fully initialised with 0.
   this->equidistLogBins(pBinsLog, kNumBins, 1., 4000.);
   hists.Add(fd.make<TH1F>("momentumLog", "|p(#mu)| (start);|p| [GeV]", kNumBins, pBinsLog));
   hists.Add(fd.make<TH1F>("phi", "start p_{#phi(#mu)};#phi", 100, -TMath::Pi(), TMath::Pi()));
   hists.Add(fd.make<TH1F>("cosPhi", "cos(p_{#phi(#mu)}) (start);cos(#phi)", 100, -1., 1.));
   hists.Add(fd.make<TH1F>("phiXz", "start p_{#phi_{xz}(#mu)};#phi_{xz}",
                           100, -TMath::Pi(), TMath::Pi()));
   hists.Add(fd.make<TH1F>("theta", "#theta(#mu) (start);#theta", 100, 0., TMath::Pi()));
   hists.Add(fd.make<TH1F>("thetaY", "#theta_{y}(#mu) (start);#theta_{y}", 100,0.,TMath::Pi()/2.));
   
   hists.Add(fd.make<TH2F>("momVsPhi", "|p(#mu)| vs #phi (start);#phi;|p| [GeV]",
                           50, -TMath::Pi(), TMath::Pi(), 50, 1.5, 1000.));
   hists.Add(fd.make<TH2F>("momVsTheta", "|p(#mu)| vs #theta (start);#theta;|p| [GeV]",
                           50, 0., TMath::Pi(), 50, 1.5, 1000.));
   hists.Add(fd.make<TH2F>("momVsThetaY", "|p(#mu)| vs #theta_{y} (start);#theta_{y};|p| [GeV]",
                           50, 0., TMath::Pi()/2., 50, 1.5, 1000.));
   hists.Add(fd.make<TH2F>("momVsZ", "|p(#mu)| vs z (start);z [cm];|p| [GeV]",
                           50, -1600., 1600., 50, 1.5, 1000.));
   hists.Add(fd.make<TH2F>("thetaVsZ", "#theta vs z (start);z [cm];#theta",
                           50, -1600., 1600., 50, 0., TMath::Pi()));
   hists.Add(fd.make<TH2F>("thetaYvsZ", "#theta_{y} vs z (start);z [cm];#theta",
                           50, -1600., 1600., 50, 0., TMath::PiOver2()));
   hists.Add(fd.make<TH2F>("yVsThetaY", "#theta_{y}(#mu) vs y (start);#theta_{y};y [cm]",
                           50, 0., TMath::Pi()/2., 50, -1000., 1000.));
   hists.Add(fd.make<TH2F>("yVsThetaYnoR", "#theta_{y}(#mu) vs y (start, barrel);#theta_{y};y [cm]",
                           50, 0., TMath::Pi()/2., 50, -1000., 1000.));
   
   hists.Add(fd.make<TH1F>("radius", "start radius;r [cm]", 100, 0., 1000.));
   hists.Add(fd.make<TH1F>("z", "start z;z [cm]", 100, -1600., 1600.));
   hists.Add(fd.make<TH2F>("xyPlane", "start xy;x [cm];y [cm]", 50, -1000., 1000.,
                           50, -1000., 1000.));
   hists.Add(fd.make<TH2F>("rzPlane",
                           "start rz (y < 0 #Rightarrow r_{#pm} = -r);z [cm];r_{#pm} [cm]",
                           50, -1600., 1600., 50, -1000., 1000.));
 }

void CosmicGenFilterHelix::endJob ( void )

overridevirtual

Reimplemented from edm::one::EDFilterBase.

Definition at line 293 of file CosmicGenFilterHelix.cc.

References mps_splice::line, theMinP2, theMinPt2, theNumPass, theNumTotal, theTargetCylinder, theTargetPlaneMax, and theTargetPlaneMin.

 {
   const char *border = "////////////////////////////////////////////////////////";
   const char *line = "\n// ";
   edm::LogInfo("Filter") << "@SUB=CosmicGenFilterHelix::endJob"
              << border << line
                          << theNumPass << " events out of " << theNumTotal
                          << ", i.e. " << theNumPass*100./theNumTotal << "%, "
              << "reached target cylinder," 
              << line << "defined by r < " 
              << theTargetCylinder->radius() << " cm and " 
              << theTargetPlaneMin->position().z() << " < z < " 
              << theTargetPlaneMax->position().z() << " cm."
              << line << "Minimal required (transverse) momentum was "
              << TMath::Sqrt(theMinP2) << " (" << TMath::Sqrt(theMinPt2) << ") GeV."
              << "\n" << border;
 }

bool CosmicGenFilterHelix::equidistLogBins	(	double *	bins,
		int	nBins,
		double	first,
		double	last
	)		const

private

Definition at line 489 of file CosmicGenFilterHelix.cc.

References plotBeamSpotDB::first, mps_fire::i, and plotBeamSpotDB::last.

Referenced by createHistsEnd(), and createHistsStart().

 {
   // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last'
   // that are equidistant when viewed in log scale,
   // so 'bins' must have length nBins+1;
   // If 'first', 'last' or 'nBins' are not positive, failure is reported.
 
   if (nBins < 1 || first <= 0. || last <= 0.) return false;
 
   bins[0] = first;
   bins[nBins] = last;
   const double firstLog = TMath::Log10(bins[0]);
   const double lastLog  = TMath::Log10(bins[nBins]);
   for (int i = 1; i < nBins; ++i) {
     bins[i] = TMath::Power(10., firstLog + i*(lastLog-firstLog)/(nBins));
   }
 
   return true;
 }

bool CosmicGenFilterHelix::filter	(	edm::Event &	event,
		const edm::EventSetup &	eventSetup
	)

override

Definition at line 84 of file CosmicGenFilterHelix.cc.

References ecalTB2006H4_GenSimDigiReco_cfg::bField, charge(), edm::Event::getByToken(), edm::HepMCProduct::GetEvent(), getMagneticField(), getPropagator(), monitorStart(), propagateToCutCylinder(), PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi::propagator, mps_fire::result, theDoMonitor, theHistsBefore, theNumPass, theNumTotal, and theSrcToken.

 {
   edm::Handle<edm::HepMCProduct> hepMCEvt;
   iEvent.getByToken(theSrcToken, hepMCEvt);
   const HepMC::GenEvent *mCEvt = hepMCEvt->GetEvent();
   const MagneticField *bField = this->getMagneticField(iSetup); // should be fast (?)
   const Propagator *propagator = this->getPropagator(iSetup);
 
   ++theNumTotal;
   bool result = false;
   for (HepMC::GenEvent::particle_const_iterator iPart = mCEvt->particles_begin(),
      endPart = mCEvt->particles_end(); iPart != endPart; ++iPart) {
     int charge = 0; // there is no method providing charge in GenParticle :-(
     if ((*iPart)->status() != 1) continue; // look only at stable particles
     if (!this->charge((*iPart)->pdg_id(), charge)) continue;
 
     // Get the position and momentum
     const HepMC::ThreeVector hepVertex((*iPart)->production_vertex()->point3d());
     const GlobalPoint vert(hepVertex.x()/10., hepVertex.y()/10., hepVertex.z()/10.); // to cm
     const HepMC::FourVector hepMomentum((*iPart)->momentum());
     const GlobalVector mom(hepMomentum.x(), hepMomentum.y(), hepMomentum.z());
 
     if (theDoMonitor) this->monitorStart(vert, mom, charge, theHistsBefore);
 
     if (this->propagateToCutCylinder(vert, mom, charge, bField, propagator)) {
       result = true;
     }
   }
 
   if (result) ++theNumPass;
   return result;
 }

const MagneticField * CosmicGenFilterHelix::getMagneticField ( const edm::EventSetup & setup ) const

private

provide magnetic field from Event Setup

Definition at line 327 of file CosmicGenFilterHelix.cc.

References fieldHandle, edm::EventSetup::get(), and edm::ESHandle< T >::product().

Referenced by filter().

 {
   edm::ESHandle<MagneticField> fieldHandle;
   setup.get<IdealMagneticFieldRecord>().get(fieldHandle); 
 
   return fieldHandle.product();
 }

const Propagator * CosmicGenFilterHelix::getPropagator ( const edm::EventSetup & setup ) const

private

Definition at line 336 of file CosmicGenFilterHelix.cc.

References edm::EventSetup::get(), edm::ESHandle< T >::product(), and thePropagatorName.

Referenced by filter().

 {
   edm::ESHandle<Propagator> propHandle;
   setup.get<TrackingComponentsRecord>().get(thePropagatorName, propHandle);
 
   const Propagator *prop = propHandle.product();
   if (!dynamic_cast<const SteppingHelixPropagator*>(prop)) {
     edm::LogWarning("BadConfig") << "@SUB=CosmicGenFilterHelix::getPropagator"
                                  << "Not a SteppingHelixPropagator!";
 
   }
   return prop;
 }

void CosmicGenFilterHelix::monitorEnd	(	const GlobalPoint &	endVert,
		const GlobalVector &	endMom,
		const GlobalPoint &	vert,
		const GlobalVector &	mom,
		double	path,
		TObjArray &	hists
	)

private

Definition at line 412 of file CosmicGenFilterHelix.cc.

References HcalObjRepresent::Fill(), PV3DBase< T, PVType, FrameType >::mag(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::theta(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by propagateToCutCylinder().

 {
   const double scalarMomStart = mom.mag();
   const double phiXzStart = TMath::ATan2(mom.z(), mom.x());
   const double thetaYStart = TMath::ATan2(TMath::Sqrt(mom.x()*mom.x()+mom.z()*mom.z()), -mom.y());
   const double scalarMomEnd = endMom.mag();
   const double ptEnd = endMom.perp();
   const double phiXzEnd = TMath::ATan2(endMom.z(), endMom.x());
   const double thetaYEnd = TMath::ATan2(TMath::Sqrt(endMom.x()*endMom.x()+endMom.z()*endMom.z()),
                     -endMom.y());
   const double thetaEnd = endMom.theta();
 
   const double diffMomRel = (scalarMomStart-scalarMomEnd)/scalarMomStart;
 
   const double zEnd = endVert.z();
   const double rEnd = endVert.perp();
   const double diffZ = zEnd - vert.z();
 
   const static int iPathEnd = hists.IndexOf(hists.FindObject("pathEnd"));
   static_cast<TH1*>(hists[iPathEnd])->Fill(path);
   const static int iMomEnd = hists.IndexOf(hists.FindObject("momEnd"));
   static_cast<TH1*>(hists[iMomEnd])->Fill(scalarMomEnd);
   const static int iMomEndLog = hists.IndexOf(hists.FindObject("momEndLog"));
   static_cast<TH1*>(hists[iMomEndLog])->Fill(scalarMomEnd);
   const static int iPtEnd = hists.IndexOf(hists.FindObject("ptEnd"));
   static_cast<TH1*>(hists[iPtEnd])->Fill(ptEnd);
   const static int iPtEndLog = hists.IndexOf(hists.FindObject("ptEndLog"));
   static_cast<TH1*>(hists[iPtEndLog])->Fill(ptEnd);
   const static int iPhiXzEnd = hists.IndexOf(hists.FindObject("phiXzEnd"));
   static_cast<TH1*>(hists[iPhiXzEnd])->Fill(phiXzEnd);
   const static int iThetaYEnd = hists.IndexOf(hists.FindObject("thetaYEnd"));
   static_cast<TH1*>(hists[iThetaYEnd])->Fill(thetaYEnd);
 
   const static int iMomStartEnd = hists.IndexOf(hists.FindObject("momStartEnd"));
   static_cast<TH1*>(hists[iMomStartEnd])->Fill(scalarMomStart-scalarMomEnd);
   const static int iMomStartEndRel = hists.IndexOf(hists.FindObject("momStartEndRel"));
   static_cast<TH1*>(hists[iMomStartEndRel])->Fill(diffMomRel);
   const static int iPhiStartEnd = hists.IndexOf(hists.FindObject("phiXzStartEnd"));
   static_cast<TH1*>(hists[iPhiStartEnd])->Fill(phiXzStart-phiXzEnd);
   const static int iThetaStartEnd = hists.IndexOf(hists.FindObject("thetaYStartEnd"));
   static_cast<TH1*>(hists[iThetaStartEnd])->Fill(thetaYStart-thetaYEnd);
 
   const static int iPhiStartVsEnd = hists.IndexOf(hists.FindObject("phiXzStartVsEnd"));
   static_cast<TH2*>(hists[iPhiStartVsEnd])->Fill(phiXzEnd, phiXzStart);
   const static int iThetaStartVsEnd = hists.IndexOf(hists.FindObject("thetaYStartVsEnd"));
   static_cast<TH2*>(hists[iThetaStartVsEnd])->Fill(thetaYEnd, thetaYStart);
 
   const static int iMomStartEndRelVsZ = hists.IndexOf(hists.FindObject("momStartEndRelVsZ"));
   static_cast<TH2*>(hists[iMomStartEndRelVsZ])->Fill(vert.z(), diffMomRel);
   const static int iPhiStartEndVsZ = hists.IndexOf(hists.FindObject("phiXzStartEndVsZ"));
   static_cast<TH2*>(hists[iPhiStartEndVsZ])->Fill(vert.z(), phiXzStart-phiXzEnd);
   const static int iThetaStartEndVsZ = hists.IndexOf(hists.FindObject("thetaYStartEndVsZ"));
   static_cast<TH2*>(hists[iThetaStartEndVsZ])->Fill(vert.z(), thetaYStart-thetaYEnd);
   const static int iMomStartEndRelVsP = hists.IndexOf(hists.FindObject("momStartEndRelVsP"));
   static_cast<TH2*>(hists[iMomStartEndRelVsP])->Fill(scalarMomStart, diffMomRel);
   const static int iPhiStartEndVsP = hists.IndexOf(hists.FindObject("phiXzStartEndVsP"));
   static_cast<TH2*>(hists[iPhiStartEndVsP])->Fill(scalarMomStart, phiXzStart-phiXzEnd);
   const static int iThetaStartEndVsP = hists.IndexOf(hists.FindObject("thetaYStartEndVsP"));
   static_cast<TH2*>(hists[iThetaStartEndVsP])->Fill(scalarMomStart, thetaYStart-thetaYEnd);
 
   const static int iRadiusEnd = hists.IndexOf(hists.FindObject("radiusEnd"));
   static_cast<TH1*>(hists[iRadiusEnd])->Fill(rEnd);
   const static int iZend = hists.IndexOf(hists.FindObject("zEnd"));
   static_cast<TH1*>(hists[iZend])->Fill(zEnd);
   const static int iZdiff = hists.IndexOf(hists.FindObject("zDiff"));
   static_cast<TH1*>(hists[iZdiff])->Fill(diffZ);
   const static int iXyPlaneEnd = hists.IndexOf(hists.FindObject("xyPlaneEnd"));
   static_cast<TH1*>(hists[iXyPlaneEnd])->Fill(endVert.x(), endVert.y());
   const static int iRzPlaneEnd = hists.IndexOf(hists.FindObject("rzPlaneEnd"));
   static_cast<TH1*>(hists[iRzPlaneEnd])->Fill(zEnd, (endVert.y() > 0 ? rEnd : -rEnd));
   const static int iThetaVsZend = hists.IndexOf(hists.FindObject("thetaVsZend"));
   static_cast<TH2*>(hists[iThetaVsZend])->Fill(zEnd, thetaEnd);
 }

void CosmicGenFilterHelix::monitorStart	(	const GlobalPoint &	vert,
		const GlobalVector &	mom,
		int	charge,
		TObjArray &	hists
	)

private

Definition at line 351 of file CosmicGenFilterHelix.cc.

References HcalObjRepresent::Fill(), PV3DBase< T, PVType, FrameType >::mag(), PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), alignCSCRings::r, PV3DBase< T, PVType, FrameType >::theta(), theta(), DOFs::thetaY, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by filter(), and propagateToCutCylinder().

 {
   const double scalarMom = mom.mag();
   const double phi = mom.phi();
   const double phiXz = TMath::ATan2(mom.z(), mom.x());
   const double theta = mom.theta();
   const double thetaY = TMath::ATan2(TMath::Sqrt(mom.x()*mom.x()+mom.z()*mom.z()), -mom.y());
 
   const double z = vert.z();
   const double r = vert.perp();
 
   const static int iMomP = hists.IndexOf(hists.FindObject("momentumP"));
   const static int iMomM = hists.IndexOf(hists.FindObject("momentumM"));
   if (charge > 0) static_cast<TH1*>(hists[iMomP])->Fill(scalarMom);
   else            static_cast<TH1*>(hists[iMomM])->Fill(scalarMom);
   const static int iMom2 = hists.IndexOf(hists.FindObject("momentum2"));
   static_cast<TH1*>(hists[iMom2])->Fill(scalarMom);
   const static int iMomLog = hists.IndexOf(hists.FindObject("momentumLog"));
   static_cast<TH1*>(hists[iMomLog])->Fill(scalarMom);
   const static int iPhi = hists.IndexOf(hists.FindObject("phi"));
   static_cast<TH1*>(hists[iPhi])->Fill(phi);
   const static int iCosPhi = hists.IndexOf(hists.FindObject("cosPhi"));
   static_cast<TH1*>(hists[iCosPhi])->Fill(TMath::Cos(phi));
   const static int iPhiXz = hists.IndexOf(hists.FindObject("phiXz"));
   static_cast<TH1*>(hists[iPhiXz])->Fill(phiXz);
   const static int iTheta = hists.IndexOf(hists.FindObject("theta"));
   static_cast<TH1*>(hists[iTheta])->Fill(theta);
   const static int iThetaY = hists.IndexOf(hists.FindObject("thetaY"));
   static_cast<TH1*>(hists[iThetaY])->Fill(thetaY);
 
   const static int iMomVsTheta = hists.IndexOf(hists.FindObject("momVsTheta"));
   static_cast<TH2*>(hists[iMomVsTheta])->Fill(theta, scalarMom);
   const static int iMomVsThetaY = hists.IndexOf(hists.FindObject("momVsThetaY"));
   static_cast<TH2*>(hists[iMomVsThetaY])->Fill(thetaY, scalarMom);
   const static int iMomVsPhi = hists.IndexOf(hists.FindObject("momVsPhi"));
   static_cast<TH2*>(hists[iMomVsPhi])->Fill(phi, scalarMom);
   const static int iMomVsZ = hists.IndexOf(hists.FindObject("momVsZ"));
   static_cast<TH2*>(hists[iMomVsZ])->Fill(z, scalarMom);
   const static int iThetaVsZ = hists.IndexOf(hists.FindObject("thetaVsZ"));
   static_cast<TH2*>(hists[iThetaVsZ])->Fill(z, theta);
   const static int iThetaYvsZ = hists.IndexOf(hists.FindObject("thetaYvsZ"));
   static_cast<TH2*>(hists[iThetaYvsZ])->Fill(z, thetaY);
   const static int iYvsThetaY = hists.IndexOf(hists.FindObject("yVsThetaY"));
   static_cast<TH2*>(hists[iYvsThetaY])->Fill(thetaY, vert.y());
   const static int iYvsThetaYnoR = hists.IndexOf(hists.FindObject("yVsThetaYnoR"));
   if (z > -1400. && z < 1400.) {
     static_cast<TH2*>(hists[iYvsThetaYnoR])->Fill(thetaY, vert.y());
   }
 
   const static int iRadius = hists.IndexOf(hists.FindObject("radius"));
   static_cast<TH1*>(hists[iRadius])->Fill(r);
   const static int iZ = hists.IndexOf(hists.FindObject("z"));
   static_cast<TH1*>(hists[iZ])->Fill(z);
   const static int iXy = hists.IndexOf(hists.FindObject("xyPlane"));
   static_cast<TH1*>(hists[iXy])->Fill(vert.x(), vert.y());
   const static int iRz = hists.IndexOf(hists.FindObject("rzPlane"));
   static_cast<TH1*>(hists[iRz])->Fill(z, (vert.y() > 0 ? r : -r));
 }

bool CosmicGenFilterHelix::propagateToCutCylinder	(	const GlobalPoint &	vertStart,
		const GlobalVector &	momStart,
		int	charge,
		const MagneticField *	field,
		const Propagator *	propagator
	)

private

actually propagate to the defined cylinder

Definition at line 118 of file CosmicGenFilterHelix.cc.

References monitorEnd(), monitorStart(), Propagator::propagateWithPath(), mps_fire::result, theDoMonitor, theHistsAfter, theMinP2, theMinPt2, theTargetCylinder, theTargetPlaneMax, and theTargetPlaneMin.

Referenced by filter().

 {
   typedef std::pair<TrajectoryStateOnSurface, double> TsosPath;
 
   const FreeTrajectoryState fts(GlobalTrajectoryParameters(vertStart, momStart, charge, field));
 
   bool result = true;
   TsosPath aTsosPath(propagator->propagateWithPath(fts, *theTargetCylinder));
   if (!aTsosPath.first.isValid()) {
     result = false;
   } else if (aTsosPath.first.globalPosition().z() < theTargetPlaneMin->position().z()) {
     // If on cylinder, but outside minimum z, try minimum z-plane:
     // (Would it be possible to miss radius on plane, but reach cylinder afterwards in z-range?
     //  No, at least not in B-field parallel to z-axis which is cylinder axis.)
     aTsosPath = propagator->propagateWithPath(fts, *theTargetPlaneMin);
     if (!aTsosPath.first.isValid()
     || aTsosPath.first.globalPosition().perp() > theTargetCylinder->radius()) {
       result = false;
     }
   } else if (aTsosPath.first.globalPosition().z() > theTargetPlaneMax->position().z()) {
     // Analog for outside maximum z:
     aTsosPath = propagator->propagateWithPath(fts, *theTargetPlaneMax);
     if (!aTsosPath.first.isValid()
     || aTsosPath.first.globalPosition().perp() > theTargetCylinder->radius()) {
       result = false;
     }
   }
 
   if (result) {
     const GlobalVector momEnd(aTsosPath.first.globalMomentum());
     if (momEnd.perp2() < theMinPt2 || momEnd.mag2() < theMinP2) {
       result = false;
     } else if (theDoMonitor) {
       const GlobalPoint vertEnd(aTsosPath.first.globalPosition());
       this->monitorStart(vertStart, momStart, charge, theHistsAfter);
       this->monitorEnd(vertEnd, momEnd, vertStart, momStart, aTsosPath.second, theHistsAfter);
     }
   }
 
   return result;
 }