Merge lp:~ma5dev/madanalysis5/v1.5beta into lp:madanalysis5

[Eric Conte]

=== modified file 'bin/ma5'
--- bin/ma5	2016-12-12 11:48:45 +0000
+++ bin/ma5	2016-12-20 15:49:11 +0000
@@ -62,8 +62,13 @@
 
 # Release version
 # Do not touch it !!!!!  
+<<<<<<< TREE
 version = "1.5"
 date = "2016/12/12"
+=======
+version = "1.5.33"
+date = "2016/12/20"
+>>>>>>> MERGE-SOURCE
 
 # Loading the MadAnalysis session
 import madanalysis.core.launcher

=== modified file 'madanalysis/UpdateNotes.txt'
--- madanalysis/UpdateNotes.txt	2016-12-12 11:48:45 +0000
+++ madanalysis/UpdateNotes.txt	2016-12-20 15:49:11 +0000
@@ -1,5 +1,6 @@
 Update notes for MadAnalysis 5 (in reverse time order)
 
+<<<<<<< TREE
 159  1.5   (2016/11/14) ma5team:    - embedded into MadGraph
                                     - FIFO mode for MadGraph
                                     - No more RecastingTools
@@ -10,6 +11,75 @@
                                     - protections added for PADMA5tune and delphesMA5tune -> not 
                                       compatible with root 6 
                                     - extend the list of recasted analyses
+=======
+172           (2016/12/20) econte:  y 
+
+171           (2016/12/20) econte:  restart in script mode 
+
+170           (2016/12/19) econte:  test 
+
+169           (2016/12/19) econte:  test 
+
+168           (2016/12/19) econte:  test 
+
+167           (2016/12/19) econte:  test 
+
+166           (2016/12/19) econte:  test 
+
+165           (2016/12/19) econte:  removed unuseless files 
+
+159  1.5.25   (2016/11/14) fuks:    Bug fix in the madgraph interface (multiparticles) 
+
+158  1.5.24   (2016/10/28) fuks:    Bug in the eflow isolation 
+
+157  1.5.23   (2016/10/27) fuks:    Fixing root detection issue 
+
+153  1.5.22   (2016/09/21) fuks:    Fixing matplotlib installation 
+
+150  1.5.21   (2016/09/03) fuks:    Fixing the mg5-ma5 interface 
+
+149  1.5.20   (2016/09/03) fuks:    fixing last few things for the mg5-ma5 interface 
+
+148  1.5.19   (2016/09/02) fuks:    ma5-mg5 interface 
+
+147  1.5.18   (2016/09/01) fuks:    small fixes for the mg5-ma5 interface 
+
+146  1.5.17   (2016/09/01) fuks:    Fixing fastjet plots for the mg5-ma5 interface 
+
+145  1.5.16   (2016/09/01) fuks:    ma5-mg5 interface 
+
+144  1.5.15   (2016/09/01) fuks:    Fixing installation options for the mg5-ma5 interface 
+
+142  1.5.14   (2016/09/01) fuks:    Protections added for PADMA5tune and delphesMA5tune -> not 
+                                   compatible with root 6 
+
+141  1.5.13   (2016/09/01) fuks:    gzip fix 
+
+140  1.5.12   (2016/09/01) fuks:    FIFO step 1 
+
+137  1.5.11   (2016/08/30) fuks:    Few bugs in the detection module 
+
+133  1.5.10   (2016/08/30) fuks:    Fixing last bugs for the ma5-mg5 interface 
+
+131  1.5.9    (2016/08/30) fuks:    Fixing ma5-mg5 interface issues 
+
+127  1.5.8    (2016/08/30) fuks:    Fixing issues related to the mg5_amc interface 
+
+121  1.5.7    (2016/08/26) fuks:    fixing a few bugs 
+
+120  1.5.6    (2016/08/26) fuks:    Fixing a few bugs 
+
+119  1.5.5    (2016/08/24) fuks:    Small bug with the colors 
+
+118  1.5.4    (2016/07/24) fuks:    Fixing bugs in the recasting module when negative weights are around 
+
+117  1.5.3    (2016/07/24) fuks:    Fixing a few bug in the recasting module (reading of the recasting 
+                                   card) 
+
+115  1.6      (2016/07/23) fuks:    No more RecastingTools 
+
+114  1.5      (2016/07/23) fuks:    Fixing the root detection with the delphes activation methods 
+>>>>>>> MERGE-SOURCE
 
 111  1.4    (2016/07/20) ma5team:   - Root is now optional and Pyroot is not a requirement anymore
                                     - matplotlib can be used for the histograms

=== modified file 'madanalysis/interpreter/interpreter.py'
--- madanalysis/interpreter/interpreter.py	2016-12-08 10:28:55 +0000
+++ madanalysis/interpreter/interpreter.py	2016-12-20 15:49:11 +0000
@@ -200,29 +200,35 @@
      # Restart
     def do_restart(self, line):
         """ sending a signal allowing to restart the interpreter """
