Example 5: Adding new module

============================



In the Isolation module that comes with the Delphes release, the isolation variable is not stored as a member of the candidate class but simply calculated to check whether the object under considered is isolated or not. 



In this example you will learn how to write a new Isolation module that will store the isolation variable as member of the candidate class. In order for this variable to be stored in the final electron collection, you will also have to modify the Electron class, where the relevant electron variables are defined, and the TreeWriter class, which translates Candidate collections into specific final state objects (electrons here).





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



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



3. In modules/NewIsolation.cc store the relative isolation variable in candidate: just after "ratio = sum/candidateMomentum.Pt();" add the following line: 



   candidate -> RelIso = ratio; 



4. Add the RelIso variable to Candidate class and Electron class in classes/DelphesClasses.h and in classes/DelphesClasses.cc



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



6. Modify TreeWriter::ProcessElectrons() in modules/TreeWriter.cc to store RelIso:



entry->RelIso = candidate->RelIso;



7. Regenerate Makefile and rebuild Delphes:

./configure

make -j 4



8. Modify the configuration file examples/delphes_card_CMS.tcl :





module NewIsolation ElectronIsolation {





 set CandidateInputArray ElectronEfficiency/electrons

 set IsolationInputArray EFlowMerger/eflow



 set OutputArray electrons

	

 set DeltaRMax 0.5

	

 set PTMin 0.5

	

 set PTRatioMax 10000	

}





9. Add new ROOT tree branch to the TreeWriter section of the conf. file:



module TreeWriter TreeWriter {

  ...

  add Branch ElectronIsolation/electrons ElectronIso Electron

  ...

}



10. Now you can rerun Delphes:



   ./DelphesSTDHEP examples/delphes_card_CMS.tcl step_3.root /projet/pth/delphes/Samples/z_ee.hep





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



root -l step_3.root

Delphes->Draw("ElectronIso.RelIso", "ElectronIso.PT > 20 && ElectronIso.RelIso < 1.0")

gPad->SetLogy()





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



Here are the new module source files:





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

// modules/NewIsolation.h

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



#ifndef NewIsolation_h

#define NewIsolation_h





#include "classes/DelphesModule.h"



class TObjArray;



class ExRootFilter;

class NewIsolationClassifier;



class NewIsolation: public DelphesModule

{

public:



  NewIsolation();

  ~NewIsolation();



  void Init();

  void Process();

  void Finish();



private:



  Double_t fDeltaRMax;



  Double_t fPTRatioMax;



  Double_t fPTSumMax;



  Bool_t fUsePTSum;



  NewIsolationClassifier *fClassifier; //!



  ExRootFilter *fFilter;



  TIterator *fItIsolationInputArray; //!



  TIterator *fItCandidateInputArray; //!



  TIterator *fItRhoInputArray; //!



  const TObjArray *fIsolationInputArray; //!



  const TObjArray *fCandidateInputArray; //!



  const TObjArray *fRhoInputArray; //!



  TObjArray *fOutputArray; //!



  ClassDef(NewIsolation, 1)

};



#endif



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

// modules/NewIsolation.cc

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



#include "modules/NewIsolation.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 NewIsolationClassifier : public ExRootClassifier

{

public:



  NewIsolationClassifier() {}



  Int_t GetCategory(TObject *object);



  Double_t fPTMin;

};



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



Int_t NewIsolationClassifier::GetCategory(TObject *object)

{

  Candidate *track = static_cast<Candidate*>(object);

  const TLorentzVector &momentum = track->Momentum;



  if(momentum.Pt() < fPTMin) return -1;



  return 0;

}



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



NewIsolation::NewIsolation() :

  fClassifier(0), fFilter(0),

  fItIsolationInputArray(0), fItCandidateInputArray(0),

  fItRhoInputArray(0)

{

  fClassifier = new NewIsolationClassifier;

}



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



NewIsolation::~NewIsolation()

{

}



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



void NewIsolation::Init()

{

  const char *rhoInputArrayName;



  fDeltaRMax = GetDouble("DeltaRMax", 0.5);



  fPTRatioMax = GetDouble("PTRatioMax", 0.1);



  fPTSumMax = GetDouble("PTSumMax", 5.0);



  fUsePTSum = GetBool("UsePTSum", false);



  fClassifier->fPTMin = GetDouble("PTMin", 0.5);



  // import input array(s)



  fIsolationInputArray = ImportArray(GetString("IsolationInputArray", "Delphes/partons"));

  fItIsolationInputArray = fIsolationInputArray->MakeIterator();



  fFilter = new ExRootFilter(fIsolationInputArray);



  fCandidateInputArray = ImportArray(GetString("CandidateInputArray", "Calorimeter/electrons"));

  fItCandidateInputArray = fCandidateInputArray->MakeIterator();



  rhoInputArrayName = GetString("RhoInputArray", "");

  if(rhoInputArrayName[0] != '\0')

  {

    fRhoInputArray = ImportArray(rhoInputArrayName);

    fItRhoInputArray = fRhoInputArray->MakeIterator();

  }

  else

  {

    fRhoInputArray = 0;

  }



  // create output array



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

}



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



void NewIsolation::Finish()

{

  if(fItRhoInputArray) delete fItRhoInputArray;

  if(fFilter) delete fFilter;

  if(fItCandidateInputArray) delete fItCandidateInputArray;

  if(fItIsolationInputArray) delete fItIsolationInputArray;

}



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



void NewIsolation::Process()

{

  Candidate *candidate, *isolation, *object;

  TObjArray *isolationArray;

  Double_t sum, ratio;

  Int_t counter;

  Double_t eta = 0.0;

  Double_t rho = 0.0;



  if(fRhoInputArray && fRhoInputArray->GetEntriesFast() > 0)

  {

    candidate = static_cast<Candidate*>(fRhoInputArray->At(0));

    rho = candidate->Momentum.Pt();

  }



  // select isolation objects

  fFilter->Reset();

  isolationArray = fFilter->GetSubArray(fClassifier, 0);



  if(isolationArray == 0) return;



  TIter itIsolationArray(isolationArray);



  // loop over all input jets

  fItCandidateInputArray->Reset();

  while((candidate = static_cast<Candidate*>(fItCandidateInputArray->Next())))

  {

    const TLorentzVector &candidateMomentum = candidate->Momentum;

    eta = TMath::Abs(candidateMomentum.Eta());



    // loop over all input tracks

    sum = 0.0;

    counter = 0;

    itIsolationArray.Reset();

    while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))

    {

      const TLorentzVector &isolationMomentum = isolation->Momentum;



      if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&

         !candidate->Overlaps(isolation))

      {

        sum += isolationMomentum.Pt();

        ++counter;

      }

    }



    // find rho

    rho = 0.0;

    if(fRhoInputArray)

    {

      fItRhoInputArray->Reset();

      while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))

      {

        if(eta >= object->Edges[0] && eta < object->Edges[1])

        {

          rho = object->Momentum.Pt();

        }

      }

    }



    // correct sum for pile-up contamination

    sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();



    ratio = sum/candidateMomentum.Pt();



    if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;



    fOutputArray->Add(candidate);

  }

}