main92.cc

// main92.cc is a part of the PYTHIA event generator.
// Copyright (C) 2021 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// Keywords: analysis; root;

// This is a simple test program.
// Modified by Rene Brun and Axel Naumann to put the Pythia::event
// into a TTree.

#include <math.h>

// Header file to access Pythia 8 program elements.
#include "Pythia8/Pythia.h"

// ROOT, for saving Pythia events as trees in a file.
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
#include "TColor.h"

#include "EventGeometry.hh"

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"

using namespace Pythia8;

static int nEvents=1;

/*
static double norm_n2_b05 = 1.89612;
static double norm_n256_b05 = 1.89612;

static double norm_n2_b10 = 1./(1. - (1/std::sqrt(3.)));
static double norm_n256_b10 = M_PI/(M_PI-2.);

static double norm_n2_b20 = 1./( (3./2.) - ( 2./std::sqrt(3) ) );
static double norm_n256_b20 = M_PI/((3.*M_PI/2.) - 4.);
*/

static double norm_n2_b05 = 1/0.433786;
static double norm_n256_b05 = 1/0.433786;

static double norm_n2_b10 = 1/0.222521;
static double norm_n256_b10 = 1/0.210397;
  
static double norm_n2_b20 = 1/0.0672669;
static double norm_n256_b20 = 1/0.0578312;

//static int the_tune = 21;

std::vector<std::vector<fastjet::PseudoJet>> generate_events(int n, bool do_backgorund=true, int the_tune=21, int the_batch=0, int do_particle=0);

using EMD = fastjet::contrib::eventgeometry::EMDFloat64<fastjet::contrib::eventgeometry::TransverseMomentum,
							fastjet::contrib::eventgeometry::DeltaR>;

// `PairwiseEMD` is a templated class accepting a fully qualified `EMD` type
//using PairwiseEMD = fastjet::contrib::eventgeometry::PairwiseEMD<EMD>;

// `EMDFloat64` uses `double` for the floating-point type
// first template parameter is an Event type, second is a PairwiseDistance type
using EMDCosPhi = fastjet::contrib::eventgeometry::EMDFloat64<fastjet::contrib::eventgeometry::TransverseMomentum,
							      fastjet::contrib::eventgeometry::CosPhi>;
// demonstrate calculating single EMDs
//EMDCosPhi        IRing_EMD_Obj_Beta05(2*M_PI, 0.5, true); // R, beta, norm
//EMDCosPhi        IRing_EMD_Obj_Beta10(2*M_PI, 1.0, true); // R, beta, norm
//EMDCosPhi        IRing_EMD_Obj_Beta20(2*M_PI, 2.0, true); // R, beta, norm

EMDCosPhi        IRing_EMD_Obj_Beta05(1., 1.0/2., true); // R, beta, norm
EMDCosPhi        IRing_EMD_Obj_Beta10(1., 2.0/2., true); // R, beta, norm
EMDCosPhi        IRing_EMD_Obj_Beta20(1., 4.0/2., true); // R, beta, norm

// empty angle brackets use the default floating-point type `double`
// using CorrelationDimension = fastjet::contrib::eventgeometry::CorrelationDimension<>;

class IsoMin {
public:
  
  using Params = std::valarray<double>;
  using FixedParams = std::valarray<std::vector<fastjet::PseudoJet>>;
  using FuncT = std::function<double(const Params&, const FixedParams&)>;
  
  IsoMin(const FuncT& fin, unsigned int ndim,
	 const FixedParams& fixpar, //const RndT & rndin,
	 unsigned int npop=20, unsigned int ngen=20,
	 double margin=0.1)
    : _f(fin), _q(fixpar), //_rnd(rndin),
      _NDim(ndim), _margin(margin),
      _pop(npop), _fit(npop, -1.0), showTrace(false) {}
  
  void guess(const Params & p) {
    _pop.push_back(p);
    limit01(_pop.back());
    _fit.push_back(f(_pop.back()));
  }
  
  double evolve(unsigned int NGen) {
    for ( unsigned n = 0; n < NGen; ++n ) {
      // Calculate the fitness.
      auto mm = minmax();
      // Always kill the fittest individual.
      if ( showTrace ) _debug();
      for ( unsigned int i = 1; i < _pop.size(); ++i ) {
	if ( _fit[i] > rnd()*(mm.second - mm.first) )
	  // Kill all individuals that have low fitness or are just unlucky.
	  _fit[i] = -1.0;
	else
	  // Improve This individual to be more like the fittest.
	  move(_pop[i],_pop[0]);
      }
    }
    return _fit[0];
  }
  
  Params fittest() const {
    return _pop[0];
  }
  
  double fit() const {
    return _fit[0];
  }
  
  double rnd() const {
    return ((double) rand() / (RAND_MAX)) + 1; // should be on 0-1
  }
  
  Params rndParams() const {
    Params ret(_NDim);
    for ( unsigned int i = 0; i < _NDim; ++i ) ret[i] = rnd();
    return ret;
  }
  
  void limit01(Params & p) const {
    for ( unsigned int i = 0; i < _NDim; ++i )
      p[i] = std::max(0.0, std::min(p[i], 1.0));
  }
  
  void move(Params & bad, const Params & better) const {
    bad += (better - bad)*(rndParams()*(1.0 + 2.0*_margin) - _margin);
    limit01(bad);
  }
  
  double f(const Params & p) const {
    return _f(p, _q);
  }
  
  std::pair<double, double> minmax() {
    std::pair<double,double> mm(std::numeric_limits<double>::max(), 0.0);
    unsigned int iwin = 0;
    for ( unsigned int i = 0; i < _pop.size(); ++i ) {
      double & v = _fit[i];
      // negative fitness value means the individual is dead, so we
      // welocme a new immigrant.
      if ( v < 0.0 ) _pop[i] = rndParams();
      
      // The calculated fitness value cannot be negative.
      v = std::max(0.0, f(_pop[i]));
      
      // Compare to the best and worst fitness so far.
      if ( v < mm.first ) iwin = i;
      mm.first = std::min(v, mm.first);
      mm.second = std::max(mm.second, v);
    }
    
    // Move the winner to the top.
    if ( iwin != 0 ) {
      std::swap(_pop[0], _pop[iwin]);
      std::swap(_fit[0], _fit[iwin]);
    }
    return mm;
  }
  
  void _debug() {
    std::cout << "GenAlgMax population status:" << std::endl;
    for ( unsigned int i = 0; i < _pop.size(); ++i ) {
      std::cout << std::setw(10) << _fit[i] << " (" << _pop[i][0];
      for ( unsigned int ip = 1; ip < _NDim; ++ip )
	std::cout << "," << _pop[i][ip];
      std::cout << ")" << std::endl;
    }
  }
  
private:
  
  const FuncT _f;
  FixedParams _q;
  //  const double _q;
  //const RndT _rnd;
  unsigned int _NDim;
  double _margin;
  std::vector<Params> _pop;
  std::vector<double> _fit;
  
public:
  
  bool showTrace;
  
};


// Defines just a ring of particles, phi values
static std::vector<fastjet::PseudoJet> ringGen(unsigned int phiSeg, double shift=0)
{
  std::vector<fastjet::PseudoJet> phiVals;
  if(phiSeg<1)
    {
      std::cout << "phiSeg must be positive and >0" << std::endl;
      return phiVals;
    }
  for (unsigned int i=0; i<phiSeg; i++)
    {
      double the_phi_value = 2*M_PI*(i+0.5)/phiSeg+shift;
      while(the_phi_value < 0) the_phi_value+=2*M_PI;      
      while(the_phi_value > (2*M_PI)) the_phi_value-=2*M_PI;
      phiVals.push_back(fastjet::PtYPhiM(1.0/phiSeg, 0, the_phi_value, 0));
    }
  return phiVals;
}
    
