#include <Lukas/TestPythiaDecays/plugins/TestPythiaDecays.cc>

Inheritance diagram for TestPythiaDecays:

Public Member Functions
	TestPythiaDecays (const edm::ParameterSet &)

	~TestPythiaDecays () override=default

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

	EDAnalyzer (const EDAnalyzer &)=delete

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

const EDAnalyzer &	operator= (const EDAnalyzer &)=delete

bool	wantsGlobalLuminosityBlocks () const noexcept final

bool	wantsGlobalRuns () const noexcept final

bool	wantsInputProcessBlocks () const noexcept final

bool	wantsProcessBlocks () const noexcept final

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

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const noexcept

bool	wantsStreamRuns () const noexcept

	~EDAnalyzerBase () 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

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

std::vector< ESResolverIndex > const &	esGetTokenIndicesVector (edm::Transition iTrans) const

std::vector< ESRecordIndex > const &	esGetTokenRecordIndicesVector (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::array< std::vector< ModuleDescription const > , NumBranchTypes > &modulesAll, std::vector< ModuleProcessName > &modulesInPreviousProcesses, 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

void	selectInputProcessBlocks (ProductRegistry const &productRegistry, ProcessBlockHelperBase const &processBlockHelperBase)

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

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

void	updateLookup (eventsetup::ESRecordsToProductResolverIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Private Member Functions
void	analyze (const edm::Event &, const edm::EventSetup &) override

Private Attributes
std::map< int, TH1D * >	h_br

std::map< int, TH1D * >	h_br_ref

std::map< int, TH1D * >	h_decayVertexRho

std::map< int, TH1D * >	h_decayVertexZ

std::map< int, TH1D * >	h_mass

std::map< int, TH1D * >	h_mass_ref

std::map< int, TH1D * >	h_originVertexRho

std::map< int, TH1D * >	h_originVertexZ

std::map< int, TH1D * >	h_p

std::map< int, TH1D * >	h_plt

std::map< int, TH1D * >	h_plt_ref

std::map< int, TH1D * >	h_v

std::map< int, std::vector< std::string > >	knownDecayModes

std::string	outputFile

std::vector< int >	pids

Pythia8::Pythia *	pythia

const edm::EDGetTokenT< std::vector< SimTrack > >	tokTrack_

const edm::EDGetTokenT< std::vector< SimVertex > >	tokVertex_

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

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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

template<BranchType B = InEvent>
EDConsumerBaseAdaptor< B >	consumes (edm::InputTag tag) noexcept

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 ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

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

template<Transition Tr = Transition::Event>
constexpr auto	esConsumes ()

template<Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag tag)

template<Transition Tr = Transition::Event>
ESGetTokenGeneric	esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
	Used with EventSetupRecord::doGet. More...

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)

void	resetItemsToGetFrom (BranchType iType)

Detailed Description

Description: [one line class summary]

Implementation: [Notes on implementation]

Definition at line 56 of file TestPythiaDecays.cc.

Constructor & Destructor Documentation

◆ TestPythiaDecays()

TestPythiaDecays::TestPythiaDecays ( const edm::ParameterSet & iConfig )

explicit

Definition at line 100 of file TestPythiaDecays.cc.