+
+        # Restart not available in script mode
+        # avoiding that the script is read from the beginning
         if self.main.script:
             logging.getLogger('MA5').warning("'restart' command is not allowed in script mode.")
+            return None
+
+        # Asking the safety question
+        YES=True
+        if not Main.forced:
+           logging.getLogger('MA5').warning("Are you sure to restart the MadAnalysis 5 session? (Y/N)")
+           allowed_answers=['n','no','y','yes']
+           answer=""
+           while answer not in  allowed_answers:
+              answer=raw_input("Answer: ")
+              answer=answer.lower()
+              if answer=="no" or answer=="n":
+                   YES=False
+                   break
+              elif answer=='yes' or answer=='y':
+                   YES=True
+                   break
+
+        # Restart?
+        if YES:
+            self.main.repeatSession=True
+            return True
         else:
-            YES=False
-            # Asking the safety question
-            if not Main.forced:
-                logging.getLogger('MA5').warning("Are you sure to restart the MadAnalysis 5 session? (Y/N)")
-                allowed_answers=['n','no','y','yes']
-                answer=""
-                while answer not in  allowed_answers:
-                   answer=raw_input("Answer: ")
-                   answer=answer.lower()
-                   if answer=="no" or answer=="n":
-                       YES=False
-                       break
-                   elif answer=='yes' or answer=='y':
-                       YES=True
-                       break
-            if YES:
-                self.main.repeatSession=True
-                return True
-            else:
-                pass
+            return None
 
     def help_restart(self):
         logging.getLogger('MA5').info("   Syntax: restart ")

=== modified file 'tools/SampleAnalyzer/Commons/Base/Configuration.cpp'
--- tools/SampleAnalyzer/Commons/Base/Configuration.cpp	2016-11-14 10:06:38 +0000
+++ tools/SampleAnalyzer/Commons/Base/Configuration.cpp	2016-12-20 15:49:11 +0000
@@ -40,8 +40,8 @@
 // Initializing static data members
 // -----------------------------------------------------------------------------
 // DO NOT TOUCH THESE LINES
-const std::string Configuration::sampleanalyzer_version_ = "1.5.25";
-const std::string Configuration::sampleanalyzer_date_    = "2016/11/14";
+const std::string Configuration::sampleanalyzer_version_ = "1.5.33";
+const std::string Configuration::sampleanalyzer_date_    = "2016/12/20";
 // DO NOT TOUCH THESE LINES
 
 // -----------------------------------------------------------------------------

