Index: Makefile
===================================================================
--- Makefile	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
+++ Makefile	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
@@ -322,4 +322,5 @@
 	modules/UniqueObjectFinder.h \
 	modules/TrackCountingBTagging.h \
+	modules/BTaggingCMS.h \
 	modules/BTagging.h \
 	modules/TauTagging.h \
@@ -339,5 +340,5 @@
 	modules/Hector.h \
 	modules/RunPUPPI.h \
-	modules/JetFlavourAssociation.h \
+	modules/JetFlavorAssociation.h \
 	modules/ExampleModule.h
 ModulesDict$(PcmSuf): \
@@ -545,4 +546,13 @@
 	external/ExRootAnalysis/ExRootFilter.h \
 	external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/BTaggingCMS.$(ObjSuf): \
+	modules/BTaggingCMS.$(SrcSuf) \
+	modules/BTaggingCMS.h \
+	classes/DelphesClasses.h \
+	classes/DelphesFactory.h \
+	classes/DelphesFormula.h \
+	external/ExRootAnalysis/ExRootResult.h \
+	external/ExRootAnalysis/ExRootFilter.h \
+	external/ExRootAnalysis/ExRootClassifier.h
 tmp/modules/Calorimeter.$(ObjSuf): \
 	modules/Calorimeter.$(SrcSuf) \
@@ -658,7 +668,7 @@
 	external/ExRootAnalysis/ExRootFilter.h \
 	external/ExRootAnalysis/ExRootClassifier.h
-tmp/modules/JetFlavourAssociation.$(ObjSuf): \
-	modules/JetFlavourAssociation.$(SrcSuf) \
-	modules/JetFlavourAssociation.h \
+tmp/modules/JetFlavorAssociation.$(ObjSuf): \
+	modules/JetFlavorAssociation.$(SrcSuf) \
+	modules/JetFlavorAssociation.h \
 	classes/DelphesClasses.h \
 	classes/DelphesFactory.h \
@@ -897,4 +907,5 @@
 	tmp/modules/AngularSmearing.$(ObjSuf) \
 	tmp/modules/BTagging.$(ObjSuf) \
+	tmp/modules/BTaggingCMS.$(ObjSuf) \
 	tmp/modules/Calorimeter.$(ObjSuf) \
 	tmp/modules/Cloner.$(ObjSuf) \
@@ -909,5 +920,5 @@
 	tmp/modules/ImpactParameterSmearing.$(ObjSuf) \
 	tmp/modules/Isolation.$(ObjSuf) \
-	tmp/modules/JetFlavourAssociation.$(ObjSuf) \
+	tmp/modules/JetFlavorAssociation.$(ObjSuf) \
 	tmp/modules/JetPileUpSubtractor.$(ObjSuf) \
 	tmp/modules/LeptonDressing.$(ObjSuf) \
@@ -1762,5 +1773,5 @@
 	@touch $@
 
-modules/JetFlavourAssociation.h: \
+modules/BTaggingCMS.h: \
 	classes/DelphesModule.h \
 	classes/DelphesClasses.h
@@ -1779,4 +1790,9 @@
 	external/fastjet/ClusterSequenceAreaBase.hh \
 	external/fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh
+	@touch $@
+
+modules/JetFlavorAssociation.h: \
+	classes/DelphesModule.h \
+	classes/DelphesClasses.h
 	@touch $@
 
