#3 dedicated ntuple analyzer and its template

GenericNtuple/interface/KDileptonTreeReducer.h

@@ -0,0 +1,36 @@
+#ifndef KrAFT_GenericNtuple_KDileptonTreeReducer_H
+#define KrAFT_GenericNtuple_KDileptonTreeReducer_H
+
+#include "KrAFT/GenericNtuple/interface/KFlatTreeReducerBase.h"
+
+class KDileptonTreeReducer : public KFlatTreeReducerBase
+{
+public:
+  KDileptonTreeReducer(const std::string modeName,
+                       const std::string inputFileName,
+                       const std::string outputFileName)
+    :KFlatTreeReducerBase(modeName, inputFileName, outputFileName) {}
+  void run() { KFlatTreeReducerBase::run(); }
+
+private:
+  virtual void init();
+  virtual bool analyze();
+
+private:
+  // Cache for Leptons, depend on decay modes
+  typedef doubles* doublesP;
+  typedef ints* intsP;
+  doublesP leptons1_pt_, leptons1_eta_, leptons1_phi_, leptons1_m_, leptons1_iso_;
+  doublesP leptons2_pt_, leptons2_eta_, leptons2_phi_, leptons2_m_, leptons2_iso_;
+  intsP leptons1_Q_, leptons2_Q_;
+  int mode_;
+
+private:
+  // Branch items for output tree
+  double lepton1_pt_, lepton1_eta_, lepton1_phi_, lepton1_iso_;
+  double lepton2_pt_, lepton2_eta_, lepton2_phi_, lepton2_iso_;
+  double z_m_, z_pt_;
+  int z_Q_;
+};
+
+#endif

GenericNtuple/interface/KFlatTreeReducerBase.h

@@ -0,0 +1,52 @@
+#ifndef KFlatTreeReducerBase_h
+#define KFlatTreeReducerBase_h
+
+#include <iostream>
+
+#include "TFile.h"
+#include "TH1.h"
+#include "TTree.h"
+#include "TDirectory.h"
+#include "TMath.h"
+
+#include "KrAFT/GenericNtuple/interface/GenericEvent.h"
+
+//#include "Math/GenVector/LorentzVector.h"
+//#include "Math/GenVector/VectorUtil.h"
+//using namespace ROOT::Math::VectorUtil;
+//typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > LorentzVector;
+
+#include <TLorentzVector.h>
+typedef TLorentzVector LorentzVector;
+typedef std::vector<LorentzVector> LorentzVectors;
+typedef std::vector<double> doubles;
+typedef std::vector<int> ints;
+
+void printEntryFraction(int i, int n);
+
+class KFlatTreeReducerBase
+{
+public:
+  KFlatTreeReducerBase(const std::string modeName,
+                       const std::string inputFileName, 
+                       const std::string outputFileName);
+  virtual ~KFlatTreeReducerBase();
+
+  void run();
+
+protected:
+  virtual void init() = 0;
+  virtual bool analyze() = 0;
+  void getP4(const doubles* pt, const doubles* eta, const doubles* phi, const doubles* m, LorentzVectors& p4);
+  std::string modeName_;
+  TFile* inputFile_;
+  TFile* outputFile_;
+  GenericEvent* event_;
+  TDirectory* outDir_;
+  TH1* hEvent_;
+  TTree* outTree_;
+
+  bool isMC_;
+};
+
+#endif

GenericNtuple/src/KDileptonTreeReducer.cc

