Index: Makefile
===================================================================
--- Makefile	(revision 8707eebb4a19f7c7f1aa2ced879f327d1bd85a2e)
+++ Makefile	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
@@ -339,4 +339,5 @@
 	modules/Hector.h \
 	modules/RunPUPPI.h \
+	modules/JetFlavourAssociation.h \
 	modules/ExampleModule.h
 ModulesDict$(PcmSuf): \
@@ -657,4 +658,13 @@
 	external/ExRootAnalysis/ExRootFilter.h \
 	external/ExRootAnalysis/ExRootClassifier.h
+tmp/modules/JetFlavourAssociation.$(ObjSuf): \
+	modules/JetFlavourAssociation.$(SrcSuf) \
+	modules/JetFlavourAssociation.h \
+	classes/DelphesClasses.h \
+	classes/DelphesFactory.h \
+	classes/DelphesFormula.h \
+	external/ExRootAnalysis/ExRootResult.h \
+	external/ExRootAnalysis/ExRootFilter.h \
+	external/ExRootAnalysis/ExRootClassifier.h
 tmp/modules/JetPileUpSubtractor.$(ObjSuf): \
 	modules/JetPileUpSubtractor.$(SrcSuf) \
@@ -899,4 +909,5 @@
 	tmp/modules/ImpactParameterSmearing.$(ObjSuf) \
 	tmp/modules/Isolation.$(ObjSuf) \
+	tmp/modules/JetFlavourAssociation.$(ObjSuf) \
 	tmp/modules/JetPileUpSubtractor.$(ObjSuf) \
 	tmp/modules/LeptonDressing.$(ObjSuf) \
@@ -1751,4 +1762,9 @@
 	@touch $@
 
+modules/JetFlavourAssociation.h: \
+	classes/DelphesModule.h \
+	classes/DelphesClasses.h
+	@touch $@
+
 modules/TrackCountingBTagging.h: \
 	classes/DelphesModule.h
Index: classes/DelphesClasses.cc
===================================================================
--- classes/DelphesClasses.cc	(revision 8707eebb4a19f7c7f1aa2ced879f327d1bd85a2e)
+++ classes/DelphesClasses.cc	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
@@ -121,5 +121,7 @@
   Charge(0), Mass(0.0),
   IsPU(0), IsRecoPU(0), IsConstituent(0),
-  BTag(0), TauTag(0), Eem(0.0), Ehad(0.0),
+  BTag(0), BTagAlgo(0), BTagDefault(0), BTagPhysics(0), BTagNearest2(0), BTagNearest3(0), BTagHeaviest(0), BTagHighestPt(0),
+  FlavorAlgo(0), FlavorDefault(0), FlavorPhysics(0), FlavorNearest2(0), FlavorNearest3(0), FlavorHeaviest(0), FlavorHighestPt(0),
+  TauTag(0), Eem(0.0), Ehad(0.0),
   DeltaEta(0.0), DeltaPhi(0.0),
   Momentum(0.0, 0.0, 0.0, 0.0),
@@ -242,4 +244,18 @@
   object.IsConstituent = IsConstituent;
   object.BTag = BTag;
+  object.BTagAlgo = BTagAlgo;
+  object.BTagDefault = BTagDefault;
+  object.BTagPhysics = BTagPhysics;
+  object.BTagNearest2 = BTagNearest2;
+  object.BTagNearest3 = BTagNearest3;
+  object.BTagHeaviest = BTagHeaviest;
+  object.BTagHighestPt = BTagHighestPt;
+  object.FlavorAlgo = FlavorAlgo;
+  object.FlavorDefault = FlavorDefault;
+  object.FlavorPhysics = FlavorPhysics;
+  object.FlavorNearest2 = FlavorNearest2;
+  object.FlavorNearest3 = FlavorNearest3;
+  object.FlavorHeaviest = FlavorHeaviest;
+  object.FlavorHighestPt = FlavorHighestPt;
   object.TauTag = TauTag;
   object.Eem = Eem;
@@ -301,6 +317,6 @@
   object.SoftDroppedP4[4] = SoftDroppedP4[4];
 