Index: modules/BTaggingCMS.cc
===================================================================
--- modules/BTaggingCMS.cc	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
+++ modules/BTaggingCMS.cc	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
@@ -0,0 +1,421 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+/** \class BTaggingCMS
+ *
+ *  Applies b-tagging efficiency (miss identification rate) formulas
+ *  and sets b-tagging flags
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/BTaggingCMS.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+BTaggingCMS::BTaggingCMS() :
+  fItJetInputArray(0)
+{
+}
+
+//------------------------------------------------------------------------------
+
+BTaggingCMS::~BTaggingCMS()
+{
+}
+
+//------------------------------------------------------------------------------
+
+void BTaggingCMS::Init()
+{
+  map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
+  ExRootConfParam param;
+  DelphesFormula *formula;
+  Int_t i, size;
+
+  // read efficiency formulas
+  param = GetParam("EfficiencyFormulaLoose");
+  size = param.GetSize();
+
+  fEfficiencyMapLoose.clear();
+  for(i = 0; i < size/2; ++i)
+  {
+    formula = new DelphesFormula;
+    formula->Compile(param[i*2 + 1].GetString());
+    fEfficiencyMapLoose[param[i*2].GetInt()] = formula;
+  }
+
+  // set default efficiency formula
+  itEfficiencyMap = fEfficiencyMapLoose.find(0);
+  if(itEfficiencyMap == fEfficiencyMapLoose.end())
+  {
+    formula = new DelphesFormula;
+    formula->Compile("0.0");
+    fEfficiencyMapLoose[0] = formula;
+  }
+
+  // read efficiency formulas
+  param = GetParam("EfficiencyFormulaMedium");
+  size = param.GetSize();
+
+  fEfficiencyMapMedium.clear();
+  for(i = 0; i < size/2; ++i)
+  {
+    formula = new DelphesFormula;
+    formula->Compile(param[i*2 + 1].GetString());
+    fEfficiencyMapMedium[param[i*2].GetInt()] = formula;
+  }
+
+  // set default efficiency formula
+  itEfficiencyMap = fEfficiencyMapMedium.find(0);
+  if(itEfficiencyMap == fEfficiencyMapMedium.end())
+  {
+    formula = new DelphesFormula;
+    formula->Compile("0.0");
+    fEfficiencyMapMedium[0] = formula;
+  }
+
+  // read efficiency formulas
+  param = GetParam("EfficiencyFormulaTight");
+  size = param.GetSize();
+
+  fEfficiencyMapTight.clear();
+  for(i = 0; i < size/2; ++i)
+  {
+    formula = new DelphesFormula;
+    formula->Compile(param[i*2 + 1].GetString());
+    fEfficiencyMapTight[param[i*2].GetInt()] = formula;
+  }
+
+  // set default efficiency formula
+  itEfficiencyMap = fEfficiencyMapTight.find(0);
+  if(itEfficiencyMap == fEfficiencyMapTight.end())
+  {
+    formula = new DelphesFormula;
+    formula->Compile("0.0");
+    fEfficiencyMapTight[0] = formula;
+  }
+
+  // import input array(s)
+  fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
+  fItJetInputArray = fJetInputArray->MakeIterator();
+}
+
+//------------------------------------------------------------------------------
+
+void BTaggingCMS::Finish()
+{
+  map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
+  DelphesFormula *formula;
+
+  if(fItJetInputArray) delete fItJetInputArray;
+
+  for(itEfficiencyMap = fEfficiencyMapLoose.begin(); itEfficiencyMap != fEfficiencyMapLoose.end(); ++itEfficiencyMap)
+  {
+    formula = itEfficiencyMap->second;
+    if(formula) delete formula;
+  }
+
+  for(itEfficiencyMap = fEfficiencyMapMedium.begin(); itEfficiencyMap != fEfficiencyMapMedium.end(); ++itEfficiencyMap)
+  {
+    formula = itEfficiencyMap->second;
+    if(formula) delete formula;
+  }
+
+  for(itEfficiencyMap = fEfficiencyMapTight.begin(); itEfficiencyMap != fEfficiencyMapTight.end(); ++itEfficiencyMap)
+  {
+    formula = itEfficiencyMap->second;
+    if(formula) delete formula;
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void BTaggingCMS::Process()
+{
+  Candidate *jet;
+  map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
+  DelphesFormula *formula;
+  float randomNumber;
+
+  // loop over all input jets
+  fItJetInputArray->Reset();
+
+  while((jet = static_cast<Candidate*>(fItJetInputArray->Next())))
+  {
+    const TLorentzVector &jetMomentum = jet->Momentum;
+
+    randomNumber = gRandom->Uniform();
+
+    // --------------------------
+    // Loose b-tagging
+    // --------------------------
+
+    // find heaviest flavor and b-tag
+    if(fEfficiencyMapLoose.size() != 0)
+    {
+      itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorHeaviest);
+      if(itEfficiencyMap == fEfficiencyMapLoose.end())
+      {
+        itEfficiencyMap = fEfficiencyMapLoose.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHeaviest = 1;
+
+      // find hisghestpt flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorHighestPt);
+      if(itEfficiencyMap == fEfficiencyMapLoose.end())
+      {
+        itEfficiencyMap = fEfficiencyMapLoose.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHighestPt = 1;
+
+      // find nearest2 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorNearest2);
+      if(itEfficiencyMap == fEfficiencyMapLoose.end())
+      {
+        itEfficiencyMap = fEfficiencyMapLoose.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest2 = 1;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorNearest3);
+      if(itEfficiencyMap == fEfficiencyMapLoose.end())
+      {
+        itEfficiencyMap = fEfficiencyMapLoose.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest3 = 1;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorAlgo);
+      if(itEfficiencyMap == fEfficiencyMapLoose.end())
+      {
+        itEfficiencyMap = fEfficiencyMapLoose.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagAlgo = 1;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorPhysics);
+      if(itEfficiencyMap == fEfficiencyMapLoose.end())
+      {
+        itEfficiencyMap = fEfficiencyMapLoose.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagPhysics = 1;
+
+      // default
+      itEfficiencyMap = fEfficiencyMapLoose.find(jet->FlavorDefault);
+      if(itEfficiencyMap == fEfficiencyMapLoose.end())
+      {
+        itEfficiencyMap = fEfficiencyMapLoose.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagDefault = 1;
+    }
+
+    // --------------------------
+    // Medium b-tagging
+    // --------------------------
+
+    // find heaviest flavor and b-tag
+    if(fEfficiencyMapMedium.size() != 0)
+    {
+      itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorHeaviest);
+      if(itEfficiencyMap == fEfficiencyMapMedium.end())
+      {
+        itEfficiencyMap = fEfficiencyMapMedium.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHeaviest = 2;
+
+      // find hisghestpt flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorHighestPt);
+      if(itEfficiencyMap == fEfficiencyMapMedium.end())
+      {
+        itEfficiencyMap = fEfficiencyMapMedium.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHighestPt = 2;
+
+      // find nearest2 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorNearest2);
+      if(itEfficiencyMap == fEfficiencyMapMedium.end())
+      {
+        itEfficiencyMap = fEfficiencyMapMedium.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest2 = 2;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorNearest3);
+      if(itEfficiencyMap == fEfficiencyMapMedium.end())
+      {
+        itEfficiencyMap = fEfficiencyMapMedium.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest3 = 2;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorAlgo);
+      if(itEfficiencyMap == fEfficiencyMapMedium.end())
+      {
+        itEfficiencyMap = fEfficiencyMapMedium.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagAlgo = 2;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorPhysics);
+      if(itEfficiencyMap == fEfficiencyMapMedium.end())
+      {
+        itEfficiencyMap = fEfficiencyMapMedium.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagPhysics = 2;
+
+      // default
+      itEfficiencyMap = fEfficiencyMapMedium.find(jet->FlavorDefault);
+      if(itEfficiencyMap == fEfficiencyMapMedium.end())
+      {
+        itEfficiencyMap = fEfficiencyMapMedium.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagDefault = 2;
+    }
+
+    // --------------------------
+    // Tight b-tagging
+    // --------------------------
+
+    // find heaviest flavor and b-tag
+    if(fEfficiencyMapTight.size() != 0)
+    {
+      itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorHeaviest);
+      if(itEfficiencyMap == fEfficiencyMapTight.end())
+      {
+        itEfficiencyMap = fEfficiencyMapTight.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHeaviest = 3;
+
+      // find hisghestpt flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorHighestPt);
+      if(itEfficiencyMap == fEfficiencyMapTight.end())
+      {
+        itEfficiencyMap = fEfficiencyMapTight.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagHighestPt = 3;
+
+      // find nearest2 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorNearest2);
+      if(itEfficiencyMap == fEfficiencyMapTight.end())
+      {
+        itEfficiencyMap = fEfficiencyMapTight.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest2 = 3;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorNearest3);
+      if(itEfficiencyMap == fEfficiencyMapTight.end())
+      {
+        itEfficiencyMap = fEfficiencyMapTight.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagNearest3 = 3;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorAlgo);
+      if(itEfficiencyMap == fEfficiencyMapTight.end())
+      {
+        itEfficiencyMap = fEfficiencyMapTight.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagAlgo = 3;
+
+      // find nearest3 flavor and b-tag
+      itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorPhysics);
+      if(itEfficiencyMap == fEfficiencyMapTight.end())
+      {
+        itEfficiencyMap = fEfficiencyMapTight.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagPhysics = 3;
+
+      // default
+      itEfficiencyMap = fEfficiencyMapTight.find(jet->FlavorDefault);
+      if(itEfficiencyMap == fEfficiencyMapTight.end())
+      {
+        itEfficiencyMap = fEfficiencyMapTight.find(0);
+      }
+      formula = itEfficiencyMap->second;
+
+      if(randomNumber <= formula->Eval(jetMomentum.Pt(), jetMomentum.Eta())) jet->BTagDefault = 3;
+    }
+  }
+}
Index: modules/BTaggingCMS.h
===================================================================
--- modules/BTaggingCMS.h	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
+++ modules/BTaggingCMS.h	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
@@ -0,0 +1,68 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef BTaggingCMS_h
+#define BTaggingCMS_h
+
+/** \class BTaggingCMS
+ *
+ *  Applies b-tagging efficiency (miss identification rate) formulas
+ *  and sets b-tagging flags
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesModule.h"
+#include "classes/DelphesClasses.h"
+
+#include <map>
+
+class TObjArray;
+class DelphesFormula;
+
+class ExRootFilter;
+
+class BTaggingCMS: public DelphesModule
+{
+
+public:
+
+  BTaggingCMS();
+  ~BTaggingCMS();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+#if !defined(__CINT__) && !defined(__CLING__)
+  std::map< Int_t, DelphesFormula * > fEfficiencyMapLoose; //!
+  std::map< Int_t, DelphesFormula * > fEfficiencyMapMedium; //!
+  std::map< Int_t, DelphesFormula * > fEfficiencyMapTight; //!
+#endif
+
+  TIterator *fItJetInputArray; //!
+
+  const TObjArray *fJetInputArray; //!
+
+  ClassDef(BTaggingCMS, 1)
+};
+
+#endif
Index: modules/JetFlavorAssociation.cc
===================================================================
--- modules/JetFlavorAssociation.cc	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
+++ modules/JetFlavorAssociation.cc	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
@@ -0,0 +1,399 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY || FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+/** \class JetFlavorAssociation
+ *
+ *  Find origin of jet && evaluate jet flavor
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/JetFlavorAssociation.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+class PartonClassifier: public ExRootClassifier
+{
+public:
+
+  PartonClassifier() {}
+  Int_t GetCategory(TObject *object);
+  Double_t fEtaMax, fPTMin;
+};
+
+//------------------------------------------------------------------------------
+// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/PartonSelector.cc
+
+Int_t PartonClassifier::GetCategory(TObject *object)
+{
+ // select parton in the parton list
+
+  Candidate *parton = static_cast<Candidate *>(object);
+  const TLorentzVector &momentum = parton->Momentum;
+  Int_t pdgCode;
+
+  // inside the eta && momentum range (be a little bit larger that the tracking coverage
+  if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
+
+  pdgCode = TMath::Abs(parton->PID);
+
+  if(parton->Status == -1) return -1;
+  if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
+  if(parton->Status == 3 || parton->Status == 2) return 0; // if status 3 return
+
+  return 0;
+}
+
+//------------------------------------------------------------------------------
+
+class PartonClassifierLHEF: public ExRootClassifier
+{
+public:
+
+  PartonClassifierLHEF() {}
+  Int_t GetCategory(TObject *object);
+  Double_t fEtaMax, fPTMin;
+};
+
+Int_t PartonClassifierLHEF::GetCategory(TObject *object)
+{
+  // select parton in the parton list
+
+  Candidate *parton = static_cast<Candidate *>(object);
+  const TLorentzVector &momentum = parton->Momentum;
+  Int_t pdgCode;
+
+  // inside the eta && momentum range (be a little bit larger that the tracking coverage
+  if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
+
+  pdgCode = TMath::Abs(parton->PID);
+  if(parton->Status == -1) return -1;
+  if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
+  if(parton->Status != 1) return -1; // if status 3 return
+
+  return 0;
+}
+
+//------------------------------------------------------------------------------
+
+JetFlavorAssociation::JetFlavorAssociation() :
+  fClassifier(0), fFilter(0),
+  fItPartonInputArray(0), fItPartonInputArrayLHEF(0),
+  fItJetInputArray(0), fItParticleInputArray(0)
+{
+  fClassifier = new PartonClassifier;
+  fClassifierLHEF = new PartonClassifierLHEF;
+}
+
+//------------------------------------------------------------------------------
+
+JetFlavorAssociation::~JetFlavorAssociation()
+{
+  if(fClassifier) delete fClassifier;
+  if(fClassifierLHEF) delete fClassifierLHEF;
+}
+
+//------------------------------------------------------------------------------
+
+void JetFlavorAssociation::Init()
+{
+  ExRootConfParam param;
+
+  fDeltaR = GetDouble("DeltaR", 0.5);
+
+  fClassifier->fPTMin = GetDouble("PartonPTMin", 0.);
+  fClassifier->fEtaMax = GetDouble("PartonEtaMax",2.5);
+
+  fClassifierLHEF->fPTMin = GetDouble("PartonPTMin", 0.);
+  fClassifierLHEF->fEtaMax = GetDouble("PartonEtaMax",2.5);
+
+  // import input array(s)
+  fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
+  fItPartonInputArray = fPartonInputArray->MakeIterator();
+
+  fPartonInputArrayLHEF = ImportArray(GetString("PartonInputArrayLHEF", "Delphes/partonsLHEF"));
+  fItPartonInputArrayLHEF = fPartonInputArrayLHEF->MakeIterator();
+
+  fFilter = new ExRootFilter(fPartonInputArray);
+  fFilterLHEF = new ExRootFilter(fPartonInputArrayLHEF);
+
+  fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
+  fItJetInputArray = fJetInputArray->MakeIterator();
+
+  fParticleInputArray = ImportArray(GetString("ParticleInputArray", "Delphes/allParticles"));
+  fItParticleInputArray = fParticleInputArray->MakeIterator();
+
+}
+
+//------------------------------------------------------------------------------
+
+void JetFlavorAssociation::Finish()
+{
+  if(fFilter) delete fFilter;
+  if(fFilterLHEF) delete fFilterLHEF;
+
+  if(fItJetInputArray) delete fItJetInputArray;
+  if(fItParticleInputArray) delete fItParticleInputArray;
+  if(fItPartonInputArray) delete fItPartonInputArray;
+  if(fItPartonInputArrayLHEF) delete fItPartonInputArrayLHEF;
+}
+
+//------------------------------------------------------------------------------
+
+void JetFlavorAssociation::Process(){
+
+  Candidate *jet;
+  TObjArray *partonArray;
+  TObjArray *partonArrayLHEF;
+
+  // select quark && gluons
+  fFilter->Reset();
+  partonArray = fFilter->GetSubArray(fClassifier, 0); // get the filtered parton array
+
+  if(partonArray == 0) return;
+  TIter itPartonArray(partonArray);
+
+  fFilterLHEF->Reset();
+  partonArrayLHEF = fFilterLHEF->GetSubArray(fClassifierLHEF, 0); // get the filtered parton array
+
+  if(partonArrayLHEF == 0) return;
+  TIter itPartonLHEFArray(partonArrayLHEF);
+
+  // loop over all input jets
+  fItJetInputArray->Reset();
+  while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
+  {
+    // get standard flavor
+    GetAlgoFlavor(jet, itPartonArray, itPartonLHEFArray);
+    GetPhysicsFlavor(jet, itPartonArray, itPartonLHEFArray);
+  }
+}
+
+//------------------------------------------------------------------------------
+// Standard definition of jet flavor in
+// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/JetPartonMatcher.cc?v=CMSSW_7_3_0_pre1
+
+void JetFlavorAssociation::GetAlgoFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
+{
+  float maxPt = 0;
+  float minDr = 1000;
+  bool isGoodParton = true;
+  int daughterCounter = 0;
+  Candidate *parton, *partonLHEF;
+  Candidate *tempParton = 0, *tempPartonNearest = 0, *tempPartonHighestPt = 0;
+  int pdgCode, pdgCodeMax = -1;
+
+  itPartonArray.Reset();
+  while((parton = static_cast<Candidate *>(itPartonArray.Next())))
+  {
+    // default delphes method
+    pdgCode = TMath::Abs(parton->PID);
+    if(TMath::Abs(parton->PID) == 21) pdgCode = 0;
+    if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
+    {
+      if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
+    }
+
+    isGoodParton = true;
+    itPartonLHEFArray.Reset();
+    while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
+    {
+      if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.001 &&
+         parton->PID == partonLHEF->PID &&
+         partonLHEF->Charge == parton->Charge)
+      {
+         isGoodParton = false;
+         break;
+      }
+
+      if(!isGoodParton) continue;
+
+      // check the daugheter
+      daughterCounter = 0;
+      if(parton->D1 != -1 || parton->D2 != -1)
+      {
+        // partons are only quarks || gluons
+        int daughterFlavor1 = -1;
+        int daughterFlavor2 = -1;
+        if(parton->D1 != -1) daughterFlavor1 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D1))->PID);
+        if(parton->D2 != -1) daughterFlavor2 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D2))->PID);
+        if((daughterFlavor1 == 1 || daughterFlavor1 == 2 || daughterFlavor1 == 3 || daughterFlavor1 == 4 || daughterFlavor1 == 5 || daughterFlavor1 == 21)) daughterCounter++;
+        if((daughterFlavor2 == 1 || daughterFlavor2 == 2 || daughterFlavor2 == 3 || daughterFlavor2 == 4 || daughterFlavor1 == 5 || daughterFlavor2 == 21)) daughterCounter++;
+      }
+      if(daughterCounter > 0) continue;
+      if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
+      {
+        if(jet->Momentum.DeltaR(parton->Momentum) < minDr)
+        {
+          minDr = jet->Momentum.DeltaR(parton->Momentum);
+          tempPartonNearest = parton;
+        }
+
+        // if not yet found && pdgId is a c, take as c
+        if(TMath::Abs(parton->PID) == 4) tempParton = parton;
+        if(TMath::Abs(parton->PID) == 5) tempParton = parton;
+        if(parton->Momentum.Pt() > maxPt)
+        {
+          maxPt = parton->Momentum.Pt();
+          tempPartonHighestPt = parton;
+        }
+      }
+    }
+  }
+
+  jet->FlavorHeaviest = tempParton ? TMath::Abs(tempParton->PID) : 0;
+  jet->FlavorHighestPt = tempPartonHighestPt ? TMath::Abs(tempPartonHighestPt->PID) : 0;
+  jet->FlavorNearest2 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
+  if(!tempParton) tempParton = tempPartonHighestPt;
+  jet->FlavorAlgo = tempParton ? TMath::Abs(tempParton->PID) : 0;
+
+  if(pdgCodeMax == 0) pdgCodeMax = 21;
+  if(pdgCodeMax == -1) pdgCodeMax = 0;
+
+  jet->FlavorDefault = pdgCodeMax;
+}
+
+//------------------------------------------------------------------------------
+
+void JetFlavorAssociation::GetPhysicsFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
+{
+  float minDr = 1000;
+  int partonCounter = 0;
+  float biggerConeSize = 0.7;
+  float dist;
+  bool isGoodCandidate;
+  int contaminatingFlavor = 0;
+  int motherCounter = 0;
+  Candidate *parton, *partonLHEF, *mother1, *mother2;
+  Candidate *tempParton = 0, *tempPartonNearest = 0;
+  vector<Candidate *> contaminations;
+  vector<Candidate *>::iterator itContaminations;
+
+  contaminations.clear();
+
+  itPartonLHEFArray.Reset();
+  while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
+  {
+    dist = jet->Momentum.DeltaR(partonLHEF->Momentum); // take the DR
+    if(partonLHEF->Status == 1 && dist < minDr)
+    {
+      tempPartonNearest = partonLHEF;
+      minDr = dist;
+    }
+
+    if(partonLHEF->Status == 1 && dist <= fDeltaR)
+    {
+      tempParton = partonLHEF;
+      partonCounter++;
+    }
+  }
+
+  itPartonArray.Reset();
+  itPartonLHEFArray.Reset();
+  while((parton = static_cast<Candidate *>(itPartonArray.Next())))
+  {
+    dist = jet->Momentum.DeltaR(parton->Momentum); // take the DR
+    isGoodCandidate = true;
+    while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
+    {
+      if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.01 &&
+         parton->PID == partonLHEF->PID &&
+         partonLHEF->Charge == parton->Charge)
+      {
+        isGoodCandidate = false;
+        break;
+      }
+    }
+
+    if(!isGoodCandidate) continue;
+
+    if(parton->D1 != -1 || parton->D2 != -1)
+    {
+      if((TMath::Abs(parton->PID) < 4 || TMath::Abs(parton->PID) == 21)) continue;
+      if(dist < biggerConeSize) contaminations.push_back(parton);
+    }
+  }
+
+  jet->FlavorNearest3 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
+
+  if(partonCounter != 1)
+  {
+    jet->FlavorPhysics = 0;
+  }
+  else if(contaminations.size() == 0)
+  {
+    jet->FlavorPhysics = TMath::Abs(tempParton->PID);
+  }
+  else if(contaminations.size() > 0)
+  {
+    jet->FlavorPhysics = TMath::Abs(tempParton->PID);
+
+    for(itContaminations = contaminations.begin(); itContaminations != contaminations.end(); ++itContaminations)
+    {
+      parton = *itContaminations;
+      contaminatingFlavor = TMath::Abs(parton->PID);
+      motherCounter = 0;
+      if(parton->M1 != -1) motherCounter++;
+      if(parton->M2 != -1) motherCounter++;
+
+      if(parton->M1 != -1)
+      {
+        mother1 = static_cast<Candidate *>(fParticleInputArray->At(parton->M1));
+        if(mother1 && motherCounter > 0 && mother1->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
+      }
+      if(parton->M2 != -1)
+      {
+        mother2 = static_cast<Candidate *>(fParticleInputArray->At(parton->M2));
+        if(mother2 && motherCounter > 0 && mother2->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
+      }
+      // mother is the initialParton --> OK
+      if(TMath::Abs(tempParton->PID) == 4)
+      {
+        // keep association --> the initialParton is a c --> the contaminated parton is a c
+        if(contaminatingFlavor == 4) continue;
+        jet->FlavorPhysics = 0; // all the other cases reject!
+        break;
+      }
+    }
+  }
+}
Index: modules/JetFlavorAssociation.h
===================================================================
--- modules/JetFlavorAssociation.h	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
+++ modules/JetFlavorAssociation.h	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
@@ -0,0 +1,78 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef JetFlavorAssociation_h
+#define JetFlavorAssociation_h
+
+/** \class JetFlavorAssociation
+ *
+ *  Find origin of jet and evaluate jet flavor
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesModule.h"
+#include "classes/DelphesClasses.h"
+#include <map>
+
+class TObjArray;
+class DelphesFormula;
+
+class ExRootFilter;
+class PartonClassifier;
+class PartonClassifierLHEF;
+
+class JetFlavorAssociation: public DelphesModule
+{
+public:
+
+  JetFlavorAssociation();
+  ~JetFlavorAssociation();
+
+  void Init();
+  void Process();
+  void Finish();
+
+  void GetAlgoFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonArrayLHEF);
+  void GetPhysicsFlavor(Candidate *jet, TIter &itPartonArray, TIter &itPartonArrayLHEF);
+
+private:
+
+  Double_t  fDeltaR;
+
+  PartonClassifier  *fClassifier; //!
+  PartonClassifierLHEF *fClassifierLHEF; //!
+
+  ExRootFilter *fFilter;
+  ExRootFilter *fFilterLHEF;
+
+  TIterator *fItPartonInputArray; //!
+  TIterator *fItPartonInputArrayLHEF; //!
+  TIterator *fItJetInputArray; //!
+  TIterator *fItParticleInputArray; //!
+
+  const TObjArray *fPartonInputArray; //!
+  const TObjArray *fPartonInputArrayLHEF; //!
+  const TObjArray *fJetInputArray; //!
+  const TObjArray *fParticleInputArray; //!
+
+  ClassDef(JetFlavorAssociation, 1)
+};
+
+#endif
Index: dules/JetFlavourAssociation.cc
===================================================================
--- modules/JetFlavourAssociation.cc	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
+++ 	(revision )
@@ -1,399 +1,0 @@
-/*
- *  Delphes: a framework for fast simulation of a generic collider experiment
- *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
- *
- *  This program is free software: you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation, either version 3 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY || FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-
-/** \class JetFlavourAssociation
- *
- *  Find origin of jet && evaluate jet flavor
- *
- *  \author P. Demin - UCL, Louvain-la-Neuve
- *
- */
-
-#include "modules/JetFlavourAssociation.h"
-
-#include "classes/DelphesClasses.h"
-#include "classes/DelphesFactory.h"
-#include "classes/DelphesFormula.h"
-
-#include "ExRootAnalysis/ExRootResult.h"
-#include "ExRootAnalysis/ExRootFilter.h"
-#include "ExRootAnalysis/ExRootClassifier.h"
-
-#include "TMath.h"
-#include "TString.h"
-#include "TFormula.h"
-#include "TRandom3.h"
-#include "TObjArray.h"
-#include "TDatabasePDG.h"
-#include "TLorentzVector.h"
-
-#include <algorithm>
-#include <stdexcept>
-#include <iostream>
-#include <sstream>
-
-using namespace std;
-
-//------------------------------------------------------------------------------
-
-class PartonClassifier: public ExRootClassifier
-{
-public:
-
-  PartonClassifier() {}
-  Int_t GetCategory(TObject *object);
-  Double_t fEtaMax, fPTMin;
-};
-
-//------------------------------------------------------------------------------
-// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/PartonSelector.cc
-
-Int_t PartonClassifier::GetCategory(TObject *object)
-{
- // select parton in the parton list
-
-  Candidate *parton = static_cast<Candidate *>(object);
-  const TLorentzVector &momentum = parton->Momentum;
-  Int_t pdgCode;
-
-  // inside the eta && momentum range (be a little bit larger that the tracking coverage
-  if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
-
-  pdgCode = TMath::Abs(parton->PID);
-
-  if(parton->Status == -1) return -1;
-  if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
-  if(parton->Status == 3 || parton->Status == 2) return 0; // if status 3 return
-
-  return 0;
-}
-
-//------------------------------------------------------------------------------
-
-class PartonClassifierLHEF: public ExRootClassifier
-{
-public:
-
-  PartonClassifierLHEF() {}
-  Int_t GetCategory(TObject *object);
-  Double_t fEtaMax, fPTMin;
-};
-
-Int_t PartonClassifierLHEF::GetCategory(TObject *object)
-{
-  // select parton in the parton list
-
-  Candidate *parton = static_cast<Candidate *>(object);
-  const TLorentzVector &momentum = parton->Momentum;
-  Int_t pdgCode;
-
-  // inside the eta && momentum range (be a little bit larger that the tracking coverage
-  if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
-
-  pdgCode = TMath::Abs(parton->PID);
-  if(parton->Status == -1) return -1;
-  if(pdgCode != 21 && pdgCode > 5) return -1; // not a parton, skip
-  if(parton->Status != 1) return -1; // if status 3 return
-
-  return 0;
-}
-
-//------------------------------------------------------------------------------
-
-JetFlavourAssociation::JetFlavourAssociation() :
-  fClassifier(0), fFilter(0),
-  fItPartonInputArray(0), fItPartonInputArrayLHEF(0),
-  fItJetInputArray(0), fItParticleInputArray(0)
-{
-  fClassifier    = new PartonClassifier;
-  fClassifierLHEF = new PartonClassifierLHEF;
-}
-
-//------------------------------------------------------------------------------
-
-JetFlavourAssociation::~JetFlavourAssociation()
-{
-  if(fClassifier) delete fClassifier;
-  if(fClassifierLHEF) delete fClassifierLHEF;
-}
-
-//------------------------------------------------------------------------------
-
-void JetFlavourAssociation::Init()
-{
-  ExRootConfParam param;
-
-  fDeltaR = GetDouble("DeltaR", 0.5);
-
-  fClassifier->fPTMin = GetDouble("PartonPTMin", 0.);
-  fClassifier->fEtaMax = GetDouble("PartonEtaMax",2.5);
-
-  fClassifierLHEF->fPTMin = GetDouble("PartonPTMin", 0.);
-  fClassifierLHEF->fEtaMax = GetDouble("PartonEtaMax",2.5);
-
-  // import input array(s)
-  fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
-  fItPartonInputArray = fPartonInputArray->MakeIterator();
-
-  fPartonInputArrayLHEF = ImportArray(GetString("PartonInputArrayLHEF", "Delphes/partonsLHEF"));
-  fItPartonInputArrayLHEF = fPartonInputArrayLHEF->MakeIterator();
-
-  fFilter = new ExRootFilter(fPartonInputArray);
-  fFilterLHEF = new ExRootFilter(fPartonInputArrayLHEF);
-
-  fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
-  fItJetInputArray = fJetInputArray->MakeIterator();
-
-  fParticleInputArray = ImportArray(GetString("ParticleInputArray", "Delphes/allParticles"));
-  fItParticleInputArray = fParticleInputArray->MakeIterator();
-
-}
-
-//------------------------------------------------------------------------------
-
-void JetFlavourAssociation::Finish()
-{
-  if(fFilter) delete fFilter;
-  if(fFilterLHEF) delete fFilterLHEF;
-
-  if(fItJetInputArray) delete fItJetInputArray;
-  if(fItParticleInputArray) delete fItParticleInputArray;
-  if(fItPartonInputArray) delete fItPartonInputArray;
-  if(fItPartonInputArrayLHEF) delete fItPartonInputArrayLHEF;
-}
-
-//------------------------------------------------------------------------------
-
-void JetFlavourAssociation::Process(){
-
-  Candidate *jet;
-  TObjArray *partonArray;
-  TObjArray *partonArrayLHEF;
-
-  // select quark && gluons
-  fFilter->Reset();
-  partonArray = fFilter->GetSubArray(fClassifier, 0); // get the filtered parton array
-
-  if(partonArray == 0) return;
-  TIter itPartonArray(partonArray);
-
-  fFilterLHEF->Reset();
-  partonArrayLHEF = fFilterLHEF->GetSubArray(fClassifierLHEF, 0); // get the filtered parton array
-
-  if(partonArrayLHEF == 0) return;
-  TIter itPartonLHEFArray(partonArrayLHEF);
-
-  // loop over all input jets
-  fItJetInputArray->Reset();
-  while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
-  {
-    // get standard flavor
-    GetAlgoFlavour(jet, itPartonArray, itPartonLHEFArray);
-    GetPhysicsFlavour(jet, itPartonArray, itPartonLHEFArray);
-  }
-}
-
-//------------------------------------------------------------------------------
-// Standard definition of jet flavor in
-// https://cmssdt.cern.ch/SDT/lxr/source/PhysicsTools/JetMCAlgos/plugins/JetPartonMatcher.cc?v=CMSSW_7_3_0_pre1
-
-void JetFlavourAssociation::GetAlgoFlavour(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
-{
-  float maxPt = 0;
-  float minDr = 1000;
-  bool isGoodParton = true;
-  int daughterCounter = 0;
-  Candidate *parton, *partonLHEF;
-  Candidate *tempParton = 0, *tempPartonNearest = 0, *tempPartonHighestPt = 0;
-  int pdgCode, pdgCodeMax = -1;
-
-  itPartonArray.Reset();
-  while((parton = static_cast<Candidate *>(itPartonArray.Next())))
-  {
-    // default delphes method
-    pdgCode = TMath::Abs(parton->PID);
-    if(TMath::Abs(parton->PID) == 21) pdgCode = 0;
-    if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
-    {
-      if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
-    }
-
-    isGoodParton = true;
-    itPartonLHEFArray.Reset();
-    while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
-    {
-      if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.001 &&
-         parton->PID == partonLHEF->PID &&
-         partonLHEF->Charge == parton->Charge)
-      {
-         isGoodParton = false;
-         break;
-      }
-
-      if(!isGoodParton) continue;
-
-      // check the daugheter
-      daughterCounter = 0;
-      if(parton->D1 != -1 || parton->D2 != -1)
-      {
-        // partons are only quarks || gluons
-        int daughterFlavour1 = -1;
-        int daughterFlavour2 = -1;
-        if(parton->D1 != -1) daughterFlavour1 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D1))->PID);
-        if(parton->D2 != -1) daughterFlavour2 = TMath::Abs(static_cast<Candidate *>(fParticleInputArray->At(parton->D2))->PID);
-        if((daughterFlavour1 == 1 || daughterFlavour1 == 2 || daughterFlavour1 == 3 || daughterFlavour1 == 4 || daughterFlavour1 == 5 || daughterFlavour1 == 21)) daughterCounter++;
-        if((daughterFlavour2 == 1 || daughterFlavour2 == 2 || daughterFlavour2 == 3 || daughterFlavour2 == 4 || daughterFlavour1 == 5 || daughterFlavour2 == 21)) daughterCounter++;
-      }
-      if(daughterCounter > 0) continue;
-      if(jet->Momentum.DeltaR(parton->Momentum) <= fDeltaR)
-      {
-        if(jet->Momentum.DeltaR(parton->Momentum) < minDr)
-        {
-          minDr = jet->Momentum.DeltaR(parton->Momentum);
-          tempPartonNearest = parton;
-        }
-
-        // if not yet found && pdgId is a c, take as c
-        if(TMath::Abs(parton->PID) == 4) tempParton = parton;
-        if(TMath::Abs(parton->PID) == 5) tempParton = parton;
-        if(parton->Momentum.Pt() > maxPt)
-        {
-          maxPt = parton->Momentum.Pt();
-          tempPartonHighestPt = parton;
-        }
-      }
-    }
-  }
-
-  jet->FlavorHeaviest = tempParton ? TMath::Abs(tempParton->PID) : 0;
-  jet->FlavorHighestPt = tempPartonHighestPt ? TMath::Abs(tempPartonHighestPt->PID) : 0;
-  jet->FlavorNearest2 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
-  if(!tempParton) tempParton = tempPartonHighestPt;
-  jet->FlavorAlgo = tempParton ? TMath::Abs(tempParton->PID) : 0;
-
-  if(pdgCodeMax == 0) pdgCodeMax = 21;
-  if(pdgCodeMax == -1) pdgCodeMax = 0;
-
-  jet->FlavorDefault = pdgCodeMax;
-}
-
-//------------------------------------------------------------------------------
-
-void JetFlavourAssociation::GetPhysicsFlavour(Candidate *jet, TIter &itPartonArray, TIter &itPartonLHEFArray)
-{
-  float minDr = 1000;
-  int partonCounter = 0;
-  float biggerConeSize = 0.7;
-  float dist;
-  bool isGoodCandidate;
-  int contaminatingFlavour = 0;
-  int motherCounter = 0;
-  Candidate *parton, *partonLHEF, *mother1, *mother2;
-  Candidate *tempParton = 0, *tempPartonNearest = 0;
-  vector<Candidate *> contaminations;
-  vector<Candidate *>::iterator itContaminations;
-
-  contaminations.clear();
-
-  itPartonLHEFArray.Reset();
-  while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
-  {
-    dist = jet->Momentum.DeltaR(partonLHEF->Momentum); // take the DR
-    if(partonLHEF->Status == 1 && dist < minDr)
-    {
-      tempPartonNearest = partonLHEF;
-      minDr = dist;
-    }
-
-    if(partonLHEF->Status == 1 && dist <= fDeltaR)
-    {
-      tempParton = partonLHEF;
-      partonCounter++;
-    }
-  }
-
-  itPartonArray.Reset();
-  itPartonLHEFArray.Reset();
-  while((parton = static_cast<Candidate *>(itPartonArray.Next())))
-  {
-    dist = jet->Momentum.DeltaR(parton->Momentum); // take the DR
-    isGoodCandidate = true;
-    while((partonLHEF = static_cast<Candidate *>(itPartonLHEFArray.Next())))
-    {
-      if(parton->Momentum.DeltaR(partonLHEF->Momentum) < 0.01 &&
-         parton->PID == partonLHEF->PID &&
-         partonLHEF->Charge == parton->Charge)
-      {
-        isGoodCandidate = false;
-        break;
-      }
-    }
-
-    if(!isGoodCandidate) continue;
-
-    if(parton->D1 != -1 || parton->D2 != -1)
-    {
-      if((TMath::Abs(parton->PID) < 4 || TMath::Abs(parton->PID) == 21)) continue;
-      if(dist < biggerConeSize) contaminations.push_back(parton);
-    }
-  }
-
-  jet->FlavorNearest3 = tempPartonNearest ? TMath::Abs(tempPartonNearest->PID) : 0;
-
-  if(partonCounter != 1)
-  {
-    jet->FlavorPhysics = 0;
-  }
-  else if(contaminations.size() == 0)
-  {
-    jet->FlavorPhysics = TMath::Abs(tempParton->PID);
-  }
-  else if(contaminations.size() > 0)
-  {
-    jet->FlavorPhysics = TMath::Abs(tempParton->PID);
-
-    for(itContaminations = contaminations.begin(); itContaminations != contaminations.end(); ++itContaminations)
-    {
-      parton = *itContaminations;
-      contaminatingFlavour = TMath::Abs(parton->PID);
-      motherCounter = 0;
-      if(parton->M1 != -1) motherCounter++;
-      if(parton->M2 != -1) motherCounter++;
-
-      if(parton->M1 != -1)
-      {
-        mother1 = static_cast<Candidate *>(fParticleInputArray->At(parton->M1));
-        if(mother1 && motherCounter > 0 && mother1->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
-      }
-      if(parton->M2 != -1)
-      {
-        mother2 = static_cast<Candidate *>(fParticleInputArray->At(parton->M2));
-        if(mother2 && motherCounter > 0 && mother2->Momentum.DeltaR(tempParton->Momentum) < 0.001) continue;
-      }
-      // mother is the initialParton --> OK
-      if(TMath::Abs(tempParton->PID) == 4)
-      {
-        // keep association --> the initialParton is a c --> the contaminated parton is a c
-        if(contaminatingFlavour == 4) continue;
-        jet->FlavorPhysics = 0; // all the other cases reject!
-        break;
-      }
-    }
-  }
-}
Index: dules/JetFlavourAssociation.h
===================================================================
--- modules/JetFlavourAssociation.h	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
+++ 	(revision )
@@ -1,78 +1,0 @@
-/*
- *  Delphes: a framework for fast simulation of a generic collider experiment
- *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
- *
- *  This program is free software: you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation, either version 3 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef JetFlavourAssociation_h
-#define JetFlavourAssociation_h
-
-/** \class JetFlavourAssociation
- *
- *  Find origin of jet and evaluate jet flavour
- *
- *  \author P. Demin - UCL, Louvain-la-Neuve
- *
- */
-
-#include "classes/DelphesModule.h"
-#include "classes/DelphesClasses.h"
-#include <map>
-
-class TObjArray;
-class DelphesFormula;
-
-class ExRootFilter;
-class PartonClassifier;
-class PartonClassifierLHEF;
-
-class JetFlavourAssociation: public DelphesModule
-{
-public:
-
-  JetFlavourAssociation();
-  ~JetFlavourAssociation();
-
-  void Init();
-  void Process();
-  void Finish();
-
-  void GetAlgoFlavour(Candidate *jet, TIter &itPartonArray, TIter &itPartonArrayLHEF);
-  void GetPhysicsFlavour(Candidate *jet, TIter &itPartonArray, TIter &itPartonArrayLHEF);
-
-private:
-
-  Double_t  fDeltaR;
-
-  PartonClassifier  *fClassifier; //!
-  PartonClassifierLHEF *fClassifierLHEF; //!
-
-  ExRootFilter *fFilter;
-  ExRootFilter *fFilterLHEF;
-
-  TIterator *fItPartonInputArray; //!
-  TIterator *fItPartonInputArrayLHEF; //!
-  TIterator *fItJetInputArray; //!
-  TIterator *fItParticleInputArray; //!
-
-  const TObjArray *fPartonInputArray; //!
-  const TObjArray *fPartonInputArrayLHEF; //!
-  const TObjArray *fJetInputArray; //!
-  const TObjArray *fParticleInputArray; //!
-
-  ClassDef(JetFlavourAssociation, 1)
-};
-
-#endif
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
+++ modules/ModulesLinkDef.h	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
@@ -42,4 +42,5 @@
 #include "modules/UniqueObjectFinder.h"
 #include "modules/TrackCountingBTagging.h"
+#include "modules/BTaggingCMS.h"
 #include "modules/BTagging.h"
 #include "modules/TauTagging.h"
@@ -59,5 +60,5 @@
 #include "modules/Hector.h"
 #include "modules/RunPUPPI.h"
-#include "modules/JetFlavourAssociation.h"
+#include "modules/JetFlavorAssociation.h"
 #include "modules/ExampleModule.h"
 
@@ -84,4 +85,5 @@
 #pragma link C++ class UniqueObjectFinder+;
 #pragma link C++ class TrackCountingBTagging+;
+#pragma link C++ class BTaggingCMS+;
 #pragma link C++ class BTagging+;
 #pragma link C++ class TauTagging+;
@@ -101,5 +103,5 @@
 #pragma link C++ class Hector+;
 #pragma link C++ class RunPUPPI+;
-#pragma link C++ class JetFlavourAssociation+;
+#pragma link C++ class JetFlavorAssociation+;
 #pragma link C++ class ExampleModule+;
 
Index: modules/TreeWriter.cc
===================================================================
--- modules/TreeWriter.cc	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
+++ modules/TreeWriter.cc	(revision 9040259ee50d946ebcb0b894874fcd5dfddab562)
@@ -291,5 +291,5 @@
     entry->ZOuter = position.Z();
     entry->TOuter = position.T()*1.0E-3/c_light;
-    
+
     entry->Dxy = candidate->Dxy;
     entry->SDxy = candidate->SDxy ;
@@ -297,5 +297,5 @@
     entry->Yd = candidate->Yd;
     entry->Zd = candidate->Zd;
-   
+
     const TLorentzVector &momentum = candidate->Momentum;
 
@@ -362,5 +362,5 @@
 
     entry->T = position.T()*1.0E-3/c_light;
-    entry->Ntimes = candidate->Ntimes;    
+    entry->Ntimes = candidate->Ntimes;
 
     FillParticles(candidate, &entry->Particles);
@@ -405,5 +405,5 @@
 
     // Isolation variables
-       
+
     entry->IsolationVar = candidate->IsolationVar;
     entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
@@ -453,5 +453,5 @@
 
     // Isolation variables
-       
+
     entry->IsolationVar = candidate->IsolationVar;
     entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
@@ -507,7 +507,7 @@
 
     entry->T = position.T()*1.0E-3/c_light;
-    
+
     // Isolation variables
-       
+
     entry->IsolationVar = candidate->IsolationVar;
     entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
@@ -568,4 +568,21 @@
 
     entry->BTag = candidate->BTag;
+
+    entry->BTagAlgo = candidate->BTagAlgo;
+    entry->BTagDefault = candidate->BTagDefault;
+    entry->BTagPhysics = candidate->BTagPhysics;
+    entry->BTagNearest2 = candidate->BTagNearest2;
+    entry->BTagNearest3 = candidate->BTagNearest3;
+    entry->BTagHeaviest = candidate->BTagHeaviest;
+    entry->BTagHighestPt = candidate->BTagHighestPt;
+
+    entry->FlavorAlgo = candidate->FlavorAlgo;
+    entry->FlavorDefault = candidate->FlavorDefault;
+    entry->FlavorPhysics = candidate->FlavorPhysics;
+    entry->FlavorNearest2 = candidate->FlavorNearest2;
+    entry->FlavorNearest3 = candidate->FlavorNearest3;
+    entry->FlavorHeaviest = candidate->FlavorHeaviest;
+    entry->FlavorHighestPt = candidate->FlavorHighestPt;
+
     entry->TauTag = candidate->TauTag;
 
@@ -593,20 +610,20 @@
     entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
     entry->PTD = candidate->PTD;
-    
+
     //--- Sub-structure variables ----
-    
-    entry->NSubJetsTrimmed     = candidate->NSubJetsTrimmed;
-    entry->NSubJetsPruned      = candidate->NSubJetsPruned; 
+
+    entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed;
+    entry->NSubJetsPruned = candidate->NSubJetsPruned;
     entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped;
-    
+
     for(i = 0; i < 5; i++)
     {
-      entry->FracPt[i]        = candidate -> FracPt[i];
-      entry->Tau[i]           = candidate -> Tau[i];
-      entry->TrimmedP4[i]     = candidate -> TrimmedP4[i];
-      entry->PrunedP4[i]      = candidate -> PrunedP4[i];
+      entry->FracPt[i] = candidate -> FracPt[i];
+      entry->Tau[i] = candidate -> Tau[i];
+      entry->TrimmedP4[i] = candidate -> TrimmedP4[i];
+      entry->PrunedP4[i] = candidate -> PrunedP4[i];
       entry->SoftDroppedP4[i] = candidate -> SoftDroppedP4[i];
     }
-   
+
     FillParticles(candidate, &entry->Particles);
   }