static double f1_IRing2_Beta05(double shift, std::vector<fastjet::PseudoJet> pjs_transverse)
{
  std::vector<fastjet::PseudoJet> Quasi_Uniform_Event_IRing2;
  Quasi_Uniform_Event_IRing2 = ringGen(2, shift);
  return IRing_EMD_Obj_Beta05(pjs_transverse, Quasi_Uniform_Event_IRing2);
}

static double f2_IRing2_Beta05(const IsoMin::Params& p, const IsoMin::FixedParams& pfix)
{
  double shift = p[0]*2*M_PI; // placeholder
  return f1_IRing2_Beta05(shift, pfix[0]);
}

static double f1_IRing256_Beta05(double shift, std::vector<fastjet::PseudoJet> pjs_transverse)
{
  std::vector<fastjet::PseudoJet> Quasi_Uniform_Event_IRing256;
  Quasi_Uniform_Event_IRing256 = ringGen(256, shift);
  return IRing_EMD_Obj_Beta05(pjs_transverse, Quasi_Uniform_Event_IRing256);
}

static double f2_IRing256_Beta05(const IsoMin::Params& p, const IsoMin::FixedParams& pfix)
{
  double shift = p[0]*2*M_PI; // placeholder
  return f1_IRing256_Beta05(shift, pfix[0]);
}

static double f1_IRing2_Beta10(double shift, std::vector<fastjet::PseudoJet> pjs_transverse)
{
  std::vector<fastjet::PseudoJet> Quasi_Uniform_Event_IRing2;
  Quasi_Uniform_Event_IRing2 = ringGen(2, shift);
  return IRing_EMD_Obj_Beta10(pjs_transverse, Quasi_Uniform_Event_IRing2);
}

static double f2_IRing2_Beta10(const IsoMin::Params& p, const IsoMin::FixedParams& pfix)
{
  double shift = p[0]*2*M_PI; // placeholder
  return f1_IRing2_Beta10(shift, pfix[0]);
}

static double f1_IRing256_Beta10(double shift, std::vector<fastjet::PseudoJet> pjs_transverse)
{
  std::vector<fastjet::PseudoJet> Quasi_Uniform_Event_IRing256;
  Quasi_Uniform_Event_IRing256 = ringGen(256, shift);
  return IRing_EMD_Obj_Beta10(pjs_transverse, Quasi_Uniform_Event_IRing256);
}

static double f2_IRing256_Beta10(const IsoMin::Params& p, const IsoMin::FixedParams& pfix)
{
  double shift = p[0]*2*M_PI; // placeholder
  return f1_IRing256_Beta10(shift, pfix[0]);
}

static double f1_IRing2_Beta20(double shift, std::vector<fastjet::PseudoJet> pjs_transverse)
{
  std::vector<fastjet::PseudoJet> Quasi_Uniform_Event_IRing2;
  Quasi_Uniform_Event_IRing2 = ringGen(2, shift);
  return IRing_EMD_Obj_Beta20(pjs_transverse, Quasi_Uniform_Event_IRing2);
}

static double f2_IRing2_Beta20(const IsoMin::Params& p, const IsoMin::FixedParams& pfix)
{
  double shift = p[0]*2*M_PI; // placeholder
  return f1_IRing2_Beta20(shift, pfix[0]);
}

static double f1_IRing256_Beta20(double shift, std::vector<fastjet::PseudoJet> pjs_transverse)
{
  std::vector<fastjet::PseudoJet> Quasi_Uniform_Event_IRing256;
  Quasi_Uniform_Event_IRing256 = ringGen(256, shift);
  return IRing_EMD_Obj_Beta20(pjs_transverse, Quasi_Uniform_Event_IRing256);
}

static double f2_IRing256_Beta20(const IsoMin::Params& p, const IsoMin::FixedParams& pfix)
{
  double shift = p[0]*2*M_PI; // placeholder
  return f1_IRing256_Beta20(shift, pfix[0]);
}

int main(int argc, char *argv[])
{
  int the_tune = -1;
  int the_batch = -1;
  int do_particle = -1;
  
  cout << "You have entered " << argc
       << " arguments:" << "\n";
  
  for (int i = 0; i < argc; ++i)
    {
      // the tune
      if(i==1)
	{	 
	  cout << "Tune: " << argv[i] << "\n";
	  the_tune = atoi(argv[i]);
	}

      // the batch
      if(i==2)
	{
          cout << "Batch: " << argv[i] << "\n";
          the_batch = atoi(argv[i]);
        }

      // do particle level EMDs (slow!!)
      if(i==3)
        {
          cout << "DoParticle: " << argv[i] << "\n";
          do_particle = atoi(argv[i]);
        }
    }
  
  // demonstrate calculating pairwise EMDs
  //PairwiseEMD pairwise_emd_obj(pow((M_PI*M_PI+(4.5*4.5)),0.5), 1.0, true);

  // preprocess events to center
  //pairwise_emd_obj.preprocess<fastjet::contrib::eventgeometry::CenterWeightedCentroid>();

  // vector of events
  std::vector<std::vector<fastjet::PseudoJet>> background_events;
  //std::vector<std::vector<fastjet::PseudoJet>> signal_events;
  
  // Generate events; skip if generation aborted.
  background_events = generate_events(nEvents, true, the_tune, the_batch, do_particle);
  //signal_events	    = generate_events(1000, false);
  
  // print description
  //std::cout << pairwise_emd_obj.description() << std::endl;

  // run EMD calculations, background
  
  //pairwise_emd_obj(background_events);
  //const std::vector<double> & bkg_emds(pairwise_emd_obj.emds(true));
  /*
  TFile *file = TFile::Open("histos.root","recreate");  
  TH1D* h1_bkg_EMDs     = new TH1D("h1_bkg_EMDs","h1_bkg_EMDs",500,0,1);
  TH1D* h1_bkg_EMDs_avg = new TH1D("h1_bkg_EMDs_avg","h1_bkg_EMDs_avg",500,0,1);
  for(unsigned int iEMD=0; iEMD<bkg_emds.size(); iEMD++)
    {
      h1_bkg_EMDs->Fill(bkg_emds.at(iEMD),1.0);
    }

  
  h1_bkg_EMDs->GetXaxis()->SetTitle("EMD [GeV]");
  h1_bkg_EMDs->GetYaxis()->SetTitle("Entries");
  h1_bkg_EMDs->SetLineWidth(2);
  h1_bkg_EMDs->SetLineColor(9);//TColor::GetColor("#FF6F61"));
  h1_bkg_EMDs->Scale(1/h1_bkg_EMDs->Integral());
  h1_bkg_EMDs->Write();
  */
  // Correlation dimension calculation
  /*
  CorrelationDimension corrdim(50, 10., 250.);
  pairwise_emd_obj.set_external_emd_handler(corrdim);

  // re-run computation
  pairwise_emd_obj(background_events);
  
  // print out correlation dimensions
  auto corrdims(corrdim.corrdims());
  auto corrdim_bins(corrdim.corrdim_bins());
  std::cout << "\nEMD         Corr. Dim.  Error\n" << std::left;
  for (unsigned i = 0; i < corrdims.first.size(); i++)
    std::cout << std::setw(12) << corrdim_bins[i]
              << std::setw(12) << corrdims.first[i]
              << std::setw(12) << corrdims.second[i]
              << '\n';
  */
 
  // Done.
  //delete file;
  return 0;
}