-  object.NSubJetsTrimmed     = NSubJetsTrimmed;
-  object.NSubJetsPruned      = NSubJetsPruned;
+  object.NSubJetsTrimmed = NSubJetsTrimmed;
+  object.NSubJetsPruned = NSubJetsPruned;
   object.NSubJetsSoftDropped = NSubJetsSoftDropped;
 
@@ -337,4 +353,18 @@
   IsConstituent = 0;
   BTag = 0;
+  BTagAlgo = 0;
+  BTagDefault = 0;
+  BTagPhysics = 0;
+  BTagNearest2 = 0;
+  BTagNearest3 = 0;
+  BTagHeaviest = 0;
+  BTagHighestPt = 0;
+  FlavorAlgo = 0;
+  FlavorDefault = 0;
+  FlavorPhysics = 0;
+  FlavorNearest2 = 0;
+  FlavorNearest3 = 0;
+  FlavorHeaviest = 0;
+  FlavorHighestPt = 0;
   TauTag = 0;
   Eem = 0.0;
@@ -389,6 +419,6 @@
   }
 
-  NSubJetsTrimmed     = 0;
-  NSubJetsPruned      = 0;
+  NSubJetsTrimmed = 0;
+  NSubJetsPruned = 0;
   NSubJetsSoftDropped = 0;
   
Index: classes/DelphesClasses.h
===================================================================
--- classes/DelphesClasses.h	(revision 8707eebb4a19f7c7f1aa2ced879f327d1bd85a2e)
+++ classes/DelphesClasses.h	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
@@ -243,5 +243,5 @@
 
   // Isolation variables
-  
+
   Float_t IsolationVar;
   Float_t IsolationVarRhoCorr;
@@ -278,5 +278,5 @@
 
   // Isolation variables
-  
+
   Float_t IsolationVar;
   Float_t IsolationVarRhoCorr;
@@ -311,5 +311,5 @@
 
    // Isolation variables
-  
+
   Float_t IsolationVar;
   Float_t IsolationVarRhoCorr;
@@ -345,4 +345,21 @@
 
   UInt_t BTag; // 0 or 1 for a jet that has been tagged as containing a heavy quark
+
+  UInt_t BTagAlgo;
+  UInt_t BTagDefault;
+  UInt_t BTagPhysics;
+  UInt_t BTagNearest2;
+  UInt_t BTagNearest3;
+  UInt_t BTagHeaviest;
+  UInt_t BTagHighestPt;
+
+  UInt_t FlavorAlgo;
+  UInt_t FlavorDefault;
+  UInt_t FlavorPhysics;
+  UInt_t FlavorNearest2;
+  UInt_t FlavorNearest3;
+  UInt_t FlavorHeaviest;
+  UInt_t FlavorHighestPt;
+
   UInt_t TauTag; // 0 or 1 for a jet that has been tagged as a tau
 
@@ -351,22 +368,22 @@
   Float_t EhadOverEem; // ratio of the hadronic versus electromagnetic energy deposited in the calorimeter
 
-  Int_t    NCharged; // number of charged constituents 
-  Int_t    NNeutrals; // number of neutral constituents 
-  Float_t  Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents) 
-  Float_t  BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents) 
-  Float_t  MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent
-  Float_t  PTD; // average pt between constituent and jet weighted by pt of constituent
-  Float_t  FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents) 
+  Int_t NCharged; // number of charged constituents
+  Int_t NNeutrals; // number of neutral constituents
+  Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents)
+  Float_t BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents)
+  Float_t MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent
+  Float_t PTD; // average pt between constituent and jet weighted by pt of constituent
+  Float_t FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents)
 
   Float_t Tau[5]; // N-subjettiness
-  
-  TLorentzVector TrimmedP4[5]; //first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta
-  TLorentzVector PrunedP4[5]; //first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
-  TLorentzVector SoftDroppedP4[5]; //first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
-  
-  Int_t    NSubJetsTrimmed; // number of subjets trimmed 
-  Int_t    NSubJetsPruned; // number of subjets pruned 
-  Int_t    NSubJetsSoftDropped; // number of subjets soft-dropped 
-   
+
+  TLorentzVector TrimmedP4[5]; // first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta
+  TLorentzVector PrunedP4[5]; // first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
+  TLorentzVector SoftDroppedP4[5]; // first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
+
+  Int_t NSubJetsTrimmed; // number of subjets trimmed
+  Int_t NSubJetsPruned; // number of subjets pruned
+  Int_t NSubJetsSoftDropped; // number of subjets soft-dropped
+
   TRefArray Constituents; // references to constituents
   TRefArray Particles; // references to generated particles