@@ -0,0 +1,116 @@
+#include "KrAFT/GenericNtuple/interface/KDileptonTreeReducer.h"
+
+using namespace std;
+
+void KDileptonTreeReducer::init()
+{
+  if ( modeName_ == "ElEl")
+  {
+    mode_ = 1;
+    leptons1_pt_  = leptons2_pt_  = event_->electrons_pt_ ;
+    leptons1_eta_ = leptons2_eta_ = event_->electrons_eta_;
+    leptons1_phi_ = leptons2_phi_ = event_->electrons_phi_;
+    leptons1_m_   = leptons2_m_   = event_->electrons_m_  ;
+    leptons1_Q_   = leptons2_Q_   = event_->electrons_Q_  ;
+    leptons1_iso_ = leptons2_iso_ = event_->electrons_relIso_;
+  }
+  else if ( modeName_ == "MuMu" )
+  {
+    mode_ = 2;
+    leptons1_pt_  = leptons2_pt_  = event_->muons_pt_ ;
+    leptons1_eta_ = leptons2_eta_ = event_->muons_eta_;
+    leptons1_phi_ = leptons2_phi_ = event_->muons_phi_;
+    leptons1_m_   = leptons2_m_   = event_->muons_m_  ;
+    leptons1_Q_   = leptons2_Q_   = event_->muons_Q_  ;
+    leptons1_iso_ = leptons2_iso_ = event_->muons_relIso_;
+  }
+  else
+  {
+    mode_ = 3;
+
+    leptons1_pt_  = event_->electrons_pt_ ;
+    leptons1_eta_ = event_->electrons_eta_;
+    leptons1_phi_ = event_->electrons_phi_;
+    leptons1_m_   = event_->electrons_m_  ;
+    leptons1_Q_   = event_->electrons_Q_  ;
+    leptons1_iso_ = event_->electrons_relIso_;
+
+    leptons2_pt_  = event_->muons_pt_ ;
+    leptons2_eta_ = event_->muons_eta_;
+    leptons2_phi_ = event_->muons_phi_;
+    leptons2_m_   = event_->muons_m_  ;
+    leptons2_Q_   = event_->muons_Q_  ;
+    leptons2_iso_ = event_->muons_relIso_;
+  }
+
+  outTree_->Branch("lepton1_pt" , &lepton1_pt_ , "lepton1_pt/D" );
+  outTree_->Branch("lepton1_eta", &lepton1_eta_, "lepton1_eta/D");
+  outTree_->Branch("lepton1_phi", &lepton1_phi_, "lepton1_phi/D");
+  outTree_->Branch("lepton1_iso", &lepton1_iso_, "lepton1_iso/D");
+
+  outTree_->Branch("lepton2_pt" , &lepton2_pt_ , "lepton2_pt/D" );
+  outTree_->Branch("lepton2_eta", &lepton2_eta_, "lepton2_eta/D");
+  outTree_->Branch("lepton2_phi", &lepton2_phi_, "lepton2_phi/D");
+  outTree_->Branch("lepton2_iso", &lepton2_iso_, "lepton2_iso/D");
+
+  outTree_->Branch("z_m" , &z_m_ , "z_m/D" );
+  outTree_->Branch("z_pt", &z_pt_, "z_pt/D");
+  outTree_->Branch("z_Q" , &z_Q_ , "z_Q/I" );
+}
+
+bool KDileptonTreeReducer::analyze()
+{
+  // Select leptons
+  LorentzVector lepton1P4, lepton2P4, zP4;
+  for ( int i=0, n=leptons1_pt_->size(); i<n; ++i )
+  {
+    if ( mode_ == 1 )
+    {
+      if ( event_->electrons_type_->at(i) < 100 ) continue; // ElEl mode
+    }
+    else if ( event_->muons_type_->at(i) == 0 ) continue; // MuMu or MuEl
+
+    lepton1_pt_  = leptons1_pt_ ->at(i);
+    lepton1_eta_ = leptons1_eta_->at(i);
+    if ( lepton1_pt_ < 20 or std::abs(lepton1_eta_) > 2.5 ) continue;
+
+    lepton1_phi_ = leptons1_phi_->at(i);
+    const double lepton1M = leptons1_m_->at(i);
+    const int lepton1Q = leptons1_Q_->at(i);
+
+    lepton1P4.SetPtEtaPhiM(lepton1_pt_, lepton1_eta_, lepton1_phi_, lepton1M);
+    lepton1_iso_ = leptons1_iso_->at(i);
+
+    for ( int j=((mode_ == 3) ? 0 : i+1); j<n; ++j )
+    {
+      if ( mode_ == 2 )
+      {
+        if ( event_->muons_type_->at(j) == 0 ) continue; // MuMu mode
+      }
+      else if ( event_->electrons_type_->at(j) < 100 ) continue; // ElEl and MuEl
+
+      lepton2_pt_  = leptons2_pt_->at(i);
+      lepton2_eta_ = leptons2_eta_->at(i);
+      if ( lepton2_pt_ < 20 or std::abs(lepton2_eta_) > 2.5 ) continue;
+
+      lepton2_phi_ = leptons2_phi_->at(j);
+      const double lepton2M = leptons2_m_->at(j);
+      const int lepton2Q = leptons2_Q_->at(j);
+
+      lepton2P4.SetPtEtaPhiM(lepton2_pt_, lepton2_eta_, lepton2_phi_, lepton2M);
+      lepton2_iso_ = leptons2_iso_->at(j);
+
+      zP4 = lepton1P4 + lepton2P4;
+      z_Q_ = lepton1Q + lepton2Q;
+
+      break;
+    }
+    break;
+  }
+  if ( z_Q_ == -999 ) return false;
+  z_pt_ = zP4.Pt();
+  z_m_  = zP4.M();
+
+  return true;
+}
+

GenericNtuple/src/KFlatTreeReducer.cc → GenericNtuple/src/KFlatTreeReducerBase.cc

@@ -1,12 +1,9 @@
-#include "KrAFT/GenericNtuple/interface/KFlatTreeReducer.h"
+#include "KrAFT/GenericNtuple/interface/KFlatTreeReducerBase.h"
 #include "TString.h"
 #include "TNamed.h"
 
 using namespace std;
 