=== removed file 'tools/SampleAnalyzer/Interfaces/fastjet/JetClusteringFastJet.bak'
--- tools/SampleAnalyzer/Interfaces/fastjet/JetClusteringFastJet.bak	2014-05-30 16:41:57 +0000
+++ tools/SampleAnalyzer/Interfaces/fastjet/JetClusteringFastJet.bak	1970-01-01 00:00:00 +0000
@@ -1,345 +0,0 @@
-////////////////////////////////////////////////////////////////////////////////
-//  
-//  Copyright (C) 2012-2013 Eric Conte, Benjamin Fuks
-//  The MadAnalysis development team, email: <[e-mail]>
-//  
-//  This file is part of MadAnalysis 5.
-//  Official website: <https://launchpad.net/madanalysis5>
-//  
-//  MadAnalysis 5 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.
-//  
-//  MadAnalysis 5 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 MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
-//  
-////////////////////////////////////////////////////////////////////////////////
-
-
-//SampleAnalyser headers
-#include "SampleAnalyzer/Interfaces/fastjet/JetClusteringFastJet.h"
-#include "SampleAnalyzer/Commons/Service/LoopService.h"
-
-//FastJet headers
-#include <fastjet/ClusterSequence.hh>
-#include <fastjet/PseudoJet.hh>
-
-using namespace MA5;
-
-
-JetClusteringFastJet::JetClusteringFastJet(std::string Algo):JetClusterAlgoBase(Algo)
-{ JetAlgorithm_=Algo; JetDefinition_=0; }
-
-JetClusteringFastJet::~JetClusteringFastJet() 
-{ if (JetDefinition_!=0) delete JetDefinition_; }
-
-
-void JetClusteringFastJet::GetFinalState(const MCParticleFormat* part, std::set<const MCParticleFormat*>& finalstates)
-{
-  for (unsigned int i=0; i<part->daughters().size(); i++)
-  {
-    if (PHYSICS->Id->IsFinalState(part->daughters()[i])) finalstates.insert(part->daughters()[i]);
-    else return GetFinalState(part->daughters()[i],finalstates);
-  }
-}
-
-Bool_t JetClusteringFastJet::IsLast(const MCParticleFormat* part, EventFormat& myEvent)
-{
-  for (unsigned int i=0; i<part->daughters().size(); i++)
-  {
-    if (part->daughters()[i]->pdgid()==part->pdgid()) return false;
-  }
-  return true;
-}
-
-
-bool JetClusteringFastJet::Execute(SampleFormat& mySample, EventFormat& myEvent)
-{
-  // Veto
-  std::vector<bool> vetos(myEvent.mc()->particles().size(),false);
-  std::set<const MCParticleFormat*> vetos2;
-
-  // Filling the dataformat with electron/muon
-  for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
-  {
-    const MCParticleFormat& part = myEvent.mc()->particles()[i];
-    UInt_t absid = std::abs(part.pdgid());
-
-    // Rejecting particle with a null pt (initial state ?)
-    if (part.pt()<1e-10) continue;
-
-    // Treating intermediate particles
-    if (PHYSICS->Id->IsInterState(part))
-    {
-      // rejecting not interesting particles
-      if (absid!=5 && absid!=4 && absid!=15) continue;
-
-      // keeping the last particle with the same id in the decay chain
-      if (!IsLast(&part, myEvent)) continue;
-
-      // looking for b quarks
-      if (absid==5)
-      {
-        bool found=false;
-        for (unsigned int j=0;j<myEvent.rec()->MCBquarks_.size();j++)
-        {
-          if (myEvent.rec()->MCBquarks_[j]==&(part)) 
-          {found=true; break;}
-        }
-        if (!found) myEvent.rec()->MCBquarks_.push_back(&(part));
-      }
-
-      // looking for c quarks
-      else if (absid==4)
-      {
-        bool found=false;
-        for (unsigned int j=0;j<myEvent.rec()->MCCquarks_.size();j++)
-        {
-          if (myEvent.rec()->MCCquarks_[j]==&(part)) 
-          {found=true; break;}
-        }
-        if (!found) myEvent.rec()->MCCquarks_.push_back(&(part));
-      }
-
-      // looking for taus
-      else if (absid==15)
-      {
-        // rejecting particle if coming from hadronization
-        if (LOOP->ComingFromHadronDecay(&part,mySample)) continue;
-
-        // Looking taus daughters id
-        bool leptonic   = true;
-        bool muonic     = false;
-        bool electronic = false;
-        for (unsigned int j=0;j<part.daughters().size();j++)
-        {
-          UInt_t pdgid = std::abs(part.daughters()[j]->pdgid());
-          if      (pdgid==13) muonic=true;
-          else if (pdgid==11) electronic=true;
-          else if (pdgid!=22 /*photons*/ &&
-                   !(pdgid>=11 && pdgid<=16) /*neutrinos*/) 
-                  leptonic=false;
-        }
-        if (!leptonic) {muonic=false; electronic=false;}
-
-        // Saving taus decaying into muons (only one copy)
-        if (muonic)
-        {
-          bool found=false;
-          for (unsigned int j=0;j<myEvent.rec()->MCMuonicTaus_.size();j++)
-          {
-            if (myEvent.rec()->MCMuonicTaus_[j]==&(part)) 
-            {found=true; break;}
-          }
-          if (!found) myEvent.rec()->MCMuonicTaus_.push_back(&(part));
-        }
-
-        // Saving taus decaying into electrons (only one copy)
-        else if (electronic)
-        {
-          bool found=false;
-          for (unsigned int j=0;j<myEvent.rec()->MCElectronicTaus_.size();j++)
-          {
-            if (myEvent.rec()->MCElectronicTaus_[j]==&(part)) 
-            {found=true; break;}
-          }
-          if (!found) myEvent.rec()->MCElectronicTaus_.push_back(&(part));
-        }
-
-        // Saving taus decaying into hadrons (only copy)
-        else
-        {
-          bool found=false;
-          for (unsigned int j=0;j<myEvent.rec()->MCHadronicTaus_.size();j++)
-          {
-            if (myEvent.rec()->MCHadronicTaus_[j]==&(part)) 
-            {found=true; break;}
-          }
-          if (!found) 
-          {
-            // Saving the hadrons in MC container
-            myEvent.rec()->MCHadronicTaus_.push_back(&(part));
-
-            // Applying efficiency
-            if (!myTAUtagger_->IsIdentified()) continue;
-
-            // Creating reco hadronic taus
-            RecTauFormat* myTau = myEvent.rec()->GetNewTau();
-            if (part.pdgid()>0) myTau->setCharge(-1);
-            else myTau->setCharge(+1);
-            myTau->setMomentum(part.momentum());
-            myTau->setMc(&part);
-            myTau->setDecayMode(PHYSICS->GetTauDecayMode(myTau->mc()));
-            if (myTau->DecayMode()<=0) myTau->setNtracks(0); // ERROR case
-            else if (myTau->DecayMode()==7 || 
-                     myTau->DecayMode()==9) myTau->setNtracks(3); // 3-Prong
-            else myTau->setNtracks(1); // 1-Prong
-
-            // Searching final state
-            GetFinalState(&part,vetos2);
-          }
-        }
-      }
-    }
-
-    // Keeping only final states
-    else if (PHYSICS->Id->IsFinalState(part))
-    {
-      // keeping only electron, muon and photon
-      if (absid!=22 && absid!=11 && absid!=13) continue;
-
-      // rejecting particle if coming from hadronization
-      if (ExclusiveId_ && LOOP->ComingFromHadronDecay(&part,mySample)) continue;
-
-      // Muons
-      if (absid==13)
-      {
-        vetos[i]=true;
-        RecLeptonFormat * muon = myEvent.rec()->GetNewMuon();
-        muon->setMomentum(part.momentum());
-        muon->setMc(&(part));
-        if (part.pdgid()==13) muon->SetCharge(-1);
-        else muon->SetCharge(+1);
-      }
-
-      // Electrons
-      else if (absid==11)
-      {
-        vetos[i]=true;
-        RecLeptonFormat * elec = myEvent.rec()->GetNewElectron();
-        elec->setMomentum(part.momentum());
-        elec->setMc(&(part));
-        if (part.pdgid()==11) elec->SetCharge(-1);
-        else elec->SetCharge(+1);
-      }
-
-      // Photons
-      else if (absid==22)
-      {
-        if (LOOP->IrrelevantPhoton(&part,mySample)) continue;
-        vetos[i]=true;
-        RecPhotonFormat * photon = myEvent.rec()->GetNewPhoton();
-        photon->setMomentum(part.momentum());
-        photon->setMc(&(part));
-      }
-    }
-  }
-
-  double & TET = myEvent.rec()->TET();
-  double & THT = myEvent.rec()->THT();
-
-  // Preparing inputs
-  std::vector<fastjet::PseudoJet> inputs;
-  for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
-  {
-    const MCParticleFormat& part = myEvent.mc()->particles()[i];
-
-    // Selecting input for jet clustering
-    if (myEvent.mc()->particles()[i].statuscode()!=1)       continue;
-    if (PHYSICS->Id->IsInvisible(myEvent.mc()->particles()[i])) continue;
-
-    // ExclusiveId mode
-    if (ExclusiveId_)
-    {
-      if (vetos[i]) continue;
-      if (vetos2.find(&part)!=vetos2.end()) continue;
-    }
-
-    // NonExclusive Id mode
-    else if (std::abs(myEvent.mc()->particles()[i].pdgid())==13) continue;
-
-    // Filling good particle for clustering
-    inputs.push_back(
-          fastjet::PseudoJet ( myEvent.mc()->particles()[i].px(), 
-                               myEvent.mc()->particles()[i].py(),
-                               myEvent.mc()->particles()[i].pz(),
-                               myEvent.mc()->particles()[i].e()   ));
-    inputs.back().set_user_index(i);
-  }
-
-  // Clustering
-  fastjet::ClusterSequence clust_seq(inputs, *JetDefinition_);
-
-  // Getting jets with PTmin = 0
-  std::vector<fastjet::PseudoJet> jets; 
-  if (Exclusive_) jets = clust_seq.exclusive_jets(0.);
-  else jets = clust_seq.inclusive_jets(0.);
-
-  // Calculating the MET  
-  ParticleBaseFormat* MET = myEvent.rec()->GetNewMet();
-  ParticleBaseFormat* MHT = myEvent.rec()->GetNewMht();
-
-  for (unsigned int i=0;i<jets.size();i++)
-  {
-    TLorentzVector q(jets[i].px(),jets[i].py(),jets[i].pz(),jets[i].e());
-    (*MET) -= q;
-    (*MHT) -= q;
-    THT += jets[i].pt();
-    TET += jets[i].pt();
-  }
-
-  // Getting jets with PTmin
-  if (Exclusive_) jets = clust_seq.exclusive_jets(Ptmin_);
-  else jets = clust_seq.inclusive_jets(Ptmin_);
-
-  // Filling the dataformat with jets
-  for (unsigned int i=0;i<jets.size();i++)
-  {
-    RecJetFormat * jet = myEvent.rec()->GetNewJet();
-    jet->setMomentum(TLorentzVector(jets[i].px(),jets[i].py(),jets[i].pz(),jets[i].e()));
-    std::vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[i]);
-    UInt_t tracks = 0;
-    for (unsigned int j=0;j<constituents.size();j++)
-    {
-      jet->AddConstituent(constituents[j].user_index());
-      //      if (std::abs(myEvent.mc()->particles()[constituents[j].user_index()].pdgid())==11) continue;
-      if (PDG->IsCharged(myEvent.mc()->particles()[constituents[j].user_index()].pdgid())) tracks++;
-    }
-    jet->ntracks_ = tracks;
-  }
-
-
-  if (ExclusiveId_)
-  {
-    for (unsigned int i=0;i<myEvent.rec()->electrons().size();i++)
-    {
-      (*MET) -= myEvent.rec()->electrons()[i].momentum();
-      TET += myEvent.rec()->electrons()[i].pt();
-    }
-    for (unsigned int i=0;i<myEvent.rec()->photons().size();i++)
-    {
-      (*MET) -= myEvent.rec()->photons()[i].momentum();
-      TET += myEvent.rec()->photons()[i].pt();
-    }
-    for (unsigned int i=0;i<myEvent.rec()->taus().size();i++)
-    {
-      (*MET) -= myEvent.rec()->taus()[i].momentum();
-      TET += myEvent.rec()->taus()[i].pt();
-    }
-  }
-
-  for (unsigned int i=0;i<myEvent.rec()->muons().size();i++)
-  {
-    (*MET) -= myEvent.rec()->muons()[i].momentum();
-    TET += myEvent.rec()->muons()[i].pt();
-  }
-
-  MET->momentum().SetPz(0.);
-  MET->momentum().SetE(MET->momentum().Pt());
-  MHT->momentum().SetPz(0.);
-  MHT->momentum().SetE(MHT->momentum().Pt());
-
-  myBtagger_->Execute(mySample,myEvent);
-  myTAUtagger_->Execute(mySample,myEvent);
-
-
-  return true;
-}
-
-

