TestPythiaDecays Class Reference

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

Inheritance diagram for TestPythiaDecays:

Public Member Functions
	TestPythiaDecays (const edm::ParameterSet &)
	~TestPythiaDecays ()
Public Member Functions inherited from edm::EDAnalyzer
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
	EDAnalyzer ()
ModuleDescription const &	moduleDescription () const
std::string	workerType () const
virtual	~EDAnalyzer ()
Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()
ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const
void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const
void	labelsForToken (EDGetToken iToken, Labels &oLabels) const
void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const
bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const
bool	registeredToConsumeMany (TypeID const &, BranchType) const
void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
virtual	~EDConsumerBase ()

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)
Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()
static void	fillDescriptions (ConfigurationDescriptions &descriptions)
static void	prevalidate (ConfigurationDescriptions &)

Private Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &) override
virtual void	beginJob () override
virtual void	endJob () 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 < string > >	knownDecayModes
std::string	outputFile
std::vector< int >	pids
Pythia8::Pythia *	pythia

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType
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 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

Description: [one line class summary]

Implementation: [Notes on implementation]

Definition at line 48 of file TestPythiaDecays.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 98 of file TestPythiaDecays.cc.

References funct::abs(), b, gather_cfg::cout, cmsRelvalreport::exit, 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, i, knownDecayModes, diffTwoXMLs::label, min(), outputFile, AlCaHLTBitMon_ParallelJobs::p, sysUtil::pid, pids, parseEventContent::prod, pythia, python.multivaluedict::sort(), AlCaHLTBitMon_QueryRunRegistry::string, and w.

{

  // 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
  for(size_t i = 0;i<pids.size();++i){

    int pid = abs(pids[i]);

    // get particle data
    if(!pdt.isParticle(pid)){
      std::cout << "ERROR: BAD PARTICLE, pythia is not aware of pid " << pid << std::endl;
      std::exit(1);
    }
    Pythia8::ParticleDataEntry * 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] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),100,mmin,mmax);    
    h_mass_ref[pid] = (TH1D*)(h_mass[pid]->Clone(strstr.str().c_str()));
    h_mass_ref[pid]->SetTitle(h_mass_ref[pid]->GetName());
    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] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),100,0,20);    

    // v histogram
    strstr.str("");
    strstr << "v_" << pid;
    h_v[pid] = new 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] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),100,0,std::min(5.*ctau0,500.));
    h_plt_ref[pid] = (TH1D*)(h_plt[pid]->Clone(strstr.str().c_str()));
    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] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),0,0,0);
    h_br[pid]->SetBit(TH1::kCanRebin);
    h_br_ref[pid] = (TH1D*)(h_br[pid]->Clone(strstr.str().c_str()));
    h_br_ref[pid]->SetTitle(h_br_ref[pid]->GetName());
    knownDecayModes[pid] = vector<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] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),100,0,200);
    strstr.str("");
    strstr << "originVertexZ_" << pid;
    h_originVertexZ[pid] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),100,0,400);
    strstr.str("");
    strstr << "decayVertexRho_" << pid;
    h_decayVertexRho[pid] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),100,0,200);
    strstr.str("");
    strstr << "decayVertexZ_" << pid;
    h_decayVertexZ[pid] = new TH1D(strstr.str().c_str(),strstr.str().c_str(),100,0,400);
  }
  
}

TestPythiaDecays::~TestPythiaDecays

(

)

Definition at line 234 of file TestPythiaDecays.cc.

References f, 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, i, outputFile, sysUtil::pid, and pids.

{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)
  TFile * f = TFile::Open(outputFile.c_str(),"RECREATE");
  f->cd();
  f->mkdir("observed");
  f->mkdir("prediction");
  for(size_t i = 0;i<pids.size();++i){
    int pid = pids[i];
    f->cd("observed");
    h_mass[pid]->Write();
    h_plt[pid]->Write();
    h_br[pid]->Write();
    h_originVertexZ[pid]->Write();
    h_originVertexRho[pid]->Write();
    h_decayVertexZ[pid]->Write();
    h_decayVertexRho[pid]->Write();
    h_p[pid]->Write();
    h_v[pid]->Write();
    f->cd("prediction");
    h_mass_ref[pid]->Write();
    h_plt_ref[pid]->Write();
    h_br_ref[pid]->Write();
  }
  f->Close();
}

Member Function Documentation

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

overrideprivatevirtual

Implements edm::EDAnalyzer.

Definition at line 270 of file TestPythiaDecays.cc.

References funct::abs(), begin, end, spr::find(), edm::Event::getByLabel(), h_br, h_br_ref, h_decayVertexRho, h_decayVertexZ, h_mass, h_originVertexRho, h_originVertexZ, h_p, h_plt, h_v, j, knownDecayModes, diffTwoXMLs::label, python.multivaluedict::map(), CoreSimTrack::momentum(), SimVertex::noParent(), SimTrack::noVertex(), AlCaHLTBitMon_ParallelJobs::p, dbtoconf::parent, SimVertex::parentIndex(), sysUtil::pid, pids, CoreSimVertex::position(), parseEventContent::prod, findQualityFiles::size, python.multivaluedict::sort(), AlCaHLTBitMon_QueryRunRegistry::string, CoreSimTrack::type(), and SimTrack::vertIndex().

{
   using namespace edm;

   Handle<std::vector<SimTrack> > simtracks;
   iEvent.getByLabel("famosSimHits",simtracks);

   Handle<std::vector<SimVertex> > simvertices;
   iEvent.getByLabel("famosSimHits",simvertices);

   // 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].size() > 0){
       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());
     vector<size_t> & childIndices = it->second;
     if(childIndices.size() == 0)
       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
   }
}

void TestPythiaDecays::beginJob ( void  )

overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 382 of file TestPythiaDecays.cc.

{
}

void TestPythiaDecays::endJob ( void  )

overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 388 of file TestPythiaDecays.cc.

{
}

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

static

Definition at line 426 of file TestPythiaDecays.cc.

References edm::ConfigurationDescriptions::addDefault(), and edm::ParameterSetDescription::setUnknown().

                                                                             {
  //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

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

private

Definition at line 78 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 79 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 75 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 76 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 68 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 71 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 73 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 74 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 69 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 72 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 77 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 70 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

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

private

Definition at line 81 of file TestPythiaDecays.cc.

Referenced by analyze(), and TestPythiaDecays().

std::string TestPythiaDecays::outputFile

private

Definition at line 84 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays(), and ~TestPythiaDecays().

std::vector<int> TestPythiaDecays::pids

private

Definition at line 67 of file TestPythiaDecays.cc.

Referenced by analyze(), TestPythiaDecays(), and ~TestPythiaDecays().

Pythia8::Pythia* TestPythiaDecays::pythia

private

Definition at line 83 of file TestPythiaDecays.cc.

Referenced by TestPythiaDecays().