-using namespace ROOT::Math::VectorUtil;
-typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > LorentzVector;
-
 void printEntryFraction(int i, int n)
 {
   static int digit = 1;
@@ -24,10 +21,9 @@ void printEntryFraction(int i, int n)
   }
 }
 
-KFlatTreeReducer::KFlatTreeReducer(const string modeName, const string inputFileName, const string outputFileName)
+KFlatTreeReducerBase::KFlatTreeReducerBase(const string modeName, const string inputFileName, const string outputFileName)
 {
-  eventScale_ = 1;
-
+  modeName_ = modeName;
   inputFile_ = TFile::Open(inputFileName.c_str());
   TTree* tree = dynamic_cast<TTree*>(inputFile_->Get(Form("%s/event", modeName.c_str())));
   const string dataType = inputFile_->Get(Form("%s/dataType", modeName.c_str()))->GetTitle();
@@ -42,37 +38,41 @@ KFlatTreeReducer::KFlatTreeReducer(const string modeName, const string inputFile
 
   outputFile_ = TFile::Open(outputFileName.c_str(), "RECREATE");
   outputFile_->cd();
-  TDirectory* outDir = outputFile_->mkdir(modeName.c_str());
+  outputFile_->mkdir(modeName.c_str());
   outTree_ = new TTree("ntuple", "ntuple");
 }
 
-void KFlatTreeReducer::setEventScale(const double scale)
-{
-  eventScale_ = scale;
-}
-
-void KFlatTreeReducer::setCrossSection(const double crossSection)
-{
-  if ( not isMC_ ) return;
-  const double nEventGen = hEvent_->GetBinContent(1);
-  eventScale_ = crossSection/nEventGen;
-}
-
-KFlatTreeReducer::~KFlatTreeReducer()
+KFlatTreeReducerBase::~KFlatTreeReducerBase()
 {
   //outputFile_->cd();
-  outTree_->Write();
+  //outTree_->Write();
   //outputFile_->Write();
   //outputFile_->Close();
 }
 
-void KFlatTreeReducer::analyze()
+void KFlatTreeReducerBase::run()
 {
+  init();
   const int nEvent = event_->tree_->GetEntries();
   for ( int i=0; i<nEvent; ++i )
   {
     printEntryFraction(i, nEvent);
     event_->tree_->GetEntry(i);
+
+    const bool isAccepted = analyze();
+    if ( !isAccepted ) continue;
+
+    outTree_->Fill();
   }
 }
 
+void KFlatTreeReducerBase::getP4(const doubles* pts, const doubles* etas, const doubles* phis, const doubles* ms, LorentzVectors& p4s)
+{
+  p4s.clear();
+  LorentzVector p4;
+  for ( int i=0, n=pts->size(); i<n; ++i )
+  {
+    p4.SetPtEtaPhiM(pts->at(i), etas->at(i), phis->at(i), ms->at(i));
+    p4s.push_back(p4);
+  }
+}

GenericNtuple/src/classes.h

@@ -1,9 +1,9 @@
-#include "KrAFT/GenericNtuple/interface/KFlatTreeReducer.h"
+#include "KrAFT/GenericNtuple/interface/KDileptonTreeReducer.h"
 #include "DataFormats/Common/interface/Wrapper.h"
 
 namespace {
   struct dictionary {
-    KFlatTreeReducer kFlatTreeReducer;
+    KDileptonTreeReducer kDileptonTreeReducer;
   };
 
 }

GenericNtuple/src/classes_def.xml

@@ -1,4 +1,4 @@
 <lcgdict>
-  <class name="KFlatTreeReducer" />
+  <class name="KDileptonTreeReducer" />
 </lcgdict>
 

GenericNtuple/test/run.py

@@ -10,12 +10,7 @@
 mode = "ElEl"
 sample = "DYJetsToLL_M-10To50filter_8TeV-madgraph"
 
-f = TFile("/pnfs/user/jhgoh/data/ntuple/20131224_1/%s.root" % sample)
-tree = f.Get("%s/event" % mode)
-#tree.Draw("nVertex", "(puWeight)*(1)")
-tree.Draw("../src/m%s.cc" % mode, "../src/cut.cc")
-
-#ana = KFlatTreeAnalyzer("ElEl",
-#    "/pnfs/user/jhgoh/data/ntuple/20131224_1/%s.root" % sample,
-#    "hist/%s__%s.root" % (sample, mode))
-#ana.analyze()
+ana = KDileptonTreeReducer("ElEl",
+    "/pnfs/user/jhgoh/data/ntuple/20131224_1/%s.root" % sample,
+    "hist/%s__%s.root" % (sample, mode))
+ana.run()