References funct::abs(), b, ztail::d, Exception, beamvalidation::exit(), geometryDiff::file, edm::ParameterSet::getParameter(), h_br, h_br_ref, h_decayVertexRho, h_decayVertexZ, h_mass, h_mass_ref, h_originVertexRho, h_originVertexZ, h_p, h_plt, h_plt_ref, h_v, mps_fire::i, knownDecayModes, TFileService::kSharedResource, label, TFileDirectory::make(), SiStripPI::min, outputFile, AlCaHLTBitMon_ParallelJobs::p, pids, dumpMFGeometry_cfg::prod, pythia, jetUpdater_cfi::sort, AlCaHLTBitMon_QueryRunRegistry::string, and w().

     : tokTrack_(consumes<std::vector<SimTrack> >(edm::InputTag("famosSimHits"))),
       tokVertex_(consumes<std::vector<SimVertex> >(edm::InputTag("famosSimHits"))) {
   usesResource(TFileService::kSharedResource);
   // output file
   outputFile = iConfig.getParameter<std::string>("outputFile");
 
   // create pythia8 instance to access particle data
   pythia = new Pythia8::Pythia();
   pythia->init();
   Pythia8::ParticleData pdt = pythia->particleData;
 
   // which particles will we study?
   pids.push_back(15);   // tau
   pids.push_back(211);  // pi+
   pids.push_back(111);  // pi0
   pids.push_back(130);  // K0L
   pids.push_back(321);  // K+
   pids.push_back(323);  // K*(392)
   pids.push_back(411);  // D+
   pids.push_back(521);  // B+
 
   // define histograms
   edm::Service<TFileService> file;
   TFileDirectory observe = file->mkdir("observed");
   TFileDirectory predict = file->mkdir("prediction");
   for (size_t i = 0; i < pids.size(); ++i) {
     int pid = std::abs(pids[i]);
 
     // get particle data
     if (!pdt.isParticle(pid)) {
       edm::LogError("TestPythiaDecays") << "ERROR: BAD PARTICLE, pythia is not aware of pid " << pid;
       throw cms::Exception("Unknown", "FastSim") << "Bad Particle Type " << pid << "\n";
       std::exit(1);
     }
     ParticleDataEntryPtr pd = pdt.particleDataEntryPtr(pid);
 
     // mass histograms
     double m0 = pd->m0();
     double w = pd->mWidth();
     double mmin, mmax;
     if (w == 0) {
       mmin = m0 - m0 / 1000.;
       mmax = m0 + m0 / 1000.;
     } else {
       mmin = m0 - 2 * w;
       mmax = m0 + 2 * w;
     }
     std::stringstream strstr;
     strstr << "mass_" << pid;
     h_mass[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, mmin, mmax);
     h_mass_ref[pid] = predict.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, mmin, mmax);
     if (w == 0)
       h_mass_ref[pid]->Fill(m0);
     else {
       for (int b = 1; b <= h_mass_ref[pid]->GetNbinsX(); ++b) {
         double _val = h_mass_ref[pid]->GetBinCenter(b);
         h_mass_ref[pid]->SetBinContent(b, TMath::BreitWigner(_val, m0, w));
       }
     }
 
     // p histogram
     strstr.str("");
     strstr << "p_" << pid;
     h_p[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 20);
 
     // v histogram
     strstr.str("");
     strstr << "v_" << pid;
     h_v[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 1.);
 
     // ctau histograms
     double ctau0 = pd->tau0() / 10.;
     strstr.str("");
     strstr << "plt_" << pid;
     h_plt[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, std::min(5. * ctau0, 500.));
     h_plt_ref[pid] = predict.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, std::min(5. * ctau0, 500.));
     h_plt_ref[pid]->SetTitle(h_plt_ref[pid]->GetName());
     for (int b = 1; b <= h_plt_ref[pid]->GetNbinsX(); ++b) {
       double _val = h_plt_ref[pid]->GetBinCenter(b);
       h_plt_ref[pid]->SetBinContent(b, TMath::Exp(-_val / ctau0));  //convert mm to cm
     }
 
     // br histograms
     strstr.str("");
     strstr << "br_" << pid;
     h_br[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 0, 0, 0);
     h_br[pid]->SetCanExtend(TH1::kAllAxes);
     h_br_ref[pid] = predict.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 0, 0, 0);
     h_br_ref[pid]->SetCanExtend(TH1::kAllAxes);
     h_br_ref[pid]->SetTitle(h_br_ref[pid]->GetName());
     knownDecayModes[pid] = std::vector<std::string>();
     for (int d = 0; d < pd->sizeChannels(); ++d) {
       Pythia8::DecayChannel& channel = pd->channel(d);
       std::vector<int> prod;
       for (int p = 0; p < channel.multiplicity(); ++p) {
         int pId = abs(channel.product(p));
         // from FastSimulation/Event/src/KineParticleFilter.cc
         bool particleCut =
             (pId > 10 && pId != 12 && pId != 14 && pId != 16 && pId != 18 && pId != 21 && (pId < 23 || pId > 40) &&
              (pId < 81 || pId > 100) && pId != 2101 && pId != 3101 && pId != 3201 && pId != 1103 && pId != 2103 &&
              pId != 2203 && pId != 3103 && pId != 3203 && pId != 3303);
         if (particleCut)
           prod.push_back(abs(channel.product(p)));
       }
       std::sort(prod.begin(), prod.end());
       strstr.str("");
       for (size_t p = 0; p < prod.size(); ++p) {
         strstr << "_" << prod[p];
       }
       std::string label = strstr.str();
       h_br[pid]->Fill(label.c_str(), 0.);
       h_br_ref[pid]->Fill(label.c_str(), channel.bRatio());
       h_br[pid]->SetEntries(0);
       knownDecayModes[pid].push_back(label);
     }
 
     // vertex plots
     strstr.str("");
     strstr << "originVertexRho_" << pid;
     h_originVertexRho[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 200);
     strstr.str("");
     strstr << "originVertexZ_" << pid;
     h_originVertexZ[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 400);
     strstr.str("");
     strstr << "decayVertexRho_" << pid;
     h_decayVertexRho[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 200);
     strstr.str("");
     strstr << "decayVertexZ_" << pid;
     h_decayVertexZ[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 400);
   }
 }

