test
CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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;
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  ,gensrc_(pset.getParameter<edm::InputTag>("gensrc"))
31  ,MotherPDGID_(pset.getUntrackedParameter("MotherPDGID",(int)(-1)))
32  ,Ipol_(pset.getUntrackedParameter("Ipol",(int)(0)))
33  ,nonSM2_(pset.getUntrackedParameter("nonSM2",(int)(0)))
34  ,nonSMN_(pset.getUntrackedParameter("nonSMN",(int)(0)))
35  ,roundOff_(pset.getUntrackedParameter("roundOff",(double)(0.01)))
36 {
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  }
49  else{
50  hepmcCollectionToken_=consumes<HepMCProduct>(gensrc_);
51  }
52 }
53 
56 
58  // Now for Tauola and TauSpinner
59  Tauolapp::Tauola::setRandomGenerator(TauSpinnerCMS::flat);
62  }
63  if(!isLHPDFConfigured_){
64  LHAPDF::initPDFSetByName(LHAPDFname_);
65  }
68  bool Ipp = true; // for pp collisions
69  // Initialize TauSpinner
70  //Ipol - polarization of input sample
71  //nonSM2 - nonstandard model calculations
72  //nonSMN
73  TauSpinner::initialize_spinner(Ipp,Ipol_,nonSM2_,nonSMN_,CMSEnergy_);
74  }
75  fInitialized=true;
76 }
77 
78 
80  RandomEngineSentry<TauSpinnerCMS> randomEngineSentry(this, e.streamID());
81  if(!fInitialized) initialize();
82 
83  Tauolapp::Tauola::setRandomGenerator(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  }
94  else{
97  //Get EVENT
98  HepMC::GenEvent *Evt = new HepMC::GenEvent(*(evt->GetEvent()));
99  stat=readParticlesFromHepMC(Evt,X,tau,tau2,tau_daughters,tau_daughters2);
100  }
101  if(MotherPDGID_<0 || abs(X.pdgid())==MotherPDGID_){
102  if(stat!=1){
103  // Determine the weight
104  if( abs(X.pdgid())==24 || abs(X.pdgid())==37 ){
105  TLorentzVector tau_1r(0,0,0,0);
106  TLorentzVector tau_1(tau.px(),tau.py(),tau.pz(),tau.e());
107  for(unsigned int i=0; i<tau_daughters.size();i++){
108  tau_1r+=TLorentzVector(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(X, tau, tau2, tau_daughters); // note that tau2 is tau neutrino
112  polSM=getTauSpin();
113  WTFlip=(2.0-WT)/WT;
114  }
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(tau_daughters.at(i).px(),tau_daughters.at(i).py(),tau_daughters.at(i).pz(),tau_daughters.at(i).e());
121  }
122  for(unsigned int i=0; i<tau_daughters2.size();i++){
123  tau_2r+=TLorentzVector(tau_daughters2.at(i).px(),tau_daughters2.at(i).py(),tau_daughters2.at(i).pz(),tau_daughters2.at(i).e());
124  }
125 
126  if(fabs(tau_1r.M()-tau_1.M())<roundOff_ && fabs(tau_2r.M()-tau_2.M())<roundOff_){
127  WT = TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters,tau_daughters2);
128  //std::cout << "WT " << WT << std::endl;
129  polSM=getTauSpin();
130  if(X.pdgid()==25 || X.pdgid()==22 || X.pdgid()==23 ){
131  if(X.pdgid()==25) X.setPdgid(23);
132  if( X.pdgid()==22 || X.pdgid()==23 ) X.setPdgid(25);
133 
134  double WTother=TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters,tau_daughters2);
135  WTFlip=WTother/WT;
136  }
137  }
138  }
139  else{
140  cout<<"TauSpinner: WARNING: Unexpected PDG for tau mother: "<<X.pdgid()<<endl;
141  }
142  }
143  }
144  bool isValid=true;
145  if(!(0<=WT && WT<10)){isValid=false; WT=1.0; WTFlip=1.0;}
146  std::auto_ptr<bool> TauSpinnerWeightisValid(new bool);
147  *TauSpinnerWeightisValid =isValid;
148  e.put(TauSpinnerWeightisValid,"TauSpinnerWTisValid");
149 
150  // regular weight
151  std::auto_ptr<double> TauSpinnerWeight(new double);
152  *TauSpinnerWeight =WT;
153  e.put(TauSpinnerWeight,"TauSpinnerWT");
154 
155  // flipped weight (ie Z->H or H->Z)
156  std::auto_ptr<double> TauSpinnerWeightFlip(new double);
157  *TauSpinnerWeightFlip =WTFlip;
158  e.put(TauSpinnerWeightFlip,"TauSpinnerWTFlip");
159 
160  // h+ polarization
161  double WThplus=WT;
162  if(polSM<0.0 && polSM!=-999 && isValid) WThplus=0;
163  std::auto_ptr<double> TauSpinnerWeighthplus(new double);
164  *TauSpinnerWeighthplus = WThplus;
165  e.put(TauSpinnerWeighthplus,"TauSpinnerWThplus");
166 
167  // h- polarization
168  double WThminus=WT;
169  if(polSM>0.0&& polSM!=-999 && isValid) WThminus=0;
170  std::auto_ptr<double> TauSpinnerWeighthminus(new double);
171  *TauSpinnerWeighthminus = WThminus;
172  e.put(TauSpinnerWeighthminus,"TauSpinnerWThminus");
173  return ;
174 }
175 
176 int TauSpinnerCMS::readParticlesfromReco(edm::Event& e,SimpleParticle &X,SimpleParticle &tau,SimpleParticle &tau2,
177  std::vector<SimpleParticle> &tau_daughters,std::vector<SimpleParticle> &tau2_daughters){
179  e.getByToken(GenParticleCollectionToken_, genParticles);
180  for(reco::GenParticleCollection::const_iterator itr = genParticles->begin(); itr!= genParticles->end(); ++itr){
181  int pdgid=abs(itr->pdgId());
182  if(pdgid==24 || pdgid==37 || pdgid ==25 || pdgid==36 || pdgid==22 || pdgid==23 ){
183  const reco::GenParticle *hx=&(*itr);
184  if(!isFirst(hx)) continue;
185  GetLastSelf(hx);
186  const reco::GenParticle *recotau1=NULL;
187  const reco::GenParticle *recotau2=NULL;
188  unsigned int ntau(0),ntauornu(0);
189  for(unsigned int i=0; i<itr->numberOfDaughters(); i++){
190  const reco::Candidate *dau=itr->daughter(i);
191  if(abs(dau->pdgId())!=pdgid){
192  if(abs(dau->pdgId())==15 || abs(dau->pdgId())==16){
193  if(ntau==0 && abs(dau->pdgId())==15){
194  recotau1=static_cast<const reco::GenParticle*>(dau);
195  GetLastSelf(recotau1);
196  ntau++;
197  }
198  else if((ntau==1 && abs(dau->pdgId())==15) || abs(dau->pdgId())==16){
199  recotau2=static_cast<const reco::GenParticle*>(dau);
200  if(abs(dau->pdgId())==15){ntau++;GetLastSelf(recotau2);}
201  }
202  ntauornu++;
203  }
204  }
205  }
206  if((ntau==2 && ntauornu==2) || (ntau==1 && ntauornu==2)){
207  X.setPx(itr->p4().Px());
208  X.setPy(itr->p4().Py());
209  X.setPz(itr->p4().Pz());
210  X.setE (itr->p4().E());
211  X.setPdgid(itr->pdgId());
212  tau.setPx(recotau1->p4().Px());
213  tau.setPy(recotau1->p4().Py());
214  tau.setPz(recotau1->p4().Pz());
215  tau.setE (recotau1->p4().E());
216  tau.setPdgid(recotau1->pdgId());
217  GetRecoDaughters(recotau1,tau_daughters,recotau1->pdgId());
218  tau2.setPx(recotau2->p4().Px());
219  tau2.setPy(recotau2->p4().Py());
220  tau2.setPz(recotau2->p4().Pz());
221  tau2.setE (recotau2->p4().E());
222  tau2.setPdgid(recotau2->pdgId());
223  if(ntau==2)GetRecoDaughters(recotau2,tau2_daughters,recotau2->pdgId());
224  return 0;
225  }
226  }
227  }
228  return 1;
229 }
230 
232  for (unsigned int i=0; i< Particle->numberOfDaughters(); i++){
233  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(Particle->daughter(i));
234  if(Particle->pdgId()==dau->pdgId()){
235  Particle=dau;
236  GetLastSelf(Particle);
237  }
238  }
239 }
240 
242  for (unsigned int i=0; i< Particle->numberOfMothers(); i++){
243  const reco::GenParticle *moth=static_cast<const reco::GenParticle*>(Particle->mother(i));
244  if(Particle->pdgId()==moth->pdgId()){
245  return false;
246  }
247  }
248  return true;
249 }
250 
251 void TauSpinnerCMS::GetRecoDaughters(const reco::GenParticle *Particle,std::vector<SimpleParticle> &daughters, int parentpdgid){
252  if( Particle->numberOfDaughters()==0 || abs(Particle->pdgId())==111){
253  SimpleParticle tp(Particle->p4().Px(), Particle->p4().Py(), Particle->p4().Pz(), Particle->p4().E(), Particle->pdgId());
254  daughters.push_back(tp);
255  return;
256  }
257  for (unsigned int i=0; i< Particle->numberOfDaughters(); i++){
258  const reco::Candidate *dau=Particle->daughter(i);
259  GetRecoDaughters(static_cast<const reco::GenParticle*>(dau),daughters,Particle->pdgId());
260  }
261 }
262 
264 {
265  if ( !fRandomEngine ) {
266  throw cms::Exception("LogicError")
267  << "TauSpinnerCMS::flat: Attempt to generate random number when engine pointer is null\n"
268  << "This might mean that the code was modified to generate a random number outside the\n"
269  << "event and beginLuminosityBlock methods, which is not allowed.\n";
270  }
271  return fRandomEngine->flat();
272 }
273 
static AlgebraicMatrix initialize()
int i
Definition: DBlmapReader.cc:9
void GetLastSelf(const reco::GenParticle *Particle)
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
edm::EDGetTokenT< edm::HepMCProduct > hepmcCollectionToken_
Definition: TauSpinnerCMS.h:92
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual void initialize()
static double flat()
bool isFirst(const reco::GenParticle *Particle)
#define X(str)
Definition: MuonsGrabber.cc:48
virtual void produce(edm::Event &, const edm::EventSetup &) overridefinal
#define NULL
Definition: scimark2.h:8
static const std::string kTauola
double roundOff_
Definition: TauSpinnerCMS.h:89
static CLHEP::HepRandomEngine * fRandomEngine
Definition: TauSpinnerCMS.h:91
edm::EDGetTokenT< reco::GenParticleCollection > GenParticleCollectionToken_
Definition: TauSpinnerCMS.h:93
void GetRecoDaughters(const reco::GenParticle *Particle, std::vector< TauSpinner::SimpleParticle > &daughters, int parentpdgid)
static bool isTauSpinnerConfigure
Definition: TauSpinnerCMS.h:75
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:121
virtual void beginJob() overridefinal
virtual size_t numberOfMothers() const
number of mothers
virtual size_t numberOfDaughters() const
number of daughters
virtual const Candidate * daughter(size_type) const
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
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)
virtual int pdgId() const =0
PDG identifier.
bool isTauolaConfigured_
Definition: TauSpinnerCMS.h:69
std::string LHAPDFname_
Definition: TauSpinnerCMS.h:71
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)
static bool fInitialized
Definition: TauSpinnerCMS.h:94
edm::InputTag gensrc_
Definition: TauSpinnerCMS.h:73
return(e1-e2)*(e1-e2)+dp *dp
virtual int pdgId() const final
PDG identifier.
TauSpinnerCMS(const edm::ParameterSet &)
StreamID streamID() const
Definition: Event.h:80
virtual void endJob() overridefinal
double CMSEnergy_
Definition: TauSpinnerCMS.h:72
tuple cout
Definition: gather_cfg.py:145
bool isLHPDFConfigured_
Definition: TauSpinnerCMS.h:70
virtual const Candidate * mother(size_type=0) const
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...
virtual const LorentzVector & p4() const final
four-momentum Lorentz vector
Definition: LeafCandidate.h:99