Tauolappand tau spinner

GeneratorInterface/ExternalDecays/BuildFile.xml

@@ -1,11 +1,12 @@
-<use   name="FWCore/ParameterSet"/>
-<use   name="GeneratorInterface/Pythia6Interface"/>
-<use   name="heppdt"/>
-<use   name="SimGeneral/HepPDTRecord"/>
-<use   name="tauolapp"/>
-<use   name="photos"/>
-<use   name="clhep"/>
-<use   name="evtgenlhc"/>
+<use name="FWCore/ParameterSet"/>
+<use name="GeneratorInterface/Pythia6Interface"/>
+<use name="heppdt"/>
+<use name="SimGeneral/HepPDTRecord"/>
+<use name="tauolapp"/>
+<use name="photos"/>
+<use name="clhep"/>
+<use name="evtgenlhc"/>
+<use name="root"/>
 <export>
   <lib   name="1"/>
 </export>

GeneratorInterface/ExternalDecays/interface/TauSpinnerCMS.h

@@ -0,0 +1,72 @@
+// -*- C++ -*-
+//
+// Package:    TauSpinnerInterface
+// Class:      TauSpinnerCMS
+// 
+/**\class TauSpinnerCMS TauSpinnerCMS.cc 
+
+*/
+//
+// Original Author:  Ian Nugent  
+//         Created:  Fri Feb 15 2013
+
+#ifndef TauSpinnerCMS_h
+#define TauSpinnerCMS_h
+
+
+#include <iostream>
+
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+
+// essentials !!!
+#include "FWCore/Framework/interface/Event.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ServiceRegistry/interface/Service.h" 
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+#include "TH1.h"
+
+#include "FWCore/Framework/interface/EDProducer.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
+#include "TauSpinner/SimpleParticle.h"
+
+using namespace edm;
+using namespace std;
+using namespace TauSpinner;
+
+class TauSpinnerCMS : public edm::EDProducer
+{
+
+ public:
+
+  //
+  explicit TauSpinnerCMS( const edm::ParameterSet& ) ;
+  virtual ~TauSpinnerCMS() {} // no need to delete ROOT stuff
+  // as it'll be deleted upon closing TFile
+
+  virtual void produce( edm::Event&, const edm::EventSetup&);
+  virtual void beginJob() ;
+  virtual void endRun( const edm::Run&, const edm::EventSetup& ) ;
+  virtual void endJob() ;
+
+ private:
+  bool isReco_;
+  bool isTauolaConfigured_;
+  bool isLHPDFConfigured_;
+  string LHAPDFname_;
+  double CMSEnergy_;
+  edm::InputTag gensrc_;
+  int MotherPDGID_,Ipol_,nonSM2_,nonSMN_;
+  static bool isTauSpinnerConfigure;
+
+  // Additional funtionms for Reco (not provided by Tauola/TauSpinner authors)
+  int readParticlesfromReco(edm::Event& e,SimpleParticle &X,SimpleParticle &tau,SimpleParticle &tau2,
+			    std::vector<SimpleParticle> &tau_daughters,std::vector<SimpleParticle> &tau2_daughters);
+  void GetLastSelf(const reco::GenParticle *Particle);
+  void GetRecoDaughters(const reco::GenParticle *Particle,std::vector<SimpleParticle> &daughters,int parentpdgid);
+  bool isFirst(const reco::GenParticle *Particle);
+  double roundOff_;
+}; 
+#endif

GeneratorInterface/ExternalDecays/interface/read_particles_from_HepMC.h

@@ -0,0 +1,45 @@
+#ifndef _READ_PARTICLES_FROM_HEPMC_H_
+#define _READ_PARTICLES_FROM_HEPMC_H_
+
+// HepMC headers
+#include "HepMC/GenEvent.h"
+#include "HepMC/GenParticle.h"
+#include "HepMC/IO_GenEvent.h"
+
+#include "TauSpinner/SimpleParticle.h"
+#include <vector>
+using namespace std;
+using namespace TauSpinner;
+
+/** Read HepMC::GenEvent.
+
+  Read HepMC event from data file
+  and return particles needed for tau spin weight calculation.
+
+  This routine is prepared for use with files generated by Pythia8.
+  Fills:
+
+  'X'              - Heavy particle (W+/-, H+/-, H, Z)
+  'tau'            - first tau
+
+  'tau_daughters'  - daughters of 'tau'
+  'tau2_daughters' - daughters of 'tau2' or empty list, if 'tau2' is nu_tau.
+
+  Returns:
+  0 - event processed correctly 
+  1 - empty event or no decay found in the event
+*/
+int readParticlesFromHepMC(const HepMC::GenEvent *event, SimpleParticle &X, SimpleParticle &tau, SimpleParticle &tau2, vector<SimpleParticle> &tau_daughters, vector<SimpleParticle> &tau2_daughters);
+
+/** Get daughters of HepMC::GenParticle
+
+  Recursively searches for final-state daughters of 'x' */
+vector<SimpleParticle> *getDaughters(HepMC::GenParticle *x);
+
+/** Find last self
+
+  Recursively finds the last particle with the same PDG ID
+  on the list of its decay products */
+HepMC::GenParticle *findLastSelf(HepMC::GenParticle *x);
+bool isFirst(HepMC::GenParticle *x);
+#endif

