ParticleDecay/plugins/TestPythiaDecays.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Lukas/TestPythiaDecays
0004 // Class:      TestPythiaDecays
0005 //
0006 /**\class TestPythiaDecays TestPythiaDecays.cc Lukas/TestPythiaDecays/plugins/TestPythiaDecays.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Lukas Vanelderen
0015 //         Created:  Tue, 13 May 2014 09:50:05 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/Exception.h"
0031
0032 #include "SimDataFormats/Track/interface/SimTrack.h"
0033 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0034 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0035
0036 #include "FWCore/ServiceRegistry/interface/Service.h"
0037 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0038 #include "CommonTools/Utils/interface/TFileDirectory.h"
0039
0040 // pythia
0041 #include <Pythia8/Pythia.h>
0042
0043 // root
0044 #include "TH1D.h"
0045 #include "TLorentzVector.h"
0046
0047 #if PYTHIA_VERSION_INTEGER < 8304
0048 typedef Pythia8::ParticleDataEntry* ParticleDataEntryPtr;
0049 #else
0050 typedef Pythia8::ParticleDataEntryPtr ParticleDataEntryPtr;
0051 #endif
0052 //
0053 // class declaration
0054 //
0055
0056 class TestPythiaDecays : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0057 public:
0058   explicit TestPythiaDecays(const edm::ParameterSet&);
0059   ~TestPythiaDecays() override = default;
0060
0061   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0062
0063 private:
0064   void analyze(const edm::Event&, const edm::EventSetup&) override;
0065
0066   // ----------member data ---------------------------
0067   const edm::EDGetTokenT<std::vector<SimTrack> > tokTrack_;
0068   const edm::EDGetTokenT<std::vector<SimVertex> > tokVertex_;
0069   std::vector<int> pids;
0070   std::map<int, TH1D*> h_mass;
0071   std::map<int, TH1D*> h_p;
0072   std::map<int, TH1D*> h_v;
0073   std::map<int, TH1D*> h_mass_ref;
0074   std::map<int, TH1D*> h_plt;  // plt: proper life time
0075   std::map<int, TH1D*> h_originVertexRho;
0076   std::map<int, TH1D*> h_originVertexZ;
0077   std::map<int, TH1D*> h_decayVertexRho;
0078   std::map<int, TH1D*> h_decayVertexZ;
0079   std::map<int, TH1D*> h_plt_ref;  // plt: proper life time
0080   std::map<int, TH1D*> h_br;
0081   std::map<int, TH1D*> h_br_ref;
0082
0083   std::map<int, std::vector<std::string> > knownDecayModes;
0084
0085   Pythia8::Pythia* pythia;
0086   std::string outputFile;
0087 };
0088
0089 //
0090 // constants, enums and typedefs
0091 //
0092
0093 //
0094 // static data member definitions
0095 //
0096
0097 //
0098 // constructors and destructor
0099 //
0100 TestPythiaDecays::TestPythiaDecays(const edm::ParameterSet& iConfig)
0101     : tokTrack_(consumes<std::vector<SimTrack> >(edm::InputTag("famosSimHits"))),
0102       tokVertex_(consumes<std::vector<SimVertex> >(edm::InputTag("famosSimHits"))) {
0103   usesResource(TFileService::kSharedResource);
0104   // output file
0105   outputFile = iConfig.getParameter<std::string>("outputFile");
0106
0107   // create pythia8 instance to access particle data
0108   pythia = new Pythia8::Pythia();
0109   pythia->init();
0110   Pythia8::ParticleData pdt = pythia->particleData;
0111
0112   // which particles will we study?
0113   pids.push_back(15);   // tau
0114   pids.push_back(211);  // pi+
0115   pids.push_back(111);  // pi0
0116   pids.push_back(130);  // K0L
0117   pids.push_back(321);  // K+
0118   pids.push_back(323);  // K*(392)
0119   pids.push_back(411);  // D+
0120   pids.push_back(521);  // B+
0121
0122   // define histograms
0123   edm::Service<TFileService> file;
0124   TFileDirectory observe = file->mkdir("observed");
0125   TFileDirectory predict = file->mkdir("prediction");
0126   for (size_t i = 0; i < pids.size(); ++i) {
0127     int pid = std::abs(pids[i]);
0128
0129     // get particle data
0130     if (!pdt.isParticle(pid)) {
0131       edm::LogError("TestPythiaDecays") << "ERROR: BAD PARTICLE, pythia is not aware of pid " << pid;
0132       throw cms::Exception("Unknown", "FastSim") << "Bad Particle Type " << pid << "\n";
0133       std::exit(1);
0134     }
0135     ParticleDataEntryPtr pd = pdt.particleDataEntryPtr(pid);
0136
0137     // mass histograms
0138     double m0 = pd->m0();
0139     double w = pd->mWidth();
0140     double mmin, mmax;
0141     if (w == 0) {
0142       mmin = m0 - m0 / 1000.;
0143       mmax = m0 + m0 / 1000.;
0144     } else {
0145       mmin = m0 - 2 * w;
0146       mmax = m0 + 2 * w;
0147     }
0148     std::stringstream strstr;
0149     strstr << "mass_" << pid;
0150     h_mass[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, mmin, mmax);
0151     h_mass_ref[pid] = predict.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, mmin, mmax);
0152     if (w == 0)
0153       h_mass_ref[pid]->Fill(m0);
0154     else {
0155       for (int b = 1; b <= h_mass_ref[pid]->GetNbinsX(); ++b) {
0156         double _val = h_mass_ref[pid]->GetBinCenter(b);
0157         h_mass_ref[pid]->SetBinContent(b, TMath::BreitWigner(_val, m0, w));
0158       }
0159     }
0160
0161     // p histogram
0162     strstr.str("");
0163     strstr << "p_" << pid;
0164     h_p[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 20);
0165
0166     // v histogram
0167     strstr.str("");
0168     strstr << "v_" << pid;
0169     h_v[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 1.);
0170
0171     // ctau histograms
0172     double ctau0 = pd->tau0() / 10.;
0173     strstr.str("");
0174     strstr << "plt_" << pid;
0175     h_plt[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, std::min(5. * ctau0, 500.));
0176     h_plt_ref[pid] = predict.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, std::min(5. * ctau0, 500.));
0177     h_plt_ref[pid]->SetTitle(h_plt_ref[pid]->GetName());
0178     for (int b = 1; b <= h_plt_ref[pid]->GetNbinsX(); ++b) {
0179       double _val = h_plt_ref[pid]->GetBinCenter(b);
0180       h_plt_ref[pid]->SetBinContent(b, TMath::Exp(-_val / ctau0));  //convert mm to cm
0181     }
0182
0183     // br histograms
0184     strstr.str("");
0185     strstr << "br_" << pid;
0186     h_br[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 0, 0, 0);
0187     h_br[pid]->SetCanExtend(TH1::kAllAxes);
0188     h_br_ref[pid] = predict.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 0, 0, 0);
0189     h_br_ref[pid]->SetCanExtend(TH1::kAllAxes);
0190     h_br_ref[pid]->SetTitle(h_br_ref[pid]->GetName());
0191     knownDecayModes[pid] = std::vector<std::string>();
0192     for (int d = 0; d < pd->sizeChannels(); ++d) {
0193       Pythia8::DecayChannel& channel = pd->channel(d);
0194       std::vector<int> prod;
0195       for (int p = 0; p < channel.multiplicity(); ++p) {
0196         int pId = abs(channel.product(p));
0197         // from FastSimulation/Event/src/KineParticleFilter.cc
0198         bool particleCut =
0199             (pId > 10 && pId != 12 && pId != 14 && pId != 16 && pId != 18 && pId != 21 && (pId < 23 || pId > 40) &&
0200              (pId < 81 || pId > 100) && pId != 2101 && pId != 3101 && pId != 3201 && pId != 1103 && pId != 2103 &&
0201              pId != 2203 && pId != 3103 && pId != 3203 && pId != 3303);
0202         if (particleCut)
0203           prod.push_back(abs(channel.product(p)));
0204       }
0205       std::sort(prod.begin(), prod.end());
0206       strstr.str("");
0207       for (size_t p = 0; p < prod.size(); ++p) {
0208         strstr << "_" << prod[p];
0209       }
0210       std::string label = strstr.str();
0211       h_br[pid]->Fill(label.c_str(), 0.);
0212       h_br_ref[pid]->Fill(label.c_str(), channel.bRatio());
0213       h_br[pid]->SetEntries(0);
0214       knownDecayModes[pid].push_back(label);
0215     }
0216
0217     // vertex plots
0218     strstr.str("");
0219     strstr << "originVertexRho_" << pid;
0220     h_originVertexRho[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 200);
0221     strstr.str("");
0222     strstr << "originVertexZ_" << pid;
0223     h_originVertexZ[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 400);
0224     strstr.str("");
0225     strstr << "decayVertexRho_" << pid;
0226     h_decayVertexRho[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 200);
0227     strstr.str("");
0228     strstr << "decayVertexZ_" << pid;
0229     h_decayVertexZ[pid] = observe.make<TH1D>(strstr.str().c_str(), strstr.str().c_str(), 100, 0, 400);
0230   }
0231 }
0232
0233 //
0234 // member functions
0235 //
0236
0237 // ------------ method called for each event  ------------
0238 void TestPythiaDecays::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0239   const edm::Handle<std::vector<SimTrack> >& simtracks = iEvent.getHandle(tokTrack_);
0240   const edm::Handle<std::vector<SimVertex> > simvertices = iEvent.getHandle(tokVertex_);
0241
0242   // create maps
0243
0244   // initialize
0245   std::map<size_t, std::vector<size_t> > childMap;  // child indices vs parent index
0246   std::map<size_t, int> parentMap;                  // parent index vs child index
0247   for (size_t j = 0; j < simtracks->size(); j++) {
0248     childMap[j] = std::vector<size_t>();
0249     parentMap[j] = -1;
0250   }
0251
0252   // do the mapping
0253   for (size_t j = 0; j < simtracks->size(); j++) {
0254     size_t childIndex = j;
0255     const SimTrack& child = simtracks->at(childIndex);
0256     if (child.noVertex())
0257       continue;
0258     const SimVertex& vertex = simvertices->at(child.vertIndex());
0259     if (vertex.noParent())
0260       continue;
0261     size_t parentIndex = vertex.parentIndex();
0262     childMap[parentIndex].push_back(childIndex);
0263     parentMap[childIndex] = int(parentIndex);
0264   }
0265
0266   for (size_t j = 0; j < simtracks->size(); j++) {
0267     const SimTrack& parent = simtracks->at(j);
0268     int pid = abs(parent.type());
0269     if (std::find(pids.begin(), pids.end(), pid) == pids.end())
0270       continue;
0271
0272     // fill mass hist
0273     double mass = parent.momentum().M();
0274     h_mass[pid]->Fill(mass);
0275
0276     // fill p hist
0277     h_p[pid]->Fill(parent.momentum().P());
0278
0279     // fill vertex position hist
0280     if (!parent.noVertex()) {
0281       const SimVertex& originVertex = simvertices->at(parent.vertIndex());
0282       h_originVertexRho[pid]->Fill(originVertex.position().Rho());
0283       h_originVertexZ[pid]->Fill(std::fabs(originVertex.position().Z()));
0284     }
0285     if (!childMap[j].empty()) {
0286       const SimTrack& child = simtracks->at(childMap[j][0]);
0287       const SimVertex& decayVertex = simvertices->at(child.vertIndex());
0288       h_decayVertexRho[pid]->Fill(decayVertex.position().Rho());
0289       h_decayVertexZ[pid]->Fill(std::fabs(decayVertex.position().Z()));
0290     }
0291   }
0292
0293   for (std::map<size_t, std::vector<size_t> >::iterator it = childMap.begin(); it != childMap.end(); ++it) {
0294     // fill ctau hist
0295     size_t parentIndex = it->first;
0296     const SimTrack& parent = simtracks->at(parentIndex);
0297     int pid = abs(parent.type());
0298     std::vector<size_t>& childIndices = it->second;
0299     if (childIndices.empty())
0300       continue;
0301
0302     if (std::find(pids.begin(), pids.end(), pid) == pids.end())
0303       continue;
0304
0305     const SimVertex& origin_vertex = simvertices->at(parent.vertIndex());
0306     const SimTrack& child0 = simtracks->at(childIndices[0]);
0307     const SimVertex& decay_vertex = simvertices->at(child0.vertIndex());
0308
0309     TLorentzVector lv_origin_vertex(origin_vertex.position().X(),
0310                                     origin_vertex.position().Y(),
0311                                     origin_vertex.position().Z(),
0312                                     origin_vertex.position().T());
0313     TLorentzVector lv_decay_vertex(decay_vertex.position().X(),
0314                                    decay_vertex.position().Y(),
0315                                    decay_vertex.position().Z(),
0316                                    decay_vertex.position().T());
0317     TLorentzVector lv_dist = lv_decay_vertex - lv_origin_vertex;
0318     TLorentzVector lv_parent(
0319         parent.momentum().Px(), parent.momentum().Py(), parent.momentum().Pz(), parent.momentum().E());
0320     TVector3 boost = lv_parent.BoostVector();
0321     lv_dist.Boost(-boost);
0322     h_v[pid]->Fill(boost.Mag());
0323     double plt = lv_dist.T();
0324     h_plt[pid]->Fill(plt);
0325
0326     // fill br hist
0327     std::vector<int> prod;
0328     for (size_t d = 0; d < childIndices.size(); ++d) {
0329       prod.push_back(abs(simtracks->at(childIndices[d]).type()));
0330     }
0331     std::sort(prod.begin(), prod.end());
0332     std::stringstream strstr;
0333     for (size_t p = 0; p < prod.size(); ++p) {
0334       strstr << "_" << prod[p];
0335     }
0336     std::string label = strstr.str();
0337     if (std::find(knownDecayModes[pid].begin(), knownDecayModes[pid].end(), label) == knownDecayModes[pid].end())
0338       label = "u" + label;
0339     h_br[pid]->Fill(label.c_str(), 1.);
0340     h_br_ref[pid]->Fill(label.c_str(), 0.);  // keep h_br and h_br_ref in sync
0341   }
0342 }
0343
0344 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0345 void TestPythiaDecays::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0346   //The following says we do not know what parameters are allowed so do no validation
0347   // Please change this to state exactly what you do use, even if it is no parameters
0348   edm::ParameterSetDescription desc;
0349   desc.setUnknown();
0350   descriptions.addDefault(desc);
0351 }
0352
0353 //define this as a plug-in
0354 DEFINE_FWK_MODULE(TestPythiaDecays);