@@ -437,5 +454,5 @@
   Float_t T; // ecal deposit time, averaged by sqrt(EM energy) over all particles, not smeared
   Int_t   Ntimes; // number of hits contributing to time measurement
-  
+
   Float_t Eem; // calorimeter tower electromagnetic energy
   Float_t Ehad; // calorimeter tower hadronic energy
@@ -497,8 +514,25 @@
   Int_t IsPU;
   Int_t IsRecoPU;
-  
+
   Int_t IsConstituent;
 
   UInt_t BTag;
+
+  UInt_t BTagAlgo;
+  UInt_t BTagDefault;
+  UInt_t BTagPhysics;
+  UInt_t BTagNearest2;
+  UInt_t BTagNearest3;
+  UInt_t BTagHeaviest;
+  UInt_t BTagHighestPt;
+
+  UInt_t FlavorAlgo;
+  UInt_t FlavorDefault;
+  UInt_t FlavorPhysics;
+  UInt_t FlavorNearest2;
+  UInt_t FlavorNearest3;
+  UInt_t FlavorHeaviest;
+  UInt_t FlavorHighestPt;
+
   UInt_t TauTag;
 
@@ -527,12 +561,12 @@
   Float_t  PTD;
   Float_t  FracPt[5];
-  
+
   //Timing information
-  
+
   Int_t    Ntimes;
-  std::vector<std::pair<Float_t,Float_t> > Ecal_E_t; 
+  std::vector<std::pair<Float_t,Float_t> > Ecal_E_t;
 
   // Isolation variables
-  
+
   Float_t IsolationVar;
   Float_t IsolationVarRhoCorr;
@@ -545,14 +579,14 @@
 
   Float_t Tau[5];
-  
+
   // Other Substructure variables
-  
-  TLorentzVector TrimmedP4[5]; //first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta
-  TLorentzVector PrunedP4[5]; //first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
-  TLorentzVector SoftDroppedP4[5]; //first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
- 
-  Int_t    NSubJetsTrimmed; // number of subjets trimmed 
-  Int_t    NSubJetsPruned; // number of subjets pruned 
-  Int_t    NSubJetsSoftDropped; // number of subjets soft-dropped 
+
+  TLorentzVector TrimmedP4[5]; // first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta
+  TLorentzVector PrunedP4[5]; // first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
+  TLorentzVector SoftDroppedP4[5]; // first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta
+
+  Int_t NSubJetsTrimmed; // number of subjets trimmed
+  Int_t NSubJetsPruned; // number of subjets pruned
+  Int_t NSubJetsSoftDropped; // number of subjets soft-dropped
 
 
Index: modules/JetFlavourAssociation.cc
===================================================================
--- modules/JetFlavourAssociation.cc	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
+++ modules/JetFlavourAssociation.cc	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
@@ -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 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: modules/JetFlavourAssociation.h
===================================================================
--- modules/JetFlavourAssociation.h	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
+++ modules/JetFlavourAssociation.h	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
@@ -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 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 8707eebb4a19f7c7f1aa2ced879f327d1bd85a2e)
+++ modules/ModulesLinkDef.h	(revision edf10ba9b647bfb1cdbdfed3f4ef281a28c3ee16)
@@ -59,4 +59,5 @@
 #include "modules/Hector.h"
 #include "modules/RunPUPPI.h"
+#include "modules/JetFlavourAssociation.h"
 #include "modules/ExampleModule.h"
 
@@ -100,4 +101,5 @@
 #pragma link C++ class Hector+;
 #pragma link C++ class RunPUPPI+;
+#pragma link C++ class JetFlavourAssociation+;
 #pragma link C++ class ExampleModule+;
 