=== removed file 'tools/SampleAnalyzer/Process/Analyzer/MergingPlots.bak'
--- tools/SampleAnalyzer/Process/Analyzer/MergingPlots.bak	2014-05-30 09:59:33 +0000
+++ tools/SampleAnalyzer/Process/Analyzer/MergingPlots.bak	1970-01-01 00:00:00 +0000
@@ -1,484 +0,0 @@
-////////////////////////////////////////////////////////////////////////////////
-//  
-//  Copyright (C) 2012-2013 Eric Conte, Benjamin Fuks
-//  The MadAnalysis development team, email: <[e-mail]>
-//  
-//  This file is part of MadAnalysis 5.
-//  Official website: <https://launchpad.net/madanalysis5>
-//  
-//  MadAnalysis 5 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.
-//  
-//  MadAnalysis 5 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 MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
-//  
-////////////////////////////////////////////////////////////////////////////////
-
-
-//SampleAnalyzer headers
-#include "SampleAnalyzer/Process/Analyzer/MergingPlots.h"
-#include "SampleAnalyzer/Commons/Base/Configuration.h"
-#include "SampleAnalyzer/Commons/Service/Physics.h"
-#include "SampleAnalyzer/Commons/Service/CompilationService.h"
-
-
-//STL headers
-#include <sstream>
-
-
-using namespace MA5;
-
-/*
-extern"C"
-{
-  void ktclus_(int *imode, double PP[512][4], int* NN, double* ECUT, double Y[512]);
-}
-*/
-
-
-bool MergingPlots::Initialize(const Configuration& cfg,
-             const std::map<std::string,std::string>& parameters)
-{
-  // Options
-  merging_njets_=0;
-  merging_nqmatch_=4;
-  merging_nosingrad_=false;
-
-  // Reading options
-  for (std::map<std::string,std::string>::const_iterator it=parameters.begin();
-       it!=parameters.end();it++)
-  {
-    if (it->first=="njets")
-    {
-      std::stringstream str;
-      str << it->second;
-      str >> merging_njets_;
-    }
-    else
-    {
-      WARNING << "parameter '" << it->first 
-              << "' is unknown and will be ignored." << endmsg;
-    }
-  }
-
-  // Initializing DJR plots
-  if (merging_njets_==0) 
-  {
-    ERROR << "number of jets requested for DJR plots is zero" << endmsg;
-    return false;
-  }
-  DJR_.resize(merging_njets_);
-  for (unsigned int i=0;i<DJR_.size();i++)
-  {
-    std::stringstream str;
-    str << "DJR" << i+1;
-    std::string title;
-    str >> title;
-    DJR_[i].Initialize(DJR_.size()+1,title);
-  }
-
-  // Initializing clustering algorithm
-  JetDefinition_ = new fastjet::JetDefinition(fastjet::kt_algorithm,1.0);
-  return true;
-}
-
-
-void MergingPlots::Execute(SampleFormat& mySample, const EventFormat& myEvent)
-{
-  // Safety
-  if (mySample.mc()==0) return;
-  if (myEvent.mc()==0) return;
-
-  // Getting number of additional jets in the event
-  UInt_t njets = ExtractJetNumber(myEvent.mc(),mySample.mc());
-  if (njets>DJR_.size()) return;
-
-  // Preparing inputs
-  std::vector<fastjet::PseudoJet> inputs;
-  SelectParticles(inputs,myEvent.mc());
-
-  // Getting DJR observables
-  std::vector<Double_t> DJRvalues(DJR_.size(),0.);
-  ExtractDJR(inputs,DJRvalues);
-
-
-  // Getting results
-  for (unsigned int i=0;i<DJR_.size();i++)
-  {
-    /*     {    double djr = std::log10(sqrt(DJRvalues[i]));
-            if (i==1 && njets==1 && djr > 1.8)
-      {
-        for (unsigned int j=0;j<myEvent->particles().size();j++)
-        {
-          bool accept=true;
-
-          // Selecting partons (but not top quark)
-          if (fabs(myEvent->particles()[j].pdgid())>5 && 
-              myEvent->particles()[j].pdgid()!=21) accept=false;
-
-          // Selecting final states
-          if (myEvent->particles()[j].statuscode()!=2) accept=false;
-
-          // Selecting states not coming from proton (beam remnant)
-          if (myEvent->particles()[j].mothup1_==1) accept=false;
-          if (myEvent->particles()[j].mothup1_==2) accept=false;
-
-          const MCParticleFormat* myPart = &(myEvent->particles()[j]);
-          bool test=true;
-          while (myPart->mother1()!=0)
-          {
-            if (myPart->mothup1_==1 || myPart->mothup1_==2)
-            { test=false; break;}
-            else if (myPart->mothup1_<=6)
-            { test=true; break;}
-            else if (myPart->mother1()->pdgid()==91 || 
-               myPart->mother1()->pdgid()==92)
-           {test=false; break;}
-            myPart = myPart->mother1();
-          }
-          if (!test) accept=false;
-
-          if (fabs(myEvent->particles()[j].momentum().Rapidity())>5.0) accept=false;
-
-          double PTJET=sqrt( myEvent->particles()[j].momentum().Px()*myEvent->particles()[j].momentum().Px() + 
-                             myEvent->particles()[j].momentum().Py()*myEvent->particles()[j].momentum().Py() ); 
-          double ETAJET=fabs(log(std::min((sqrt(PTJET*PTJET+myEvent->particles()[j].momentum().Pz()*myEvent->particles()[j].momentum().Pz())
-                                    +fabs( myEvent->particles()[j].momentum().Pz() ))/PTJET,1e5)));
-
-
-          std::cout << "i=" << j+1 
-                    << " ; id=" << myEvent->particles()[j].pdgid()
-                    << " ; s="  << myEvent->particles()[j].statuscode()
-                    << " ; m1=" << myEvent->particles()[j].mothup1_ 
-                    << " ; m2=" << myEvent->particles()[j].mothup2_
-                    << " ; eta=" << ETAJET
-                    << " ; rap=" << myEvent->particles()[j].momentum().Rapidity();
-
-          if (accept) std::cout << " - ACCEPT";
-          std::cout << std::endl;
-        }
-        //        exit(1);
-      }
-      }
-    */
-    double djr = 0;
-    if (DJRvalues[i]>0) djr = std::log10(sqrt(DJRvalues[i]));
-    DJR_[i].total->Fill(djr);
-    DJR_[i].contribution[njets]->Fill(djr);
-  }
-
-}
-
-void MergingPlots::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files)
-{
-  // Saving plots into file
-  Write_TextFormat(out());
-
-  // Deleting plots
-  for (unsigned int i=0;i<DJR_.size();i++)
-  {
-    DJR_[i].Finalize();
-  }
-  DJR_.clear();
-
-  // Free memory allocation
-  if (JetDefinition_==0) delete JetDefinition_;
-}
-
-
-
-
-Double_t MergingPlots::rapidity(Double_t px, Double_t py, Double_t pz)
-{
-  double PTJET = sqrt( px*px + py*py);
-  return fabs(log(std::min((sqrt(PTJET*PTJET+pz*pz)+fabs(pz ))/PTJET,1e5)));
-}
-
-
-void MergingPlots::ExtractDJRwithFortran(const std::vector<fastjet::PseudoJet>& inputs,std::vector<Double_t>& DJRvalues)
-{
-  double PP[512][4];
-  UNUSED(PP);
-  for (unsigned int i=0;i<inputs.size();i++)
-  {
-    PP[i][0]=inputs[i].px();
-    PP[i][1]=inputs[i].py();
-    PP[i][2]=inputs[i].pz();
-    PP[i][3]=inputs[i].e();
-  }
-  //  int IMODE = 4313;
-  //  int NN = inputs.size();
-  //  if (NN>512) NN=512;
-  //  double ECUT=1.;
-  double Y[512];
-  for (unsigned int i=0;i<512;i++) Y[i]=0;
-  //if (NN!=0) ktclus_(&IMODE,PP,&NN,&ECUT,Y);
-
-  for (unsigned int i=0;i<DJRvalues.size();i++)
-  DJRvalues[i]=Y[i];
-}
-
-void MergingPlots::ExtractDJR(const std::vector<fastjet::PseudoJet>& inputs,std::vector<Double_t>& DJRvalues)
-{
-  // JetDefinition_
-  fastjet::ClusterSequence sequence(inputs, *JetDefinition_);
-  for (unsigned int i=0;i<DJRvalues.size();i++)
-  DJRvalues[i]=sequence.exclusive_dmerge(i);
-}
-
-
-/// Saving merging plots in the text output file
-void MergingPlots::Write_TextFormat(SAFWriter& output)
-{
-  *output.GetStream() << "<MergingPlots>" << std::endl;
-  for (unsigned int i=0;i<DJR_.size();i++)
-  {
-    DJR_[i].total->Write_TextFormat(output.GetStream());
-    for (unsigned int j=0;j<DJR_[i].contribution.size();j++)
-    {
-      DJR_[i].contribution[j]->Write_TextFormat(output.GetStream());
-    }
-  }
-  *output.GetStream() << "</MergingPlots>" << std::endl;
-}
-
-
-Bool_t MergingPlots::SavePlots(const std::string& filename)
-{
-  // Opening output file
-  TFile* output = TFile::Open(filename.c_str(),"UPDATE");
-  if (!output->IsOpen())
-  {
-    ERROR << "Output file called '" << filename 
-          << "' is not found" << endmsg;
-    return false;
-  }
-
-  // Writing ROOT file
-  Write_RootFormat(output);
-  return true;
-}
-
-/// Saving merging plots in the ROOT output file
-void MergingPlots::Write_RootFormat(TFile* output)
-{
-  // Creating folder
-  output->mkdir("merging");
-  output->cd("merging");
-  
-  // Saving plots
-  std::pair<TH1F*,TH1F*> histo;
-  for (unsigned int i=0;i<DJR_.size();i++)
-  {
-    // Contribution
-    for (unsigned int j=0;j<DJR_[i].contribution.size();j++)
-    {
-      std::stringstream str;
-      str << "DJR" << i+1;
-      str << "_" << j << "jet";
-      std::string title;
-      str >> title;
-      histo.first  = new TH1F((title+"_pos").c_str(),"",100,0,100);
-      histo.second = new TH1F((title+"_neg").c_str(),"",100,0,100);
-      DJR_[i].contribution[j]->Write_RootFormat(histo);
-      histo.first->Write();
-      histo.second->Write();
-    }
-
-    // total
-    std::stringstream str;
-    str << "DJR" << i+1;
-    str << "_total";
-    std::string title;
-    str >> title;
-    histo.first  = new TH1F((title+"_pos").c_str(),"",100,0,100);
-    histo.second = new TH1F((title+"_neg").c_str(),"",100,0,100);
-    DJR_[i].total->Write_RootFormat(histo);
-    histo.first->Write();
-    histo.second->Write();
-  }
-  output->cd("");
-}
-
-/// Selecting particles for non-hadronized events
-void MergingPlots::SelectParticles_NonHadronization(std::vector<fastjet::PseudoJet>& inputs, const MCEventFormat* myEvent)
-{
-  for (unsigned int i=6;i<myEvent->particles().size();i++)
-  {
-    // Selecting partons (but not top quark)
-    if (fabs(myEvent->particles()[i].pdgid())>5 && 
-        myEvent->particles()[i].pdgid()!=21) continue;
-
-    // Selecting final states
-    if (myEvent->particles()[i].statuscode()==3) continue;
-
-    // Selecting states not coming from initial proton (beam remnant) 
-    // or hadronization
-    const MCParticleFormat* myPart = &(myEvent->particles()[i]);
-    bool test=true;
-    while (myPart->mother1()!=0)
-    {
-      if (myPart->mothup1_==1 || myPart->mothup1_==2)
-      { test=false; break;}
-      else if (myPart->mothup1_<=6)
-      { test=true; break;}
-      else if (myPart->mother1()->pdgid()==91 || 
-               myPart->mother1()->pdgid()==92)
-      {test=false; break;}
-      myPart = myPart->mother1();
-    }
-    if (!test) continue;
-
-    // Cut on the rapidity
-    double ETAJET = rapidity(myEvent->particles()[i].momentum().Px(),
-                             myEvent->particles()[i].momentum().Py(),
-                             myEvent->particles()[i].momentum().Pz());
-    if (fabs(ETAJET)>5) continue;
-    
-    // add the particle
-    inputs.push_back(fastjet::PseudoJet ( myEvent->particles()[i].px(), 
-                                          myEvent->particles()[i].py(), 
-                                          myEvent->particles()[i].pz(), 
-                                          myEvent->particles()[i].e() ) );
-
-    // labeling the particle
-    inputs.back().set_user_index(i);
-  }
-}
-
-
-/// Selecting particles for non-hadronized events
-void MergingPlots::SelectParticles(std::vector<fastjet::PseudoJet>& inputs,
-                                    const MCEventFormat* myEvent)
-{
-  for (unsigned int i=6;i<myEvent->particles().size();i++)
-  {
-    // Selecting partons (but not top quark)
-    if (fabs(myEvent->particles()[i].pdgid())>5 && 
-        myEvent->particles()[i].pdgid()!=21) continue;
-
-    // Selecting final states
-    if (myEvent->particles()[i].statuscode()!=2) continue;
-
-    // Selecting states not coming from initial proton (beam remnant) 
-    // or hadronization
-    const MCParticleFormat* myPart = &(myEvent->particles()[i]);
-    bool test=true;
-    while (myPart->mother1()!=0)
-    {
-      if (myPart->mothup1_==1 || myPart->mothup1_==2)
-      { test=false; break;}
-      else if (myPart->mothup1_<=6)
-      { test=true; break;}
-      else if (myPart->mother1()->pdgid()==91 || 
-               myPart->mother1()->pdgid()==92)
-      {test=false; break;}
-      myPart = myPart->mother1();
-    }
-    if (!test) continue;
-
-    // Cut on the rapidity
-    double ETAJET = rapidity(myEvent->particles()[i].momentum().Px(),
-                             myEvent->particles()[i].momentum().Py(),
-                             myEvent->particles()[i].momentum().Pz());
-    if (fabs(ETAJET)>5) continue;
-
-    // Remove double counting
-    if (myEvent->particles()[i].mother1()!=0 && myPart->mother2()==0)
-    {
-      if (myEvent->particles()[i].pdgid()==myEvent->particles()[i].mother1()->pdgid() &&
-          myEvent->particles()[i].statuscode()==myEvent->particles()[i].mother1()->statuscode() &&
-          fabs(myEvent->particles()[i].px()-myEvent->particles()[i].mother1()->px())<1e-04 &&
-          fabs(myEvent->particles()[i].py()-myEvent->particles()[i].mother1()->py())<1e-04 &&
-          fabs(myEvent->particles()[i].pz()-myEvent->particles()[i].mother1()->pz())<1e-04 )
-        continue;
-    }
-    
-    // add the particle
-    inputs.push_back(fastjet::PseudoJet ( myEvent->particles()[i].px(), 
-                                          myEvent->particles()[i].py(), 
-                                          myEvent->particles()[i].pz(), 
-                                          myEvent->particles()[i].e() ) );
-
-    // labeling the particle
-    inputs.back().set_user_index(i);
-  }
-}
-
-
-
-/// Number of jets
-UInt_t MergingPlots::ExtractJetNumber( const MCEventFormat* myEvent, 
-                                      MCSampleFormat* mySample)
-{
-  UInt_t njets=0;
-  for (unsigned int i=6;i<myEvent->particles().size();i++)
-  {
-    const MCParticleFormat* myPart = &myEvent->particles()[i];
-
-    // keep particles generated during the matrix element calculation
-    if (myPart->statuscode()!=3) continue;
-
-    // keep only partons
-    if (abs(myPart->pdgid())>merging_nqmatch_ && myPart->pdgid()!=21) continue;
-
-    // keep only jets whose mother is one of the initial parton
-    if (myPart->mother1()==0) continue;
-
-    // coming from initial state ?
-    if (myPart->mothup1_>6 && (myPart->mothup1_==0 || myPart->mothup2_==0)) continue;
-
-    // removing color singlet
-    /*
-    if (merging_nosingrad_)
-    {
-      for (unsigned int j=0;j<myEvent->particles().size();j++)
-      {
-        if (i!=j) continue;
-
-        const MCParticleFormat* myPart2 = &myEvent->particles()[j];
-        
-        // keep particles generated during the matrix element calculation
-        if (myPart2->statuscode()!=3) continue;
-
-        // keep only partons 
-        if ( myPart2->pdgid()!=-myPart->pdgid() && 
-             (myPart2->pdgid()!=21 && myPart->pdgid()!=21)) continue;
-
-        // only final states
-
-      }
-    }
-    */
-
-    // count particle
-    njets++;    
-  }
-  /*
-  if (njets==3)  
-  {
-    for (unsigned int i=0;i<myEvent->particles().size();i++)
-    {
-      if (fabs(myEvent->particles()[i].pdgid())>5 && 
-          myEvent->particles()[i].pdgid()!=21) continue;
-      if (myEvent->particles()[i].statuscode()==1) continue;
-
-      std::cout << "i=" << i+1 
-                << " ; id=" << myEvent->particles()[i].pdgid()
-                << " ; s="  << myEvent->particles()[i].statuscode()
-                << " ; m1=" << myEvent->particles()[i].mothup1_ 
-                << " ; m2=" << myEvent->particles()[i].mothup2_ << std::endl;
-    }
-    exit(1);
-    }
-  */
-  return njets;  
-}
-