GeneratorInterface/ExternalDecays/plugins/BuildFile.xml

@@ -0,0 +1,9 @@
+<use name="FWCore/Framework"/>
+<use name="FWCore/ParameterSet"/>
+<use name="lhapdf"/>
+<use name="SimDataFormats/GeneratorProducts"/> 
+<use name="tauolapp"/>
+<use name="GeneratorInterface/ExternalDecays"/>
+<library   file="*.cc" name="GeneratorInterfaceExternalDecaysPlugins">
+  <flags   EDM_PLUGIN="1"/>
+</library>

GeneratorInterface/ExternalDecays/plugins/TauSpinnerCMS.cc

@@ -0,0 +1,234 @@
+#include "GeneratorInterface/ExternalDecays/interface/TauSpinnerCMS.h"
+
+//MC-TESTER header files
+#include "Tauola/Tauola.h"
+#include "TauSpinner/tau_reweight_lib.h"
+#include "TauSpinner/Tauola_wrapper.h"
+#include "GeneratorInterface/ExternalDecays/interface/read_particles_from_HepMC.h"
+#include "TLorentzVector.h"
+
+bool TauSpinnerCMS::isTauSpinnerConfigure=false;
+
+TauSpinnerCMS::TauSpinnerCMS( const ParameterSet& pset ) :
+  isReco_(pset.getParameter<bool>("isReco"))
+  ,isTauolaConfigured_(pset.getParameter<bool>("isTauolaConfigured" ))
+  ,isLHPDFConfigured_(pset.getParameter<bool>("isLHPDFConfigured" ))
+  ,LHAPDFname_(pset.getUntrackedParameter("LHAPDFname",(string)("MSTW2008nnlo90cl.LHgrid")))
+  ,CMSEnergy_(pset.getParameter<double>("CMSEnergy"))//GeV
+  ,gensrc_(pset.getParameter<edm::InputTag>("gensrc"))
+  ,MotherPDGID_(pset.getUntrackedParameter("MotherPDGID",(int)(-1)))
+  ,Ipol_(pset.getUntrackedParameter("Ipol",(int)(0)))
+  ,nonSM2_(pset.getUntrackedParameter("nonSM2",(int)(0)))
+  ,nonSMN_(pset.getUntrackedParameter("nonSMN",(int)(0)))
+  ,roundOff_(pset.getUntrackedParameter("roundOff",(double)(0.01)))
+{
+  produces<bool>("TauSpinnerWTisValid").setBranchAlias("TauSpinnerWTisValid");
+  produces<double>("TauSpinnerWT").setBranchAlias("TauSpinnerWT");
+  produces<double>("TauSpinnerWTFlip").setBranchAlias("TauSpinnerWTFlip");
+  produces<double>("TauSpinnerWThplus").setBranchAlias("TauSpinnerWThplus");
+  produces<double>("TauSpinnerWThminus").setBranchAlias("TauSpinnerWThminus");
+}
+
+void TauSpinnerCMS::beginJob()
+{
+  if(!isTauolaConfigured_){
+    Tauolapp::Tauola::initialize();
+  }
+  if(!isLHPDFConfigured_){
+    LHAPDF::initPDFSetByName(LHAPDFname_);   
+  }
+  if(!isTauSpinnerConfigure){
+    isTauSpinnerConfigure=true;
+    bool Ipp = true;  // for pp collisions 
+    // Initialize TauSpinner
+    //Ipol - polarization of input sample
+    //nonSM2 - nonstandard model calculations
+    //nonSMN
+    TauSpinner::initialize_spinner(Ipp,Ipol_,nonSM2_,nonSMN_,CMSEnergy_);
+  }
+}
+
+void TauSpinnerCMS::produce( edm::Event& e, const edm::EventSetup& iSetup){
+  double WT=1.0;
+  double WTFlip=1.0;
+  double polSM=-999; //range [-1,1]
+  SimpleParticle X, tau, tau2;
+  vector<SimpleParticle> tau_daughters, tau_daughters2;
+  int stat(0);
+  if(isReco_){
+    stat=readParticlesfromReco(e,X,tau,tau2,tau_daughters,tau_daughters2);
+  }
+  else{
+    Handle< HepMCProduct > EvtHandle ;
+    e.getByLabel( "generator", EvtHandle ) ;
+    const HepMC::GenEvent* Evt = EvtHandle->GetEvent() ;
+    stat=readParticlesFromHepMC(Evt,X,tau,tau2,tau_daughters,tau_daughters2);
+  }  
+  if(MotherPDGID_<0 || abs(X.pdgid())==MotherPDGID_){
+    if(stat!=1){
+      // Determine the weight      
+      if( abs(X.pdgid())==24 ||  abs(X.pdgid())==37 ){
+        TLorentzVector tau_1r(0,0,0,0);
+        TLorentzVector tau_1(tau.px(),tau.py(),tau.pz(),tau.e());
+        for(unsigned int i=0; i<tau_daughters.size();i++){
+          tau_1r+=TLorentzVector(tau_daughters.at(i).px(),tau_daughters.at(i).py(),tau_daughters.at(i).pz(),tau_daughters.at(i).e());
+        }
+	if(fabs(tau_1r.M()-tau_1.M())<roundOff_){
+	  WT = TauSpinner::calculateWeightFromParticlesWorHpn(X, tau, tau2, tau_daughters); // note that tau2 is tau neutrino
+	  polSM=getTauSpin();
+	  WTFlip=(2.0-WT)/WT;
+	}
+      }
+      else if( X.pdgid()==25 || X.pdgid()==36 || X.pdgid()==22 || X.pdgid()==23 ){
+	TLorentzVector tau_1r(0,0,0,0), tau_2r(0,0,0,0);
+	TLorentzVector tau_1(tau.px(),tau.py(),tau.pz(),tau.e()), tau_2(tau2.px(),tau2.py(),tau2.pz(),tau2.e());
+	for(unsigned int i=0; i<tau_daughters.size();i++){
+	  tau_1r+=TLorentzVector(tau_daughters.at(i).px(),tau_daughters.at(i).py(),tau_daughters.at(i).pz(),tau_daughters.at(i).e());
+	}
+	for(unsigned int i=0; i<tau_daughters2.size();i++){
+	  tau_2r+=TLorentzVector(tau_daughters2.at(i).px(),tau_daughters2.at(i).py(),tau_daughters2.at(i).pz(),tau_daughters2.at(i).e());
+	}
+
+        if(fabs(tau_1r.M()-tau_1.M())<roundOff_ && fabs(tau_2r.M()-tau_2.M())<roundOff_){
+          WT = TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters,tau_daughters2);
+	  //std::cout << "WT " << WT << std::endl;
+          polSM=getTauSpin();
+          if(X.pdgid()==25 || X.pdgid()==22 || X.pdgid()==23 ){
+            if(X.pdgid()==25) X.setPdgid(23);
+            if( X.pdgid()==22 || X.pdgid()==23 ) X.setPdgid(25);
+
+            double WTother=TauSpinner::calculateWeightFromParticlesH(X, tau, tau2, tau_daughters,tau_daughters2);
+            WTFlip=WTother/WT;
+          }
+        }
+      }
+      else{
+	cout<<"TauSpinner: WARNING: Unexpected PDG for tau mother: "<<X.pdgid()<<endl;
+      }
+    }
+  }
+  bool isValid=true;
+  if(!(0<=WT && WT<10)){isValid=false; WT=1.0; WTFlip=1.0;}
+  std::auto_ptr<bool> TauSpinnerWeightisValid(new bool);
+  *TauSpinnerWeightisValid =isValid;
+  e.put(TauSpinnerWeightisValid,"TauSpinnerWTisValid");
+
+  // regular weight
+  std::auto_ptr<double> TauSpinnerWeight(new double);
+  *TauSpinnerWeight =WT;    
+  e.put(TauSpinnerWeight,"TauSpinnerWT");  
+
+  // flipped weight (ie Z->H or H->Z)
+  std::auto_ptr<double> TauSpinnerWeightFlip(new double);
+  *TauSpinnerWeightFlip =WTFlip;
+  e.put(TauSpinnerWeightFlip,"TauSpinnerWTFlip");
+
+  // h+ polarization
+  double WThplus=WT;
+  if(polSM<0.0 && polSM!=-999 && isValid) WThplus=0; 
+  std::auto_ptr<double> TauSpinnerWeighthplus(new double);
+  *TauSpinnerWeighthplus = WThplus;
+  e.put(TauSpinnerWeighthplus,"TauSpinnerWThplus");
+
+  // h- polarization
+  double WThminus=WT;
+  if(polSM>0.0&& polSM!=-999 && isValid) WThminus=0;
+  std::auto_ptr<double> TauSpinnerWeighthminus(new double);
+  *TauSpinnerWeighthminus = WThminus;
+  e.put(TauSpinnerWeighthminus,"TauSpinnerWThminus");
+
+  return ;
+}  
+
+void TauSpinnerCMS::endRun( const edm::Run& r, const edm::EventSetup& ){}
+
+void TauSpinnerCMS::endJob(){}
+
+int TauSpinnerCMS::readParticlesfromReco(edm::Event& e,SimpleParticle &X,SimpleParticle &tau,SimpleParticle &tau2, 
+					 std::vector<SimpleParticle> &tau_daughters,std::vector<SimpleParticle> &tau2_daughters){
+  edm::Handle<reco::GenParticleCollection> genParticles;
+  e.getByLabel(gensrc_, genParticles);
+  for(reco::GenParticleCollection::const_iterator itr = genParticles->begin(); itr!= genParticles->end(); ++itr){
+    int pdgid=abs(itr->pdgId());
+    if(pdgid==24 || pdgid==37 || pdgid ==25 || pdgid==36 || pdgid==22 || pdgid==23 ){
+      const reco::GenParticle *hx=&(*itr);
+      if(!isFirst(hx)) continue;
+      GetLastSelf(hx);
+      const reco::GenParticle *recotau1=NULL;
+      const reco::GenParticle *recotau2=NULL;
+      unsigned int ntau(0),ntauornu(0);
+      for(unsigned int i=0; i<itr->numberOfDaughters(); i++){
+	const reco::Candidate *dau=itr->daughter(i);
+	if(abs(dau->pdgId())!=pdgid){
+	  if(abs(dau->pdgId())==15 || abs(dau->pdgId())==16){
+	    if(ntau==0 && abs(dau->pdgId())==15){
+	      recotau1=static_cast<const reco::GenParticle*>(dau);
+	      GetLastSelf(recotau1);
+	      ntau++;
+	    }
+	    else if((ntau==1 && abs(dau->pdgId())==15) || abs(dau->pdgId())==16){
+	      recotau2=static_cast<const reco::GenParticle*>(dau);
+	      if(abs(dau->pdgId())==15){ntau++;GetLastSelf(recotau2);}
+	    }
+	    ntauornu++;
+	  }
+	}
+      }
+      if((ntau==2 && ntauornu==2) || (ntau==1 && ntauornu==2)){
+	X.setPx(itr->p4().Px());
+	X.setPy(itr->p4().Py());
+	X.setPz(itr->p4().Pz());
+	X.setE (itr->p4().E());
+	X.setPdgid(itr->pdgId());
+	tau.setPx(recotau1->p4().Px());
+	tau.setPy(recotau1->p4().Py());
+	tau.setPz(recotau1->p4().Pz());
+	tau.setE (recotau1->p4().E());
+        tau.setPdgid(recotau1->pdgId());
+	GetRecoDaughters(recotau1,tau_daughters,recotau1->pdgId());
+	tau2.setPx(recotau2->p4().Px());
+        tau2.setPy(recotau2->p4().Py());
+        tau2.setPz(recotau2->p4().Pz());
+        tau2.setE (recotau2->p4().E());
+        tau2.setPdgid(recotau2->pdgId());
+	if(ntau==2)GetRecoDaughters(recotau2,tau2_daughters,recotau2->pdgId());
+	return 0;
+      }
+    }
+  }
+  return 1;
+}
+
+void TauSpinnerCMS::GetLastSelf(const reco::GenParticle *Particle){
+  for (unsigned int i=0; i< Particle->numberOfDaughters(); i++){
+    const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(Particle->daughter(i));
+    if(Particle->pdgId()==dau->pdgId()){
+      Particle=dau;
+      GetLastSelf(Particle);  
+    }
+  }
+}
+
+bool TauSpinnerCMS::isFirst(const reco::GenParticle *Particle){
+  for (unsigned int i=0; i< Particle->numberOfMothers(); i++){
+    const reco::GenParticle *moth=static_cast<const reco::GenParticle*>(Particle->mother(i));
+    if(Particle->pdgId()==moth->pdgId()){
+      return false;
+    }
+  }
+  return true;
+}
+
+void TauSpinnerCMS::GetRecoDaughters(const reco::GenParticle *Particle,std::vector<SimpleParticle> &daughters, int parentpdgid){
+  if( Particle->numberOfDaughters()==0 || abs(Particle->pdgId())==111){
+    SimpleParticle tp(Particle->p4().Px(), Particle->p4().Py(), Particle->p4().Pz(), Particle->p4().E(), Particle->pdgId());
+    daughters.push_back(tp);
+    return;
+  }
+  for (unsigned int i=0; i< Particle->numberOfDaughters(); i++){
+    const reco::Candidate *dau=Particle->daughter(i);
+    GetRecoDaughters(static_cast<const reco::GenParticle*>(dau),daughters,Particle->pdgId());
+  }
+}
+
+DEFINE_FWK_MODULE(TauSpinnerCMS);