CMS 3D CMS Logo

TauSpinnerCMS.cc
Go to the documentation of this file.
2 
3 //MC-TESTER header files
4 #include "Tauola/Tauola.h"
5 #include "TauSpinner/tau_reweight_lib.h"
6 #include "TauSpinner/Tauola_wrapper.h"
8 #include "TLorentzVector.h"
9 
10 #include "CLHEP/Random/RandomEngine.h"
15 
16 using namespace edm;
17 using namespace TauSpinner;
18 
19 CLHEP::HepRandomEngine *TauSpinnerCMS::fRandomEngine = nullptr;
21 bool TauSpinnerCMS::fInitialized = false;
22 
24  : EDProducer(),
25  isReco_(pset.getParameter<bool>("isReco")),
26  isTauolaConfigured_(pset.getParameter<bool>("isTauolaConfigured")),
27  isLHPDFConfigured_(pset.getParameter<bool>("isLHPDFConfigured")),
28  LHAPDFname_(pset.getUntrackedParameter("LHAPDFname", (string)("MSTW2008nnlo90cl.LHgrid"))),
29  CMSEnergy_(pset.getParameter<double>("CMSEnergy")) //GeV
30  ,
31  gensrc_(pset.getParameter<edm::InputTag>("gensrc")),
32  MotherPDGID_(pset.getUntrackedParameter("MotherPDGID", (int)(-1))),
33  Ipol_(pset.getUntrackedParameter("Ipol", (int)(0))),
34  nonSM2_(pset.getUntrackedParameter("nonSM2", (int)(0))),
35  nonSMN_(pset.getUntrackedParameter("nonSMN", (int)(0))),
36  roundOff_(pset.getUntrackedParameter("roundOff", (double)(0.01))) {
37  //usesResource(edm::uniqueSharedResourceName());
39 
40  produces<bool>("TauSpinnerWTisValid").setBranchAlias("TauSpinnerWTisValid");
41  produces<double>("TauSpinnerWT").setBranchAlias("TauSpinnerWT");
42  produces<double>("TauSpinnerWTFlip").setBranchAlias("TauSpinnerWTFlip");
43  produces<double>("TauSpinnerWThplus").setBranchAlias("TauSpinnerWThplus");
44  produces<double>("TauSpinnerWThminus").setBranchAlias("TauSpinnerWThminus");
45 
46  if (isReco_) {
47  GenParticleCollectionToken_ = consumes<reco::GenParticleCollection>(gensrc_);
48  } else {
49  hepmcCollectionToken_ = consumes<HepMCProduct>(gensrc_);
50  }
51 }
52 
55 
57  // Now for Tauola and TauSpinner
58  Tauolapp::Tauola::setRandomGenerator(TauSpinnerCMS::flat);
59  if (!isTauolaConfigured_) {
61  }
62  if (!isLHPDFConfigured_) {
63  LHAPDF::initPDFSetByName(LHAPDFname_);
64  }
65  if (!isTauSpinnerConfigure) {
66  isTauSpinnerConfigure = true;
67  bool Ipp = true; // for pp collisions
68  // Initialize TauSpinner
69  //Ipol - polarization of input sample
70  //nonSM2 - nonstandard model calculations
71  //nonSMN
72  TauSpinner::initialize_spinner(Ipp, Ipol_, nonSM2_, nonSMN_, CMSEnergy_);
73  }
74  fInitialized = true;
75 }
76 
78  RandomEngineSentry<TauSpinnerCMS> randomEngineSentry(this, e.streamID());
79  if (!fInitialized)
80  initialize();
81 
82  Tauolapp::Tauola::setRandomGenerator(
83  TauSpinnerCMS::flat); // rest tauola++ random number incase other modules use tauola++
84  Tauolapp::jaki_.ktom = 1; // rest for when you run after tauola
85  double WT = 1.0;
86  double WTFlip = 1.0;
87  double polSM = -999; //range [-1,1]
88  SimpleParticle X, tau, tau2;
89  std::vector<SimpleParticle> tau_daughters, tau_daughters2;
90  int stat(0);
91  if (isReco_) {
92  stat = readParticlesfromReco(e, X, tau, tau2, tau_daughters, tau_daughters2);
93  } else {
95  e.getByToken(hepmcCollectionToken_, evt);
96  //Get EVENT
97  HepMC::GenEvent *Evt = new HepMC::GenEvent(*(evt->GetEvent()));
98  stat = readParticlesFromHepMC(Evt, X, tau, tau2, tau_daughters, tau_daughters2);
99  }
100  if (MotherPDGID_ < 0 || abs(X.pdgid()) == MotherPDGID_) {
101  if (stat != 1) {
102  // Determine the weight
103  if (abs(X.pdgid()) == 24 || abs(X.pdgid()) == 37) {
104  TLorentzVector tau_1r(0, 0, 0, 0);
105  TLorentzVector tau_1(tau.px(), tau.py(), tau.pz(), tau.e());
106  for (unsigned int i = 0; i < tau_daughters.size(); i++) {
107  tau_1r += TLorentzVector(
108  tau_daughters.at(i).px(), tau_daughters.at(i).py(), tau_daughters.at(i).pz(), tau_daughters.at(i).e());
109  }
110  if (fabs(tau_1r.M() - tau_1.M()) < roundOff_) {
111  WT = TauSpinner::calculateWeightFromParticlesWorHpn(
112  X, tau, tau2, tau_daughters); // note that tau2 is tau neutrino
113  polSM = getTauSpin();
114  WTFlip = (2.0 - WT) / WT;
115  }
116  } else if (X.pdgid() == 25 || X.pdgid() == 36 || X.pdgid() == 22 || X.pdgid() == 23) {
117  TLorentzVector tau_1r(0, 0, 0, 0), tau_2r(0, 0, 0, 0);
118  TLorentzVector tau_1(tau.px(), tau.py(), tau.pz(), tau.e()), tau_2(tau2.px(), tau2.py(), tau2.pz(), tau2.e());
119  for (unsigned int i = 0; i < tau_daughters.size(); i++) {
120  tau_1r += TLorentzVector(
121  tau_daughters.at(i).px(), tau_daughters.at(i).py(), tau_daughters.at(i).pz(), tau_daughters.at(i).e());
122  }
123  for (unsigned int i = 0; i < tau_daughters2.size(); i++) {
124  tau_2r += TLorentzVector(tau_daughters2.at(i).px(),
125  tau_daughters2.at(i).py(),
126  tau_daughters2.at(i).pz(),
127  tau_daughters2.at(i).e());
128  }
129 
130  if (fabs(tau_1r.M() - tau_1.M()) < roundOff_ && fabs(tau_2r.M() - tau_2.M()) < roundOff_) {
131  WT = TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters, tau_daughters2);
132  //std::cout << "WT " << WT << std::endl;
133  polSM = getTauSpin();
134  if (X.pdgid() == 25 || X.pdgid() == 22 || X.pdgid() == 23) {
135  if (X.pdgid() == 25)
136  X.setPdgid(23);
137  if (X.pdgid() == 22 || X.pdgid() == 23)
138  X.setPdgid(25);
139 
140  double WTother = TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters, tau_daughters2);
141  WTFlip = WTother / WT;
142  }
143  }
144  } else {
145  cout << "TauSpinner: WARNING: Unexpected PDG for tau mother: " << X.pdgid() << endl;
146  }
147  }
148  }
149  bool isValid = true;
150  if (!(0 <= WT && WT < 10)) {
151  isValid = false;
152  WT = 1.0;
153  WTFlip = 1.0;
154  }
155  std::unique_ptr<bool> TauSpinnerWeightisValid(new bool);
156  *TauSpinnerWeightisValid = isValid;
157  e.put(std::move(TauSpinnerWeightisValid), "TauSpinnerWTisValid");
158 
159  // regular weight
160  std::unique_ptr<double> TauSpinnerWeight(new double);
161  *TauSpinnerWeight = WT;
162  e.put(std::move(TauSpinnerWeight), "TauSpinnerWT");
163 
164  // flipped weight (ie Z->H or H->Z)
165  std::unique_ptr<double> TauSpinnerWeightFlip(new double);
166  *TauSpinnerWeightFlip = WTFlip;
167  e.put(std::move(TauSpinnerWeightFlip), "TauSpinnerWTFlip");
168 
169  // h+ polarization
170  double WThplus = WT;
171  if (polSM < 0.0 && polSM != -999 && isValid)
172  WThplus = 0;
173  std::unique_ptr<double> TauSpinnerWeighthplus(new double);
174  *TauSpinnerWeighthplus = WThplus;
175  e.put(std::move(TauSpinnerWeighthplus), "TauSpinnerWThplus");
176 
177  // h- polarization
178  double WThminus = WT;
179  if (polSM > 0.0 && polSM != -999 && isValid)
180  WThminus = 0;
181  std::unique_ptr<double> TauSpinnerWeighthminus(new double);
182  *TauSpinnerWeighthminus = WThminus;
183  e.put(std::move(TauSpinnerWeighthminus), "TauSpinnerWThminus");
184  return;
185 }
186 
188  SimpleParticle &X,
189  SimpleParticle &tau,
190  SimpleParticle &tau2,
191  std::vector<SimpleParticle> &tau_daughters,
192  std::vector<SimpleParticle> &tau2_daughters) {
195  for (reco::GenParticleCollection::const_iterator itr = genParticles->begin(); itr != genParticles->end(); ++itr) {
196  int pdgid = abs(itr->pdgId());
197  if (pdgid == 24 || pdgid == 37 || pdgid == 25 || pdgid == 36 || pdgid == 22 || pdgid == 23) {
198  const reco::GenParticle *hx = &(*itr);
199  if (!isFirst(hx))
200  continue;
201  GetLastSelf(hx);
202  const reco::GenParticle *recotau1 = nullptr;
203  const reco::GenParticle *recotau2 = nullptr;
204  unsigned int ntau(0), ntauornu(0);
205  for (unsigned int i = 0; i < itr->numberOfDaughters(); i++) {
206  const reco::Candidate *dau = itr->daughter(i);
207  if (abs(dau->pdgId()) != pdgid) {
208  if (abs(dau->pdgId()) == 15 || abs(dau->pdgId()) == 16) {
209  if (ntau == 0 && abs(dau->pdgId()) == 15) {
210  recotau1 = static_cast<const reco::GenParticle *>(dau);
211  GetLastSelf(recotau1);
212  ntau++;
213  } else if ((ntau == 1 && abs(dau->pdgId()) == 15) || abs(dau->pdgId()) == 16) {
214  recotau2 = static_cast<const reco::GenParticle *>(dau);
215  if (abs(dau->pdgId()) == 15) {
216  ntau++;
217  GetLastSelf(recotau2);
218  }
219  }
220  ntauornu++;
221  }
222  }
223  }
224  if ((ntau == 2 && ntauornu == 2) || (ntau == 1 && ntauornu == 2)) {
225  X.setPx(itr->p4().Px());
226  X.setPy(itr->p4().Py());
227  X.setPz(itr->p4().Pz());
228  X.setE(itr->p4().E());
229  X.setPdgid(itr->pdgId());
230  tau.setPx(recotau1->p4().Px());
231  tau.setPy(recotau1->p4().Py());
232  tau.setPz(recotau1->p4().Pz());
233  tau.setE(recotau1->p4().E());
234  tau.setPdgid(recotau1->pdgId());
235  GetRecoDaughters(recotau1, tau_daughters, recotau1->pdgId());
236  tau2.setPx(recotau2->p4().Px());
237  tau2.setPy(recotau2->p4().Py());
238  tau2.setPz(recotau2->p4().Pz());
239  tau2.setE(recotau2->p4().E());
240  tau2.setPdgid(recotau2->pdgId());
241  if (ntau == 2)
242  GetRecoDaughters(recotau2, tau2_daughters, recotau2->pdgId());
243  return 0;
244  }
245  }
246  }
247  return 1;
248 }
249 
251  for (unsigned int i = 0; i < Particle->numberOfDaughters(); i++) {
252  const reco::GenParticle *dau = static_cast<const reco::GenParticle *>(Particle->daughter(i));
253  if (Particle->pdgId() == dau->pdgId()) {
254  Particle = dau;
256  }
257  }
258 }
259 
261  for (unsigned int i = 0; i < Particle->numberOfMothers(); i++) {
262  const reco::GenParticle *moth = static_cast<const reco::GenParticle *>(Particle->mother(i));
263  if (Particle->pdgId() == moth->pdgId()) {
264  return false;
265  }
266  }
267  return true;
268 }
269 
271  std::vector<SimpleParticle> &daughters,
272  int parentpdgid) {
273  if (Particle->numberOfDaughters() == 0 || abs(Particle->pdgId()) == 111) {
274  SimpleParticle tp(
275  Particle->p4().Px(), Particle->p4().Py(), Particle->p4().Pz(), Particle->p4().E(), Particle->pdgId());
276  daughters.push_back(tp);
277  return;
278  }
279  for (unsigned int i = 0; i < Particle->numberOfDaughters(); i++) {
280  const reco::Candidate *dau = Particle->daughter(i);
281  GetRecoDaughters(static_cast<const reco::GenParticle *>(dau), daughters, Particle->pdgId());
282  }
283 }
284 
286  if (!fRandomEngine) {
287  throw cms::Exception("LogicError")
288  << "TauSpinnerCMS::flat: Attempt to generate random number when engine pointer is null\n"
289  << "This might mean that the code was modified to generate a random number outside the\n"
290  << "event and beginLuminosityBlock methods, which is not allowed.\n";
291  }
292  return fRandomEngine->flat();
293 }
294 
reco::Candidate::daughter
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode)
electrons_cff.bool
bool
Definition: electrons_cff.py:372
mps_fire.i
i
Definition: mps_fire.py:355
TauSpinnerCMS.h
genParticles2HepMC_cfi.genParticles
genParticles
Definition: genParticles2HepMC_cfi.py:4
TauSpinnerCMS::isTauSpinnerConfigure
static bool isTauSpinnerConfigure
Definition: TauSpinnerCMS.h:68
read_particles_from_HepMC.h
TauSpinnerCMS::endJob
void endJob() final
Definition: TauSpinnerCMS.cc:54
metsig::tau
Definition: SignAlgoResolutions.h:49
reco::GenParticle
Definition: GenParticle.h:21
X
#define X(str)
Definition: MuonsGrabber.cc:38
edm
HLT enums.
Definition: AlignableModifier.h:19
RandomNumberGenerator.h
TauSpinnerCMS::GetLastSelf
void GetLastSelf(const reco::GenParticle *Particle)
Definition: TauSpinnerCMS.cc:250
gather_cfg.cout
cout
Definition: gather_cfg.py:144
TauSpinnerCMS::roundOff_
double roundOff_
Definition: TauSpinnerCMS.h:82
TauSpinnerCMS::hepmcCollectionToken_
edm::EDGetTokenT< edm::HepMCProduct > hepmcCollectionToken_
Definition: TauSpinnerCMS.h:85
TauSpinnerCMS::isReco_
bool isReco_
Definition: TauSpinnerCMS.h:61
TauSpinnerCMS::initialize
virtual void initialize()
Definition: TauSpinnerCMS.cc:56
edm::Handle
Definition: AssociativeIterator.h:50
HepMC::GenEvent
Definition: hepmc_rootio.cc:9
TauSpinnerCMS::TauSpinnerCMS
TauSpinnerCMS(const edm::ParameterSet &)
Definition: TauSpinnerCMS.cc:23
TauSpinnerCMS::isTauolaConfigured_
bool isTauolaConfigured_
Definition: TauSpinnerCMS.h:62
TauSpinnerCMS::beginJob
void beginJob() final
Definition: TauSpinnerCMS.cc:53
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
hgcalPlots.stat
stat
Definition: hgcalPlots.py:1111
BPhysicsValidation_cfi.daughters
daughters
Definition: BPhysicsValidation_cfi.py:11
Service.h
TauSpinnerCMS::isLHPDFConfigured_
bool isLHPDFConfigured_
Definition: TauSpinnerCMS.h:63
TauSpinnerCMS::nonSMN_
int nonSMN_
Definition: TauSpinnerCMS.h:67
TauSpinnerCMS::nonSM2_
int nonSM2_
Definition: TauSpinnerCMS.h:67
TauSpinnerCMS::GenParticleCollectionToken_
edm::EDGetTokenT< reco::GenParticleCollection > GenParticleCollectionToken_
Definition: TauSpinnerCMS.h:86
TauSpinnerCMS::GetRecoDaughters
void GetRecoDaughters(const reco::GenParticle *Particle, std::vector< TauSpinner::SimpleParticle > &daughters, int parentpdgid)
Definition: TauSpinnerCMS.cc:270
cmsswSequenceInfo.tp
tp
Definition: cmsswSequenceInfo.py:17
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::ParameterSet
Definition: ParameterSet.h:36
TauSpinnerCMS::fInitialized
static bool fInitialized
Definition: TauSpinnerCMS.h:87
reco::LeafCandidate::pdgId
int pdgId() const final
PDG identifier.
Definition: LeafCandidate.h:176
createfilelist.int
int
Definition: createfilelist.py:10
TauSpinnerCMS::gensrc_
edm::InputTag gensrc_
Definition: TauSpinnerCMS.h:66
TauSpinnerCMS::fRandomEngine
static CLHEP::HepRandomEngine * fRandomEngine
Definition: TauSpinnerCMS.h:84
reco::LeafCandidate::p4
const LorentzVector & p4() const final
four-momentum Lorentz vector
Definition: LeafCandidate.h:114
edm::EventSetup
Definition: EventSetup.h:57
edm::HepMCProduct::GetEvent
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:34
itr
std::vector< std::pair< float, float > >::iterator itr
Definition: HGCDigitizer.cc:28
reco::Candidate::pdgId
virtual int pdgId() const =0
PDG identifier.
readParticlesFromHepMC
int readParticlesFromHepMC(const HepMC::GenEvent *event, TauSpinner::SimpleParticle &X, TauSpinner::SimpleParticle &tau, TauSpinner::SimpleParticle &tau2, std::vector< TauSpinner::SimpleParticle > &tau_daughters, std::vector< TauSpinner::SimpleParticle > &tau2_daughters)
reco::Candidate
Definition: Candidate.h:27
eostools.move
def move(src, dest)
Definition: eostools.py:511
RandomEngineSentry.h
edm::RandomEngineSentry
Definition: RandomEngineSentry.h:28
TauSpinnerCMS::MotherPDGID_
int MotherPDGID_
Definition: TauSpinnerCMS.h:67
Exception
Definition: hltDiff.cc:246
jets_cff.tau2
tau2
Definition: jets_cff.py:399
Particle
Definition: Particle.py:1
edm::EDProducer
Definition: EDProducer.h:36
TauSpinnerCMS::readParticlesfromReco
int readParticlesfromReco(edm::Event &e, TauSpinner::SimpleParticle &X, TauSpinner::SimpleParticle &tau, TauSpinner::SimpleParticle &tau2, std::vector< TauSpinner::SimpleParticle > &tau_daughters, std::vector< TauSpinner::SimpleParticle > &tau2_daughters)
Definition: TauSpinnerCMS.cc:187
TauSpinnerCMS::Ipol_
int Ipol_
Definition: TauSpinnerCMS.h:67
Exception.h
TauSpinnerCMS::produce
void produce(edm::Event &, const edm::EventSetup &) final
Definition: TauSpinnerCMS.cc:77
initialize
static AlgebraicMatrix initialize()
Definition: BeamSpotTransientTrackingRecHit.cc:24
TauSpinnerCMS::isFirst
bool isFirst(const reco::GenParticle *Particle)
Definition: TauSpinnerCMS.cc:260
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TauSpinnerCMS::flat
static double flat()
Definition: TauSpinnerCMS.cc:285
edm::SharedResourceNames::kTauola
static const std::string kTauola
Definition: SharedResourceNames.h:29
TauSpinnerCMS
Definition: TauSpinnerCMS.h:48
edm::Event
Definition: Event.h:73
EgammaValidation_cff.pdgid
pdgid
Definition: EgammaValidation_cff.py:30
TauSpinnerCMS::CMSEnergy_
double CMSEnergy_
Definition: TauSpinnerCMS.h:65
TauSpinnerCMS::LHAPDFname_
std::string LHAPDFname_
Definition: TauSpinnerCMS.h:64
edm::InputTag
Definition: InputTag.h:15
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37