Example 5: Adding new module
============================

1. Create file modules/IsolationLEP.h with the following content (see below)

2. Create file modules/IsolationLEP.cc with the following content (see below)

3. Add new module to modules/ModulesLinkDef.h:

#include "modules/IsolationLEP.h"
...
#pragma link C++ class IsolationLEP+;

4. Add DeltaTheta variable to Candidate class and Muon class in classes/DelphesClasses.h
and in classes/DelphesClasses.cc

5. Modify TreeWriter::ProcessMuons() in modules/TreeWriter.cc to store DeltaTheta:

entry->DeltaTheta = candidate->DeltaTheta

6. Regenerate Makefile and rebuild Delphes:
./configure
make -j 4

7. Add new module to the configuration file:

set ExecutionPath {
  ...
  GenJetFinder
  FastJetFinder

  MuonIsolation
  ...
}
  
module IsolationLEP MuonIsolation {

  set CandidateInputArray MuonMomentumSmearing/muons
  set JetInputArray FastJetFinder/jets

  set OutputArray muons

  set DeltaThetaMin 0.0
}

8. Add new ROOT tree branch to the configuration module:

module TreeWriter TreeWriter {
  ...
  add Branch MuonMomentumSmearing/muons Muon Muon
  add Branch MuonIsolation/muons MuonIso Muon
  ...
}

9. Now you can rerun Delphes:

  ./DelphesSTDHEP examples/delphes_card_FCC_basic.tcl step_3.root /afs/cern.ch/user/s/selvaggi/public/delphes_tuto/ee_zh_mmbb.hep


10. Have a look at the newly stored isolation variable:

root -l step_3.root
Delphes->Draw("MuonIso.DeltaTheta", "MuonIso.PT > 20 && MuonIso.DeltaTheta < 3.14")
gPad->SetLogy()

***************************************************************************************

Here are the new module source files:

//-----------------------------------------------------------------------------
// modules/IsolationLEP.h
//-----------------------------------------------------------------------------

#ifndef IsolationLEP_h
#define IsolationLEP_h

#include "classes/DelphesModule.h"

class TObjArray;

class ExRootFilter;

class IsolationLEP: public DelphesModule
{
public:

  IsolationLEP();
  ~IsolationLEP();

  void Init();
  void Process();
  void Finish();

private:

  Double_t fDeltaThetaMin;

  TIterator *fItCandidateInputArray; //!

  TIterator *fItJetInputArray; //!

  const TObjArray *fCandidateInputArray; //!

  const TObjArray *fJetInputArray; //!

  TObjArray *fOutputArray; //!

  ClassDef(IsolationLEP, 1)
};

#endif

//-----------------------------------------------------------------------------
// modules/IsolationLEP.cc
//-----------------------------------------------------------------------------

#include "modules/IsolationLEP.h"

#include "classes/DelphesClasses.h"
#include "classes/DelphesFactory.h"
#include "classes/DelphesFormula.h"

#include "ExRootAnalysis/ExRootResult.h"
#include "ExRootAnalysis/ExRootFilter.h"

#include "TMath.h"
#include "TString.h"
#include "TFormula.h"
#include "TRandom3.h"
#include "TObjArray.h"
#include "TDatabasePDG.h"
#include "TLorentzVector.h"
#include "TVector3.h"

#include <algorithm>
#include <stdexcept>
#include <iostream>
#include <sstream>

using namespace std;


//------------------------------------------------------------------------------

IsolationLEP::IsolationLEP() :
  fItCandidateInputArray(0), fItJetInputArray(0)
{
}

//------------------------------------------------------------------------------

IsolationLEP::~IsolationLEP()
{
}

//------------------------------------------------------------------------------

void IsolationLEP::Init()
{
  fDeltaThetaMin = GetDouble("DeltaThetaMin", 20);

  // import input array(s)

  fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "MuonMomentumSmearing/muons"));
  fItCandidateInputArray = fCandidateInputArray->MakeIterator();

  fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
  fItJetInputArray = fJetInputArray->MakeIterator();

  // create output array

  fOutputArray = ExportArray(GetString("OutputArray", "electrons"));
}

//------------------------------------------------------------------------------

void IsolationLEP::Finish()
{
  if(fItJetInputArray) delete fItJetInputArray;
  if(fItCandidateInputArray) delete fItCandidateInputArray;
}

//------------------------------------------------------------------------------

void IsolationLEP::Process()
{
  Candidate *jet, *candidate, *constituent;
  Double_t angle, deltaTheta = 4.0;
  //Int_t counter;
  Int_t nConstituents;

  // loop over all input leptons

  fItCandidateInputArray->Reset();
  while((candidate = static_cast<Candidate*>(fItCandidateInputArray->Next())))
  {
    const TLorentzVector &candidateMomentum = candidate->Momentum;
    const TVector3 &candidateBoost = candidateMomentum.BoostVector();

    // loop over all input jets
    fItJetInputArray->Reset();

    while((jet = static_cast<Candidate*>(fItJetInputArray->Next())))
    {
      const TLorentzVector &jetMomentum = jet->Momentum;
      const TVector3 &jetBoost = jetMomentum.BoostVector();

      nConstituents = 0;

      // compute number of jet constituents
      TIter itConstituents(jet->GetCandidates());
      itConstituents.Reset();

      while((constituent = static_cast<Candidate*>(itConstituents.Next())))
      {
        if (constituent->Charge != 0) nConstituents ++;
      }
      angle = TMath::Abs(jetBoost.Angle(candidateBoost));

      if(angle < deltaTheta && nConstituents > 1) deltaTheta = angle;
    }

    candidate->DeltaTheta = deltaTheta;

    if(deltaTheta < fDeltaThetaMin) continue;

    fOutputArray->Add(candidate);
  }
}


//-----------------------------------------------------------------------------