std::vector<std::vector<fastjet::PseudoJet>> generate_events(int n, bool do_background, int the_tune, int the_batch, int do_particle)
{
  // Create Pythia instance and set it up to generate hard QCD processes
  // above pTHat = 20 GeV for pp collisions at 14 TeV.
  Pythia pythia;
  if(do_background)
    {
      pythia.readString("HardQCD:all = on");
      //pythia.readString("SoftQCD:all = on");
      //pythia.readString("Top:gg2ttbar = on");
      //pythia.readString("SLHA:file=sps1aNarrowStopGluinoRPV.spc");
    }
  else // different signals ...
    {
      //pythia.readString("Top:gg2ttbar = on");
      //pythia.readString("ExtraDimensionsG*:all = on");
      //pythia.readString("SUSY:gg2squarkantisquark = on");
      //pythia.readString("SUSY:idA = 1000006");
      //pythia.readString("SUSY:idB = 1000006");
      //pythia.readString("SLHA:file=mc15_13TeV.375855.MGPy8EG_A14N_GG_ttn1_1600_5000_1200.spc");
      //pythia.readString("SLHA:file=sps1aNarrowStopGluinoRPV.spc");
      //pythia.readString("SUSY:gg2squarkantisquark = on");
    }

  // Settings for any process
  pythia.readString("PhaseSpace:pTHatMin = 500.");
  pythia.readString("PhaseSpace:bias2Selection = on");
  pythia.readString("Beams:eCM = 13000.");
  pythia.readString("Tune:pp = "+std::to_string(the_tune));  
  pythia.init();

  // Fastjet analysis - select algorithm and parameters
  double Rparam = 0.4;
  fastjet::JetDefinition         *jetDef = NULL;
  jetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, Rparam);
  
  // Set up the ROOT TFile and TTree.
  TString tree_str = "";
  if(do_background) tree_str = "bkg";
  else tree_str = "sig";
  TFile *file = TFile::Open("out_20220818/emdtree_"+tree_str+"_TUNE"+std::to_string(the_tune)+"_BATCH"+std::to_string(the_batch)+"_PARTICLE"+std::to_string(do_particle)+".root",
			    "recreate");
  //Event *event = &pythia.event;

  TTree *T = new TTree("EMDTree","emd_tree");
  //T->Branch("event",&event);

  int m_njets=-1;
  T->Branch("njets", &m_njets, "njets/I");

  int m_nparticles=-1;
  T->Branch("nparticles", &m_nparticles, "nparticles/I");

  int m_njets_PartonLevel=-1;
  T->Branch("njets_PartonLevel", &m_njets_PartonLevel, "njets_PartonLevel/I");

  int m_nparticles_PartonLevel=-1;
  T->Branch("nparticles_PartonLevel", &m_nparticles_PartonLevel, "nparticles_PartonLevel/I");
  
  double m_jet_0_pt=-1;
  T->Branch("jet_0_pt", &m_jet_0_pt, "jet_0_pt/D");
  double m_jet_1_pt=-1;
  T->Branch("jet_1_pt", &m_jet_1_pt, "jet_1_pt/D");
  double m_jet_2_pt=-1;
  T->Branch("jet_2_pt", &m_jet_2_pt, "jet_2_pt/D");

  double m_jet_0_eta=-1;
  T->Branch("jet_0_eta", &m_jet_0_eta, "jet_0_eta/D");
  double m_jet_1_eta=-1;
  T->Branch("jet_1_eta", &m_jet_1_eta, "jet_1_eta/D");
  double m_jet_2_eta=-1;
  T->Branch("jet_2_eta", &m_jet_2_eta, "jet_2_eta/D");

  double m_j1j2_dphi=-1;
  T->Branch("j1j2_dphi", &m_j1j2_dphi, "m_j1j2_dphi/D");
  double m_j1j2_dr=-1;
  T->Branch("j1j2_dr", &m_j1j2_dr, "m_j1j2_dr/D");

  double m_jet_PartonLevel_0_pt=-1;
  T->Branch("jet_PartonLevel_0_pt", &m_jet_PartonLevel_0_pt, "jet_PartonLevel_0_pt/D");
  double m_jet_PartonLevel_1_pt=-1;
  T->Branch("jet_PartonLevel_1_pt", &m_jet_PartonLevel_1_pt, "jet_PartonLevel_1_pt/D");
  double m_jet_PartonLevel_2_pt=-1;
  T->Branch("jet_PartonLevel_2_pt", &m_jet_PartonLevel_2_pt, "jet_PartonLevel_2_pt/D");

  double m_jet_PartonLevel_0_eta=-1;
  T->Branch("jet_PartonLevel_0_eta", &m_jet_PartonLevel_0_eta, "jet_PartonLevel_0_eta/D");
  double m_jet_PartonLevel_1_eta=-1;
  T->Branch("jet_PartonLevel_1_eta", &m_jet_PartonLevel_1_eta, "jet_PartonLevel_1_eta/D");
  double m_jet_PartonLevel_2_eta=-1;
  T->Branch("jet_PartonLevel_2_eta", &m_jet_PartonLevel_2_eta, "jet_PartonLevel_2_eta/D");

  double m_j1j2_PartonLevel_dphi=-1;
  T->Branch("j1j2_PartonLevel_dphi", &m_j1j2_PartonLevel_dphi, "m_j1j2_PartonLevel_dphi/D");
  double m_j1j2_PartonLevel_dr=-1;
  T->Branch("j1j2_PartonLevel_dr", &m_j1j2_PartonLevel_dr, "m_j1j2_PartonLevel_dr/D");
  
  double m_emdParticles_n2_b05=-1;
  T->Branch("emdParticles_n2_b05",&m_emdParticles_n2_b05, "emdParticles_n2_b05/D");
  double m_emdParticles_n2_b10=-1;
  T->Branch("emdParticles_n2_b10",&m_emdParticles_n2_b10, "emdParticles_n2_b10/D");
  double m_emdParticles_n2_b20=-1;
  T->Branch("emdParticles_n2_b20",&m_emdParticles_n2_b20, "emdParticles_n2_b20/D");
  
  double m_emdParticles_n256_b05=-1;
  T->Branch("emdParticles_n256_b05",&m_emdParticles_n256_b05, "emdParticles_n256_b05/D");
  double m_emdParticles_n256_b10=-1;
  T->Branch("emdParticles_n256_b10",&m_emdParticles_n256_b10, "emdParticles_n256_b10/D");
  double m_emdParticles_n256_b20=-1;
  T->Branch("emdParticles_n256_b20",&m_emdParticles_n256_b20, "emdParticles_n256_b20/D");
  
  double m_emdJets_n2_b05=-1;
  T->Branch("emdJets_n2_b05",&m_emdJets_n2_b05, "emdJets_n2_b05/D");
  double m_emdJets_n2_b10=-1;
  T->Branch("emdJets_n2_b10",&m_emdJets_n2_b10, "emdJets_n2_b10/D");
  double m_emdJets_n2_b20=-1;
  T->Branch("emdJets_n2_b20",&m_emdJets_n2_b20, "emdJets_n2_b20/D");
  
  double m_emdJets_n256_b05=-1;
  T->Branch("emdJets_n256_b05",&m_emdJets_n256_b05, "emdJets_n256_b05/D");
  double m_emdJets_n256_b10=-1;
  T->Branch("emdJets_n256_b10",&m_emdJets_n256_b10, "emdJets_n256_b10/D");
  double m_emdJets_n256_b20=-1;
  T->Branch("emdJets_n256_b20",&m_emdJets_n256_b20, "emdJets_n256_b20/D");  

  double m_emdParticlesPartonLevel_n2_b05=-1;
  T->Branch("emdParticlesPartonLevel_n2_b05", &m_emdParticlesPartonLevel_n2_b05, "emdParticlesPartonLevel_n2_b05/D");
  double m_emdParticlesPartonLevel_n2_b10=-1;
  T->Branch("emdParticlesPartonLevel_n2_b10", &m_emdParticlesPartonLevel_n2_b10, "emdParticlesPartonLevel_n2_b10/D");
  double m_emdParticlesPartonLevel_n2_b20=-1;
  T->Branch("emdParticlesPartonLevel_n2_b20", &m_emdParticlesPartonLevel_n2_b20, "emdParticlesPartonLevel_n2_b20/D");
  
  double m_emdParticlesPartonLevel_n256_b05=-1;
  T->Branch("emdParticlesPartonLevel_n256_b05", &m_emdParticlesPartonLevel_n256_b05, "emdParticlesPartonLevel_n256_b05/D");
  double m_emdParticlesPartonLevel_n256_b10=-1;
  T->Branch("emdParticlesPartonLevel_n256_b10", &m_emdParticlesPartonLevel_n256_b10, "emdParticlesPartonLevel_n256_b10/D");
  double m_emdParticlesPartonLevel_n256_b20=-1;
  T->Branch("emdParticlesPartonLevel_n256_b20", &m_emdParticlesPartonLevel_n256_b20, "emdParticlesPartonLevel_n256_b20/D");
  
  double m_emdJetsPartonLevel_n2_b05=-1;
  T->Branch("emdJetsPartonLevel_n2_b05", &m_emdJetsPartonLevel_n2_b05, "emdJetsPartonLevel_n2_b05/D");
  double m_emdJetsPartonLevel_n2_b10=-1;
  T->Branch("emdJetsPartonLevel_n2_b10", &m_emdJetsPartonLevel_n2_b10, "emdJetsPartonLevel_n2_b10/D");
  double m_emdJetsPartonLevel_n2_b20=-1;
  T->Branch("emdJetsPartonLevel_n2_b20", &m_emdJetsPartonLevel_n2_b20, "emdJetsPartonLevel_n2_b20/D");
  
  double m_emdJetsPartonLevel_n256_b05=-1;
  T->Branch("emdJetsPartonLevel_n256_b05", &m_emdJetsPartonLevel_n256_b05, "emdJetsPartonLevel_n256_b05/D");
  double m_emdJetsPartonLevel_n256_b10=-1;
  T->Branch("emdJetsPartonLevel_n256_b10", &m_emdJetsPartonLevel_n256_b10, "emdJetsPartonLevel_n256_b10/D");
  double m_emdJetsPartonLevel_n256_b20=-1;
  T->Branch("emdJetsPartonLevel_n256_b20", &m_emdJetsPartonLevel_n256_b20, "emdJetsPartonLevel_n256_b20/D");
  
  std::vector<std::vector<fastjet::PseudoJet>> events;  
  
  std::vector<fastjet::PseudoJet> the_ring = ringGen(256, 0);
  std::valarray<std::vector<fastjet::PseudoJet>> the_ring_pfix = {the_ring};

  std::vector<fastjet::PseudoJet> the_triangle = ringGen(3, 0);
  std::valarray<std::vector<fastjet::PseudoJet>> the_triangle_pfix = {the_triangle};

  std::vector<fastjet::PseudoJet> the_square = ringGen(4, 0);
  std::valarray<std::vector<fastjet::PseudoJet>> the_square_pfix = {the_square};
  
  IsoMin mmIRing2_IRing256_Beta05(f2_IRing2_Beta05, 1, the_ring_pfix);
  const double best_IRing2_IRing256_Beta05 = mmIRing2_IRing256_Beta05.evolve(20);
  std::valarray<double> fittest_IRing2_IRing256_Beta05 = mmIRing2_IRing256_Beta05.fittest();
  float m_emdJets_n2_n256_b05 = f2_IRing2_Beta05(fittest_IRing2_IRing256_Beta05, the_ring_pfix);
  
  IsoMin mmIRing2_IRing256_Beta10(f2_IRing2_Beta10, 1, the_ring_pfix);
  const double best_IRing2_IRing256_Beta10 = mmIRing2_IRing256_Beta10.evolve(20);
  std::valarray<double> fittest_IRing2_IRing256_Beta10 = mmIRing2_IRing256_Beta10.fittest();
  float m_emdJets_n2_n256_b10 = f2_IRing2_Beta10(fittest_IRing2_IRing256_Beta10, the_ring_pfix);
  
  IsoMin mmIRing2_IRing256_Beta20(f2_IRing2_Beta20, 1, the_ring_pfix);
  const double best_IRing2_IRing256_Beta20 = mmIRing2_IRing256_Beta20.evolve(20);
  std::valarray<double> fittest_IRing2_IRing256_Beta20 = mmIRing2_IRing256_Beta20.fittest();
  float m_emdJets_n2_n256_b20 = f2_IRing2_Beta20(fittest_IRing2_IRing256_Beta20, the_ring_pfix);


  IsoMin mmIRing2_IRing3_Beta05(f2_IRing2_Beta05, 1, the_triangle_pfix);
  const double best_IRing2_IRing3_Beta05 = mmIRing2_IRing3_Beta05.evolve(20);
  std::valarray<double> fittest_IRing2_IRing3_Beta05 = mmIRing2_IRing3_Beta05.fittest();
  float m_emdJets_n2_n3_b05 = f2_IRing2_Beta05(fittest_IRing2_IRing3_Beta05, the_triangle_pfix);

  IsoMin mmIRing2_IRing3_Beta10(f2_IRing2_Beta10, 1, the_triangle_pfix);
  const double best_IRing2_IRing3_Beta10 = mmIRing2_IRing3_Beta10.evolve(20);
  std::valarray<double> fittest_IRing2_IRing3_Beta10 = mmIRing2_IRing3_Beta10.fittest();
  float m_emdJets_n2_n3_b10 = f2_IRing2_Beta10(fittest_IRing2_IRing3_Beta10, the_triangle_pfix);

  IsoMin mmIRing2_IRing3_Beta20(f2_IRing2_Beta20, 1, the_triangle_pfix);
  const double best_IRing2_IRing3_Beta20 = mmIRing2_IRing3_Beta20.evolve(20);
  std::valarray<double> fittest_IRing2_IRing3_Beta20 = mmIRing2_IRing3_Beta20.fittest();
  float m_emdJets_n2_n3_b20 = f2_IRing2_Beta20(fittest_IRing2_IRing3_Beta20, the_triangle_pfix);

  IsoMin mmIRing2_IRing4_Beta05(f2_IRing2_Beta05, 1, the_square_pfix);
  const double best_IRing2_IRing4_Beta05 = mmIRing2_IRing4_Beta05.evolve(20);
  std::valarray<double> fittest_IRing2_IRing4_Beta05 = mmIRing2_IRing4_Beta05.fittest();
  float m_emdJets_n2_n4_b05 = f2_IRing2_Beta05(fittest_IRing2_IRing4_Beta05, the_square_pfix);

  IsoMin mmIRing2_IRing4_Beta10(f2_IRing2_Beta10, 1, the_square_pfix);
  const double best_IRing2_IRing4_Beta10 = mmIRing2_IRing4_Beta10.evolve(20);
  std::valarray<double> fittest_IRing2_IRing4_Beta10 = mmIRing2_IRing4_Beta10.fittest();
  float m_emdJets_n2_n4_b10 = f2_IRing2_Beta10(fittest_IRing2_IRing4_Beta10, the_square_pfix);

  IsoMin mmIRing2_IRing4_Beta20(f2_IRing2_Beta20, 1, the_square_pfix);
  const double best_IRing2_IRing4_Beta20 = mmIRing2_IRing4_Beta20.evolve(20);
  std::valarray<double> fittest_IRing2_IRing4_Beta20 = mmIRing2_IRing4_Beta20.fittest();
  float m_emdJets_n2_n4_b20 = f2_IRing2_Beta20(fittest_IRing2_IRing4_Beta20, the_square_pfix);

  
  
  /*
  std::cout << "\nDistance between references is"
	    << "\n\tBeta=0.5: " << m_emdJets_n2_n256_b05
	    << "\n\tBeta=1.0: " << m_emdJets_n2_n256_b10
	    << "\n\tBeta=2.0: "	<< m_emdJets_n2_n256_b20
	    << std::endl;
  */
  
  // Begin event loop. Generate event; skip if generation aborted.
  bool firstEvent = true;
  for (int iEvent = 0; iEvent < n; ++iEvent)
    {
      //std::cout << "Processing event " << iEvent << "/" << nEvents << std::endl;
            
      if (!pythia.next()) continue;

      std::vector <fastjet::PseudoJet> fjInputs_PartonLevel;
      std::vector <fastjet::PseudoJet> fjInputs_PartonLevel_transverse;

      std::vector <fastjet::PseudoJet> fjInputs;
      std::vector <fastjet::PseudoJet> fjInputs_transverse;
      
      // Loop over event record to decide what to pass to FastJet
      for (int i = 0; i < pythia.event.size(); ++i)
	{
	  // Final state only
	  if (! (pythia.event[i].isFinal() || pythia.event[i].isFinalPartonLevel()) )
	    continue;
	  
	  // No neutrinos
	  if (pythia.event[i].idAbs() == 12
	      || pythia.event[i].idAbs() == 14
	      || pythia.event[i].idAbs() == 16)
	    continue;
	  
	  // Only |eta| < 4.9
	  //std::cout << pythia.event.size() << " " << pythia.event[i].pT() << std::endl;
	  if (pythia.event[i].pT() < 0.5)
	  if (abs(pythia.event[i].eta()) > 4.9)
	    continue;
	  
	  // Store as input to Fastjet
	  if(pythia.event[i].isFinalPartonLevel())
	    fjInputs_PartonLevel.push_back( fastjet::PseudoJet( pythia.event[i].px(),
								pythia.event[i].py(),
								pythia.event[i].pz(),
								pythia.event[i].e() ) );

	  if(pythia.event[i].isFinal())
	    fjInputs.push_back( fastjet::PseudoJet( pythia.event[i].px(),
						    pythia.event[i].py(),
						    pythia.event[i].pz(),
						    pythia.event[i].e() ) );
	}

      // particle multiplicity
      m_nparticles = fjInputs.size();
      
      if (fjInputs.size() == 0 && fjInputs_PartonLevel.size() == 0)
	{
	  cout << "Error: event with no final state particles" << endl;
	  continue;
	}
  
      // Run Fastjet algorithm
      std::vector <fastjet::PseudoJet> inclusiveJets, sortedJets, pjs_transverse, pjs_transverse_particles;
      fastjet::PseudoJet vectorial_sum(fastjet::PtYPhiM(0,0,0,0));
      fastjet::PseudoJet vectorial_sum_particles(fastjet::PtYPhiM(0,0,0,0));      
      fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);

      /*
      std::vector <fastjet::PseudoJet> inclusiveJets_PartonLevel, sortedJets_PartonLevel, pjs_transverse_PartonLevel, pjs_transverse_particles_PartonLevel;
      fastjet::PseudoJet vectorial_sum_PartonLevel(fastjet::PtYPhiM(0,0,0,0));
      fastjet::PseudoJet vectorial_sum_particles_PartonLevel(fastjet::PtYPhiM(0,0,0,0));
      fastjet::ClusterSequence clustSeq_PartonLevel(fjInputs_PartonLevel, *jetDef);
      */
      
      // For the first event, print the FastJet details
      if (firstEvent)
	{
	  cout << "Ran " << jetDef->description() << endl;
	  //cout << "Strategy adopted by FastJet was "
	  //     << clustSeq.strategy_string() << endl << endl;
	  firstEvent = false;
	}
  
      // Extract inclusive jets sorted by pT (note minimum pT of 20.0 GeV)
      inclusiveJets = clustSeq.inclusive_jets(60.0);
      sortedJets    = sorted_by_pt(inclusiveJets);

      //inclusiveJets_PartonLevel = clustSeq_PartonLevel.inclusive_jets(20.0);
      //sortedJets_PartonLevel    = sorted_by_pt(inclusiveJets_PartonLevel);
      
      // get the individual particle vectors
      /*
      for (unsigned int i=0; i < fjInputs.size(); i++)
	{
	  vectorial_sum_particles-=fastjet::PtYPhiM(fjInputs[i].pt(),
						    fjInputs[i].rapidity(),
						    fjInputs[i].phi(),
						    fjInputs[i].m());
	  pjs_transverse_particles.push_back(fastjet::PtYPhiM(fjInputs[i].perp(),
							      0,
							      fjInputs[i].phi(),
							      0));
	}
      pjs_transverse_particles.push_back(fastjet::PtYPhiM(vectorial_sum_particles.pt(), 0, vectorial_sum_particles.phi(), 0));


      for (unsigned int i=0; i < fjInputs_PartonLevel.size(); i++)
        {
          vectorial_sum_particles_PartonLevel-=fastjet::PtYPhiM(fjInputs_PartonLevel[i].pt(),
								fjInputs_PartonLevel[i].rapidity(),
								fjInputs_PartonLevel[i].phi(),
								fjInputs_PartonLevel[i].m());
          pjs_transverse_particles_PartonLevel.push_back(fastjet::PtYPhiM(fjInputs_PartonLevel[i].perp(),
									  0,
									  fjInputs_PartonLevel[i].phi(),
									  0));
        }
      pjs_transverse_particles_PartonLevel.push_back(fastjet::PtYPhiM(vectorial_sum_particles_PartonLevel.pt(),
								      0,
								      vectorial_sum_particles_PartonLevel.phi(),
								      0));
      */
      
      // Keep track of jets with pT > 60 GeV
      int  jetCount60 = 0;
      for (unsigned int i = 0; i < sortedJets.size(); i++)
	{
	  // Only count jets that have |eta| < 2.0
	  if (abs(sortedJets[i].rap()) > 4.5) continue;
	  // Fill histograms and count jets with ET > 20.0
	  if (sortedJets[i].perp() > 60.0)
	    {
	      jetCount60++;
	      pjs_transverse.push_back(fastjet::PtYPhiM(sortedJets[i].perp(), 0, sortedJets[i].phi(), 0));
	      vectorial_sum-=fastjet::PtYPhiM(sortedJets[i].pt(), sortedJets[i].rapidity(), sortedJets[i].phi(), sortedJets[i].m());
	    }
	}
      pjs_transverse.push_back(fastjet::PtYPhiM(vectorial_sum.pt(), 0, vectorial_sum.phi(), 0));


      /*
     int  jetCount60_PartonLevel = 0;
      for (unsigned int i = 0; i < sortedJets_PartonLevel.size(); i++)
        {
          // Only count jets that have |eta| < 2.0
          if (abs(sortedJets_PartonLevel[i].rap()) > 4.5) continue;
          // Fill histograms and count jets with ET > 20.0
          if (sortedJets_PartonLevel[i].perp() > 60.0)
            {
              jetCount60_PartonLevel++;
              pjs_transverse_PartonLevel.push_back(fastjet::PtYPhiM(sortedJets_PartonLevel[i].perp(), 0, sortedJets_PartonLevel[i].phi(), 0));
              vectorial_sum_PartonLevel-=fastjet::PtYPhiM(sortedJets_PartonLevel[i].pt(), sortedJets_PartonLevel[i].rapidity(), sortedJets_PartonLevel[i].phi(), sortedJets_PartonLevel[i].m());
            }
        }
      pjs_transverse_PartonLevel.push_back(fastjet::PtYPhiM(vectorial_sum_PartonLevel.pt(), 0, vectorial_sum_PartonLevel.phi(), 0));
      */
      
      if( (jetCount60<2) ) //&& (jetCount60_PartonLevel<2) ) continue;
      // Fill the pythia event into the TTree.
      // Warning: the files will rapidly become large if all events
      // are saved. In some cases it may be convenient to do some
      // processing of events and only save those that appear
      // interesting for future analyses.

      //std::cout << "A good event!" << std::endl;
      
      //events.push_back(sortedJets);  
      m_njets = jetCount60;
      //m_njets_PartonLevel = jetCount60_PartonLevel;
      
      //std::cout << "nParticles " << m_nparticles << " nJets " << m_njets << std::endl;
      
      // after hadronisation
      if(sortedJets.size()>0)
	{
	  m_jet_0_pt = sortedJets[0].pt();
	  m_jet_0_eta = sortedJets[0].eta();
	}
      if(sortedJets.size()>1)
	{
	  m_jet_1_pt = sortedJets[1].pt();
	  m_jet_1_eta = sortedJets[1].eta();

	  m_j1j2_dphi = sortedJets[0].delta_phi_to(sortedJets[1]);
	  m_j1j2_dr = sortedJets[0].delta_R(sortedJets[1]);
	}
      if(sortedJets.size()>2)
	{
	  m_jet_2_pt = sortedJets[2].pt(); // third jet
	  m_jet_2_eta = sortedJets[2].eta();
	}

      // parton level
      /*
      if(sortedJets_PartonLevel.size()>0)
	{
	  m_jet_PartonLevel_0_pt = sortedJets_PartonLevel[0].pt();
	  m_jet_PartonLevel_0_eta = sortedJets_PartonLevel[0].eta();
	}
      if(sortedJets_PartonLevel.size()>1)
	{
	  m_jet_PartonLevel_1_pt = sortedJets_PartonLevel[1].pt();
	  m_jet_PartonLevel_1_eta = sortedJets_PartonLevel[1].eta();
	  
	  m_j1j2_PartonLevel_dphi = sortedJets_PartonLevel[0].delta_phi_to(sortedJets_PartonLevel[1]);
	  m_j1j2_PartonLevel_dr = sortedJets_PartonLevel[0].delta_R(sortedJets_PartonLevel[1]);
	}
      if(sortedJets_PartonLevel.size()>2)
	{
	  m_jet_PartonLevel_2_pt = sortedJets_PartonLevel[2].pt(); // third jet
	  m_jet_PartonLevel_2_eta = sortedJets_PartonLevel[2].eta();
	}
      */
	
      float shift = sortedJets[0].phi();      
      std::valarray<std::vector<fastjet::PseudoJet>> pfix = {pjs_transverse};

      //float shift_PartonLevel = sortedJets_PartonLevel[0].phi();
      //std::valarray<std::vector<fastjet::PseudoJet>> pfix_PartonLevel = {pjs_transverse_PartonLevel};
      
      // after hadronization
      
      IsoMin mmIRing2_Beta05(f2_IRing2_Beta05, 1, pfix);
      const double best_IRing2_Beta05 = mmIRing2_Beta05.evolve(20);
      std::valarray<double> fittest_IRing2_Beta05 = mmIRing2_Beta05.fittest();
      m_emdJets_n2_b05 = f2_IRing2_Beta05(fittest_IRing2_Beta05, pfix);

      IsoMin mmIRing2_Beta10(f2_IRing2_Beta10, 1, pfix);
      const double best_IRing2_Beta10 = mmIRing2_Beta10.evolve(20);
      std::valarray<double> fittest_IRing2_Beta10 = mmIRing2_Beta10.fittest();
      m_emdJets_n2_b10 = f2_IRing2_Beta10(fittest_IRing2_Beta10, pfix);

      IsoMin mmIRing2_Beta20(f2_IRing2_Beta20, 1, pfix);
      const double best_IRing2_Beta20 = mmIRing2_Beta20.evolve(20);
      std::valarray<double> fittest_IRing2_Beta20 = mmIRing2_Beta20.fittest();
      m_emdJets_n2_b20 = f2_IRing2_Beta20(fittest_IRing2_Beta20, pfix);
      
      IsoMin mmIRing256_Beta05(f2_IRing256_Beta05, 1, pfix);
      const double best_IRing256_Beta05 = mmIRing256_Beta05.evolve(20);
      std::valarray<double> fittest_IRing256_Beta05 = mmIRing256_Beta05.fittest();
      m_emdJets_n256_b05 = f2_IRing256_Beta05(fittest_IRing256_Beta05, pfix);
      
      IsoMin mmIRing256_Beta10(f2_IRing256_Beta10, 1, pfix);
      const double best_IRing256_Beta10 = mmIRing256_Beta10.evolve(20);
      std::valarray<double> fittest_IRing256_Beta10 = mmIRing256_Beta10.fittest();
      m_emdJets_n256_b10 = f2_IRing256_Beta10(fittest_IRing256_Beta10, pfix);

      IsoMin mmIRing256_Beta20(f2_IRing256_Beta20, 1, pfix);
      const double best_IRing256_Beta20 = mmIRing256_Beta20.evolve(20);
      std::valarray<double> fittest_IRing256_Beta20 = mmIRing256_Beta20.fittest();
      m_emdJets_n256_b20 = f2_IRing256_Beta20(fittest_IRing256_Beta20, pfix);

      // parton level stuff

      /*
      
      IsoMin mmIRing2_Beta05_PartonLevel(f2_IRing2_Beta05, 1, pfix_PartonLevel);
      const double best_IRing2_Beta05_PartonLevel = mmIRing2_Beta05_PartonLevel.evolve(20);
      std::valarray<double> fittest_IRing2_Beta05_PartonLevel = mmIRing2_Beta05_PartonLevel.fittest();
      m_emdJetsPartonLevel_n2_b05 = f2_IRing2_Beta05(fittest_IRing2_Beta05, pfix);

      IsoMin mmIRing2_Beta10_PartonLevel(f2_IRing2_Beta10, 1, pfix_PartonLevel);
      const double best_IRing2_Beta10_PartonLevel = mmIRing2_Beta10_PartonLevel.evolve(20);
      std::valarray<double> fittest_IRing2_Beta10_PartonLevel = mmIRing2_Beta10_PartonLevel.fittest();
      m_emdJetsPartonLevel_n2_b10 = f2_IRing2_Beta10(fittest_IRing2_Beta10, pfix_PartonLevel);

      IsoMin mmIRing2_Beta20_PartonLevel(f2_IRing2_Beta20, 1, pfix_PartonLevel);
      const double best_IRing2_Beta20_PartonLevel = mmIRing2_Beta20_PartonLevel.evolve(20);
      std::valarray<double> fittest_IRing2_Beta20_PartonLevel = mmIRing2_Beta20_PartonLevel.fittest();
      m_emdJetsPartonLevel_n2_b20 = f2_IRing2_Beta20(fittest_IRing2_Beta20, pfix_PartonLevel);

      IsoMin mmIRing256_Beta05_PartonLevel(f2_IRing256_Beta05, 1, pfix_PartonLevel);
      const double best_IRing256_Beta05_PartonLevel = mmIRing256_Beta05_PartonLevel.evolve(20);
      std::valarray<double> fittest_IRing256_Beta05_PartonLevel = mmIRing256_Beta05_PartonLevel.fittest();
      m_emdJetsPartonLevel_n256_b05 = f2_IRing256_Beta05(fittest_IRing256_Beta05, pfix_PartonLevel);

      IsoMin mmIRing256_Beta10_PartonLevel(f2_IRing256_Beta10, 1, pfix_PartonLevel);
      const double best_IRing256_Beta10_PartonLevel = mmIRing256_Beta10_PartonLevel.evolve(20);
      std::valarray<double> fittest_IRing256_Beta10_PartonLevel = mmIRing256_Beta10_PartonLevel.fittest();
      m_emdJetsPartonLevel_n256_b10 = f2_IRing256_Beta10(fittest_IRing256_Beta10, pfix_PartonLevel);

      IsoMin mmIRing256_Beta20_PartonLevel(f2_IRing256_Beta20, 1, pfix_PartonLevel);
      const double best_IRing256_Beta20_PartonLevel = mmIRing256_Beta20_PartonLevel.evolve(20);
      std::valarray<double> fittest_IRing256_Beta20_PartonLevel = mmIRing256_Beta20_PartonLevel.fittest();
      m_emdJetsPartonLevel_n256_b20 = f2_IRing256_Beta20(fittest_IRing256_Beta20_PartonLevel, pfix_PartonLevel);

      */
      
      if(iEvent%100==0) // debug info
	{
	  std::cout << "\nEvent " << iEvent << "/" << nEvents << std::endl;

	  std::cout << "There are " << fjInputs_PartonLevel.size() << " particles at parton level, and "
		    << fjInputs.size() << " after hadronisation." << std::endl;
	  
	  std::cout << "There are " << sortedJets.size() << " jets." << std::endl;
	  for (int iJet=0; iJet<sortedJets.size(); iJet++)
	    {
	      std::cout << "\tjet " << iJet
			<< " pt: " <<  sortedJets[iJet].pt()
			<< " eta: " <<  sortedJets[iJet].eta()
			<< " phi: " <<  sortedJets[iJet].phi()
			<< " E: " <<  sortedJets[iJet].E()
			<< std::endl;     
	    }

	  /*
	  std::cout << "There are " << sortedJets_PartonLevel.size() << " parton-level jets." << std::endl;
	  for (int iJet=0; iJet<sortedJets_PartonLevel.size(); iJet++)
            {
              std::cout << "\tjet " << iJet
                        << " pt: " <<  sortedJets_PartonLevel[iJet].pt()
                        << " eta: " <<  sortedJets_PartonLevel[iJet].eta()
			<< " phi: " <<  sortedJets_PartonLevel[iJet].phi()
                        << " E: " <<  sortedJets_PartonLevel[iJet].E()
                        << std::endl;
            }
	  */
	  
	  std::cout << "EMDs:" << std::endl;
	  std::cout << "\tN=2 Beta=0.5: " << m_emdJets_n2_b05*norm_n2_b05 << std::endl;
	  std::cout << "\tN=2 Beta=1.0: " << m_emdJets_n2_b10*norm_n2_b10 << std::endl;
	  std::cout << "\tN=2 Beta=2.0: " << m_emdJets_n2_b20*norm_n2_b20 << std::endl;
	  std::cout << "\tN=256 Beta=0.5: " << m_emdJets_n256_b05*norm_n256_b05 << std::endl;
          std::cout << "\tN=256 Beta=1.0: " << m_emdJets_n256_b10*norm_n256_b10 << std::endl;
          std::cout << "\tN=256 Beta=2.0: " << m_emdJets_n256_b20*norm_n256_b20 << std::endl;

	  /*
	  std::cout << "Parton Level EMDs:" << std::endl;
          std::cout << "\tN=2 Beta=0.5: " << m_emdJetsPartonLevel_n2_b05*norm_n2_b05 << std::endl;
          std::cout << "\tN=2 Beta=1.0: " << m_emdJetsPartonLevel_n2_b10*norm_n2_b10 << std::endl;
          std::cout << "\tN=2 Beta=2.0: " << m_emdJetsPartonLevel_n2_b20*norm_n2_b20 << std::endl;
          std::cout << "\tN=256 Beta=0.5: " << m_emdJetsPartonLevel_n256_b05*norm_n256_b05 << std::endl;
          std::cout << "\tN=256 Beta=1.0: " << m_emdJetsPartonLevel_n256_b10*norm_n256_b10 << std::endl;
          std::cout << "\tN=256 Beta=2.0: " << m_emdJetsPartonLevel_n256_b20*norm_n256_b20 << std::endl;
	  */
	}

      /*
      
      float shift_particles = fjInputs[0].phi();
      std::valarray<std::vector<fastjet::PseudoJet>> pfix_particles = {pjs_transverse_particles};

      IsoMin particle_mmIRing2_Beta05(f2_IRing2_Beta05, 1, pfix_particles);
      const double best_IRing2_Beta05_Particles = particle_mmIRing2_Beta05.evolve(20);
      std::valarray<double> fittest_particles_IRing2_Beta05 = particle_mmIRing2_Beta05.fittest();
      m_emdParticles_n2_b05 = f2_IRing2_Beta05(fittest_particles_IRing2_Beta05, pfix_particles); // GET NORM FACTOR

      IsoMin particle_mmIRing2_Beta10(f2_IRing2_Beta10, 1, pfix_particles);
      const double best_IRing2_Beta10_Particles = particle_mmIRing2_Beta10.evolve(20);
      std::valarray<double> fittest_particles_IRing2_Beta10 = particle_mmIRing2_Beta10.fittest();
      m_emdParticles_n2_b10 = f2_IRing2_Beta10(fittest_particles_IRing2_Beta10, pfix_particles); // GET NORM FACTOR

      IsoMin particle_mmIRing2_Beta20(f2_IRing2_Beta20, 1, pfix_particles);
      const double best_IRing2_Beta20_Particles = particle_mmIRing2_Beta20.evolve(20);
      std::valarray<double> fittest_particles_IRing2_Beta20 = particle_mmIRing2_Beta20.fittest();
      m_emdParticles_n2_b20 = f2_IRing2_Beta20(fittest_particles_IRing2_Beta20, pfix_particles);  // GET NORM FACTOR
      
      IsoMin particle_mmIRing256_Beta05(f2_IRing256_Beta05, 1, pfix_particles);
      const double best_IRing256_Beta05_Particles = particle_mmIRing256_Beta05.evolve(20);
      std::valarray<double> fittest_particles_IRing256_Beta05 = particle_mmIRing256_Beta05.fittest();
      m_emdParticles_n256_b05 = f2_IRing256_Beta05(fittest_particles_IRing256_Beta05, pfix_particles); // GET NORM FACTOR

      IsoMin particle_mmIRing256_Beta10(f2_IRing256_Beta10, 1, pfix_particles);
      const double best_IRing256_Beta10_Particles = particle_mmIRing256_Beta10.evolve(20);
      std::valarray<double> fittest_particles_IRing256_Beta10 = particle_mmIRing256_Beta10.fittest();
      m_emdParticles_n256_b10 = f2_IRing256_Beta10(fittest_particles_IRing256_Beta10, pfix_particles); // GET NORM FACTOR
      
      IsoMin particle_mmIRing256_Beta20(f2_IRing256_Beta20, 1, pfix_particles);
      const double best_IRing256_Beta20_Particles = particle_mmIRing256_Beta20.evolve(20);
      std::valarray<double> fittest_particles_IRing256_Beta20 = particle_mmIRing256_Beta20.fittest();
      m_emdParticles_n256_b20 = f2_IRing256_Beta20(fittest_particles_IRing256_Beta20, pfix_particles);
     
      
      float shift_particles_PartonLevel = fjInputs[0].phi();
      std::valarray<std::vector<fastjet::PseudoJet>> pfix_particles_PartonLevel = {pjs_transverse_particles_PartonLevel};

      IsoMin particle_PartonLevel_mmIRing2_Beta05(f2_IRing2_Beta05, 1, pfix_particles_PartonLevel);
      const double best_IRing2_Beta05_Particles_PartonLevel = particle_PartonLevel_mmIRing2_Beta05.evolve(20);
      std::valarray<double> fittest_particles_PartonLevel_IRing2_Beta05 = particle_PartonLevel_mmIRing2_Beta05.fittest();
      m_emdParticlesPartonLevel_n2_b05 = f2_IRing2_Beta05(fittest_particles_PartonLevel_IRing2_Beta05, pfix_particles_PartonLevel); 

      IsoMin particle_PartonLevel_mmIRing2_Beta10(f2_IRing2_Beta10, 1, pfix_particles_PartonLevel);
      const double best_IRing2_Beta10_Particles_PartonLevel = particle_PartonLevel_mmIRing2_Beta10.evolve(20);
      std::valarray<double> fittest_particles_PartonLevel_IRing2_Beta10 = particle_PartonLevel_mmIRing2_Beta10.fittest();
      m_emdParticlesPartonLevel_n2_b10 = f2_IRing2_Beta10(fittest_particles_PartonLevel_IRing2_Beta10, pfix_particles_PartonLevel);

      IsoMin particle_PartonLevel_mmIRing2_Beta20(f2_IRing2_Beta20, 1, pfix_particles_PartonLevel);
      const double best_IRing2_Beta20_Particles_PartonLevel = particle_PartonLevel_mmIRing2_Beta20.evolve(20);
      std::valarray<double> fittest_particles_PartonLevel_IRing2_Beta20 = particle_PartonLevel_mmIRing2_Beta20.fittest();
      m_emdParticlesPartonLevel_n2_b20 = f2_IRing2_Beta20(fittest_particles_PartonLevel_IRing2_Beta20, pfix_particles_PartonLevel);  

      IsoMin particle_PartonLevel_mmIRing256_Beta05(f2_IRing256_Beta05, 1, pfix_particles_PartonLevel);
      const double best_IRing256_Beta05_Particles_PartonLevel = particle_PartonLevel_mmIRing256_Beta05.evolve(20);
      std::valarray<double> fittest_particles_PartonLevel_IRing256_Beta05 = particle_PartonLevel_mmIRing256_Beta05.fittest();
      m_emdParticlesPartonLevel_n256_b05 = f2_IRing256_Beta05(fittest_particles_PartonLevel_IRing256_Beta05, pfix_particles_PartonLevel); 
      
      IsoMin particle_PartonLevel_mmIRing256_Beta10(f2_IRing256_Beta10, 1, pfix_particles_PartonLevel);
      const double best_IRing256_Beta10_Particles_PartonLevel = particle_PartonLevel_mmIRing256_Beta10.evolve(20);
      std::valarray<double> fittest_particles_PartonLevel_IRing256_Beta10 = particle_PartonLevel_mmIRing256_Beta10.fittest();
      m_emdParticlesPartonLevel_n256_b10 = f2_IRing256_Beta10(fittest_particles_PartonLevel_IRing256_Beta10, pfix_particles_PartonLevel);
      
      IsoMin particle_PartonLevel_mmIRing256_Beta20(f2_IRing256_Beta20, 1, pfix_particles_PartonLevel);
      const double best_IRing256_Beta20_Particles_PartonLevel = particle_PartonLevel_mmIRing256_Beta20.evolve(20);
      std::valarray<double> fittest_particles_PartonLevel_IRing256_Beta20 = particle_PartonLevel_mmIRing256_Beta20.fittest();
      m_emdParticlesPartonLevel_n256_b20 = f2_IRing256_Beta20(fittest_particles_PartonLevel_IRing256_Beta20, pfix_particles_PartonLevel);
      
      
       if(iEvent%100==0) // debug info
        {
	  std::cout << "Particle EMDs:" << std::endl;
          std::cout << "\tN=2 Beta=0.5: " << m_emdParticles_n2_b05*norm_n2_b05 << std::endl;
          std::cout << "\tN=2 Beta=1.0: " << m_emdParticles_n2_b10*norm_n2_b10 << std::endl;
          std::cout << "\tN=2 Beta=2.0: " << m_emdParticles_n2_b20*norm_n2_b20 << std::endl;
          std::cout << "\tN=256 Beta=0.5: " << m_emdParticles_n256_b05*norm_n256_b05 << std::endl;
          std::cout << "\tN=256 Beta=1.0: " << m_emdParticles_n256_b10*norm_n256_b10 << std::endl;
          std::cout << "\tN=256 Beta=2.0: " << m_emdParticles_n256_b20*norm_n256_b20 << std::endl;

          std::cout << "Parton Level Particle EMDs:" << std::endl;
          std::cout << "\tN=2 Beta=0.5: " << m_emdParticlesPartonLevel_n2_b05*norm_n2_b05 << std::endl;
          std::cout << "\tN=2 Beta=1.0: " << m_emdParticlesPartonLevel_n2_b10*norm_n2_b10 << std::endl;
          std::cout << "\tN=2 Beta=2.0: " << m_emdParticlesPartonLevel_n2_b20*norm_n2_b20 << std::endl;
          std::cout << "\tN=256 Beta=0.5: " << m_emdParticlesPartonLevel_n256_b05*norm_n256_b05 << std::endl;
          std::cout << "\tN=256 Beta=1.0: " << m_emdParticlesPartonLevel_n256_b10*norm_n256_b10 << std::endl;
          std::cout << "\tN=256 Beta=2.0: " << m_emdParticlesPartonLevel_n256_b20*norm_n256_b20 << std::endl;
        }
       */
      T->Fill();
      firstEvent=false;
    }   // End event loop.

  // Statistics on event generation.
  pythia.stat();

  //  Write tree.
  //T->Print();
  T->Write();

  std::cout << "\nDistance between references is"
            << "\n\tBeta=0.5, N=2 vs N=256: " << m_emdJets_n2_n256_b05
            << "\n\tBeta=1.0, N=2 vs N=256: " << m_emdJets_n2_n256_b10
            << "\n\tBeta=2.0, N=2 vs N=256: " << m_emdJets_n2_n256_b20
            << std::endl;

   std::cout<< "\n\tBeta=0.5, N=2 vs N=3: " << m_emdJets_n2_n3_b05
            << "\n\tBeta=1.0, N=2 vs N=3: " << m_emdJets_n2_n3_b10
            << "\n\tBeta=2.0, N=2 vs N=3: " << m_emdJets_n2_n3_b20
            << std::endl;

 std::cout<< "\n\tBeta=0.5, N=2 vs N=4: " << m_emdJets_n2_n4_b20
	  << "\n\tBeta=1.0, N=2 vs N=4: " << m_emdJets_n2_n4_b10
	  << "\n\tBeta=2.0, N=2 vs N=4: " << m_emdJets_n2_n4_b20
	  << std::endl;
   
  delete file;
  return events;
}