◆ ~TestPythiaDecays()

TestPythiaDecays::~TestPythiaDecays ( )

overridedefault

Member Function Documentation

◆ analyze()

void TestPythiaDecays::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprivatevirtual

Implements edm::one::EDAnalyzerBase.

Definition at line 238 of file TestPythiaDecays.cc.

References funct::abs(), ztail::d, relativeConstraints::empty, mps_fire::end, spr::find(), h_br, h_br_ref, h_decayVertexRho, h_decayVertexZ, h_mass, h_originVertexRho, h_originVertexZ, h_p, h_plt, h_v, iEvent, createfilelist::int, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, dqmiolumiharvest::j, knownDecayModes, label, genParticles_cff::map, EgHLTOffHistBins_cfi::mass, AlCaHLTBitMon_ParallelJobs::p, class-composition::parent, pids, CoreSimVertex::position(), dumpMFGeometry_cfg::prod, jetUpdater_cfi::sort, AlCaHLTBitMon_QueryRunRegistry::string, tokTrack_, tokVertex_, bphysicsOniaDQM_cfi::vertex, and SimTrack::vertIndex().

                                                                                   {
   const edm::Handle<std::vector<SimTrack> >& simtracks = iEvent.getHandle(tokTrack_);
   const edm::Handle<std::vector<SimVertex> > simvertices = iEvent.getHandle(tokVertex_);
 
   // create maps
 
   // initialize
   std::map<size_t, std::vector<size_t> > childMap;  // child indices vs parent index
   std::map<size_t, int> parentMap;                  // parent index vs child index
   for (size_t j = 0; j < simtracks->size(); j++) {
     childMap[j] = std::vector<size_t>();
     parentMap[j] = -1;
   }
 
   // do the mapping
   for (size_t j = 0; j < simtracks->size(); j++) {
     size_t childIndex = j;
     const SimTrack& child = simtracks->at(childIndex);
     if (child.noVertex())
       continue;
     const SimVertex& vertex = simvertices->at(child.vertIndex());
     if (vertex.noParent())
       continue;
     size_t parentIndex = vertex.parentIndex();
     childMap[parentIndex].push_back(childIndex);
     parentMap[childIndex] = int(parentIndex);
   }
 
   for (size_t j = 0; j < simtracks->size(); j++) {
     const SimTrack& parent = simtracks->at(j);
     int pid = abs(parent.type());
     if (std::find(pids.begin(), pids.end(), pid) == pids.end())
       continue;
 
     // fill mass hist
     double mass = parent.momentum().M();
     h_mass[pid]->Fill(mass);
 
     // fill p hist
     h_p[pid]->Fill(parent.momentum().P());
 
     // fill vertex position hist
     if (!parent.noVertex()) {
       const SimVertex& originVertex = simvertices->at(parent.vertIndex());
       h_originVertexRho[pid]->Fill(originVertex.position().Rho());
       h_originVertexZ[pid]->Fill(std::fabs(originVertex.position().Z()));
     }
     if (!childMap[j].empty()) {
       const SimTrack& child = simtracks->at(childMap[j][0]);
       const SimVertex& decayVertex = simvertices->at(child.vertIndex());
       h_decayVertexRho[pid]->Fill(decayVertex.position().Rho());
       h_decayVertexZ[pid]->Fill(std::fabs(decayVertex.position().Z()));
     }
   }
 
   for (std::map<size_t, std::vector<size_t> >::iterator it = childMap.begin(); it != childMap.end(); ++it) {
     // fill ctau hist
     size_t parentIndex = it->first;
     const SimTrack& parent = simtracks->at(parentIndex);
     int pid = abs(parent.type());
     std::vector<size_t>& childIndices = it->second;
     if (childIndices.empty())
       continue;
 
     if (std::find(pids.begin(), pids.end(), pid) == pids.end())
       continue;
 
     const SimVertex& origin_vertex = simvertices->at(parent.vertIndex());
     const SimTrack& child0 = simtracks->at(childIndices[0]);
     const SimVertex& decay_vertex = simvertices->at(child0.vertIndex());
 
     TLorentzVector lv_origin_vertex(origin_vertex.position().X(),
                                     origin_vertex.position().Y(),
                                     origin_vertex.position().Z(),
                                     origin_vertex.position().T());
     TLorentzVector lv_decay_vertex(decay_vertex.position().X(),
                                    decay_vertex.position().Y(),
                                    decay_vertex.position().Z(),
                                    decay_vertex.position().T());
     TLorentzVector lv_dist = lv_decay_vertex - lv_origin_vertex;
     TLorentzVector lv_parent(
         parent.momentum().Px(), parent.momentum().Py(), parent.momentum().Pz(), parent.momentum().E());
     TVector3 boost = lv_parent.BoostVector();
     lv_dist.Boost(-boost);
     h_v[pid]->Fill(boost.Mag());
     double plt = lv_dist.T();
     h_plt[pid]->Fill(plt);
 
     // fill br hist
     std::vector<int> prod;
     for (size_t d = 0; d < childIndices.size(); ++d) {
       prod.push_back(abs(simtracks->at(childIndices[d]).type()));
     }
     std::sort(prod.begin(), prod.end());
     std::stringstream strstr;
     for (size_t p = 0; p < prod.size(); ++p) {
       strstr << "_" << prod[p];
     }
     std::string label = strstr.str();
     if (std::find(knownDecayModes[pid].begin(), knownDecayModes[pid].end(), label) == knownDecayModes[pid].end())
       label = "u" + label;
     h_br[pid]->Fill(label.c_str(), 1.);
     h_br_ref[pid]->Fill(label.c_str(), 0.);  // keep h_br and h_br_ref in sync
   }
 }

◆ fillDescriptions()

void TestPythiaDecays::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

static

Definition at line 345 of file TestPythiaDecays.cc.

References edm::ConfigurationDescriptions::addDefault(), and submitPVResolutionJobs::desc.

                                                                                   {
   //The following says we do not know what parameters are allowed so do no validation
   // Please change this to state exactly what you do use, even if it is no parameters
   edm::ParameterSetDescription desc;
   desc.setUnknown();
   descriptions.addDefault(desc);
 }

Member Data Documentation

◆ h_br

std::map<int, TH1D*> TestPythiaDecays::h_br

private

Definition at line 80 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_br_ref

std::map<int, TH1D*> TestPythiaDecays::h_br_ref

private

Definition at line 81 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_decayVertexRho

std::map<int, TH1D*> TestPythiaDecays::h_decayVertexRho

private

Definition at line 77 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_decayVertexZ

std::map<int, TH1D*> TestPythiaDecays::h_decayVertexZ

private

Definition at line 78 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_mass

std::map<int, TH1D*> TestPythiaDecays::h_mass

private

Definition at line 70 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_mass_ref

std::map<int, TH1D*> TestPythiaDecays::h_mass_ref

private

Definition at line 73 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays().

◆ h_originVertexRho

std::map<int, TH1D*> TestPythiaDecays::h_originVertexRho

private

Definition at line 75 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_originVertexZ

std::map<int, TH1D*> TestPythiaDecays::h_originVertexZ

private

Definition at line 76 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_p

std::map<int, TH1D*> TestPythiaDecays::h_p

private

Definition at line 71 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_plt

std::map<int, TH1D*> TestPythiaDecays::h_plt

private

Definition at line 74 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ h_plt_ref

std::map<int, TH1D*> TestPythiaDecays::h_plt_ref

private

Definition at line 79 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays().

◆ h_v

std::map<int, TH1D*> TestPythiaDecays::h_v

private

Definition at line 72 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ knownDecayModes

std::map<int, std::vector<std::string> > TestPythiaDecays::knownDecayModes

private

Definition at line 83 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ outputFile

std::string TestPythiaDecays::outputFile

private

Definition at line 86 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays().

◆ pids

std::vector<int> TestPythiaDecays::pids

private

Definition at line 69 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

◆ pythia

Pythia8::Pythia* TestPythiaDecays::pythia

private

Definition at line 85 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays().

◆ tokTrack_

const edm::EDGetTokenT<std::vector<SimTrack> > TestPythiaDecays::tokTrack_

private

Definition at line 67 of file TestPythiaDecays.cc.

Referenced by analyze().

◆ tokVertex_

const edm::EDGetTokenT<std::vector<SimVertex> > TestPythiaDecays::tokVertex_

private

Definition at line 68 of file TestPythiaDecays.cc.

Referenced by analyze().

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ TestPythiaDecays()

◆ ~TestPythiaDecays()

Member Function Documentation

◆ analyze()

◆ fillDescriptions()

Member Data Documentation

◆ h_br

◆ h_br_ref

◆ h_decayVertexRho

◆ h_decayVertexZ

◆ h_mass

◆ h_mass_ref

◆ h_originVertexRho

◆ h_originVertexZ

◆ h_p

◆ h_plt

◆ h_plt_ref

◆ h_v

◆ knownDecayModes

◆ outputFile

◆ pids

◆ pythia

◆ tokTrack_

◆ tokVertex_