cp -r /atlas/aad/workarea/Higgs/ttHAna2015/rootAna/Mig/NNLOReweighter .

cp -r /afs/cern.ch/user/a/aad/public/NNLOReweighter .

compile the code

Now we have all needed packages. Lets compile:

rc find_packages
rc compile

Run the code

The code is made to be very flexible and based on "JobOptions" written in C++ (live with it). The executable is called run (even if you can easily create your own). It includes other C++ files (with .h extension, arbitrary choice) that contain functions to configure the code. Run can read options also from commend line and from a configuration text file (default in share/config.txt). For now lets just run the default. The input file can be a single root file (should end with the .root extension) or a list of rootfiles in a text file (should end with a .txt). Wildcards are only supported when used in the text file following what is accepted by the ROOT TChain class.

echo /quark4/ttHbb/TopAnaOutput/00-02-03/mc/1l/group.phys-higgs.mc15_13TeV.343366.aMCNloP8_ttH125_slep.e4706s2726r7772p2952.TTHbbL022.1l.v2_out.root/group.phys-higgs.11138029._000001.out.root > input.txt
## or for people at cern: echo  /afs/cern.ch/work/a/aad/public/group.phys-higgs.11138029._000001.out.root  > input.txt
## run input.txt
## equivalent to run /quark4/ttHbb/TopAnaOutput/00-02-03/mc/1l/group.phys-higgs.mc15_13TeV.343366.aMCNloP8_ttH125_slep.e4706s2726r7772p2952.TTHbbL022.1l.v2_out.root/group.phys-higgs.11138029._000001.out.root
## lets limit the number of events to run faster
run input.txt maxEvents=10000

This will create an output file called tth_ana_output.root. You can specify the output file name by giving it as a second argument after the output. The output file name should end with .root extension (otherwise the code will switch to another output format assuming outputs will be written in independent root files and not in a single one, not covered in this tutorial).

run input.txt  output.root maxEvents=10000

Running your own stuff

share/run.cxx contains the main function. It loops over the list of systematics (configurable, for now we are only using nominal) and create a sequence of algorithms to be executed by the ApplicationMgr. You basically have to add your algorithm to a sequence to be executed. It will have access to I/O and the DataStore through the AlgoBase class. To make things easier the flag analysisName= is created and used inside run.cxx. It calls a predefined function in a separated JO file for a specific study. Default is "ttHDefaultAna". Some of you have already their JO with their name. For this tutorial I created analysisName=CPPMTuTo which calls the function runTuto in util/CPPMTuTo.h. The code to do this rudimentary:

#include "CPPMTuTo.h"
...
if("CPPMTuTo"==config.analysisName){
   runTuto(theapp, config, runOption);
}

So lets run with CPPMTuTo:

run input.txt  output.root maxEvents=10000 analysisName=CPPMTuTo

std::string localName = "cppmTuto";
  ///create our sequence and add it to application manager.
  AlgSequence* basic_seq=new AlgSequence("AlgSequence_"+localName);
  theapp.addAlgo(basic_seq);

  TuToAlgo* tutoalgo = new TuToAlgo("TuToAlgo_"+localName);
  tutoalgo->m_jetcontainername="selected_jets";
  basic_seq->addAlgo(tutoalgo);

Now lets create a new algo to the sequence. First we have to create it. Lets call it TuToAlgo. So you have to create TuToAlgo.h class in TTHbbNtupleAnalysis and TuToAlgo.cxx in Root. TuToAlgo.h will contain (have to inherit from AlgoBase):

#ifndef TuToAlgo_HH
#define TuToAlgo_HH

#include "TTHbbNtupleAnalysis/AlgoBase.h"

class TuToAlgo : public AlgoBase{

 public:

  explicit TuToAlgo(std::string);
  virtual ~TuToAlgo();

  virtual bool initialize();
  virtual bool execute();
  virtual bool finalize();
  
  std::string m_jetcontainername;

 private:
  double m_eventWeight;
  TH1* h_jet_pt;
};
#endif

And TuToAlgo.cxx contains:

#include "TTHbbNtupleAnalysis/TuToAlgo.h"

TuToAlgo::TuToAlgo(std::string name) : AlgoBase(name){
  m_jetcontainername="";
}

TuToAlgo::~TuToAlgo(){
  ////
}


bool TuToAlgo::initialize(){

  std::cout << m_name << " initialize " << endl;
  m_histoSvc->registerHisto(h_jet_pt,"h_jet_pt",60,25,250);
  m_histoSvc->sumw2();
  m_cuts->setEventWeight(m_eventWeight);

  return true;

}

bool TuToAlgo::execute(){

  m_cuts->ocount("in execute "+m_name,0);

  MetaData *meta = m_dataStore->getMetaData();
  if(!meta){
    std::cout << "WARNING " << m_name << " meta data not found " << std::endl;
    return false;
  }

  WeightData *weightData = m_dataStore->getWeightData();
  if(!weightData){
    std::cout << "WARNING " << m_name << " weight data not found " << std::endl;
    return false;
  }
  m_eventWeight = weightData->getEventWeight("lumi*weight_mc"); ///default "lumi*weight_mc*weight_leptonSF" set in run.cxx, will make default configurable from config.txt

  JetContainer* jets = m_dataStore->getJetContainer(m_jetcontainername);
  if(!jets){
    std::cout << "WARNING " << m_name << " JetContainer " << m_jetcontainername << "  not found " << std::endl;
    return false;
  }

  for(auto jet: *jets){
    h_jet_pt->Fill(jet->Pt()/1000., m_eventWeight);
  }
 
  return true;
}

bool TuToAlgo::finalize(){

  std::cout << m_name << " finalize " << endl;
  return true;

}

add in utils/run.cxx at the beginning with all other algorithm includes:

#include "TTHbbNtupleAnalysis/TuToAlgo.h"

Now compile go to the run directory that you created and run the code again. You should have your histo in the output file. Browse the output to find it. Notice that it is under a directory with your algo name.

adding histos in different regions

Now lets add several histos to the same algo with different region cuts. First we have to make the algo inherit from RegionInterface. Add to TuToAlgo.h:

#include "TTHbbMVAInterfaces/RegionInterface.h"
#include "MVAVariables/JetOrderingTool.h"

TH1* h_jetpt1_pt[20]; /// max allowed regions is 20

Then in TuToAlgo.cxx: Add "RegionInterface()," to the constructor. Before "m_histoSvc->sumw2();" add the following to the Initialize function:

for(auto r : m_regions){
    m_histoSvc->registerHisto(h_jetpt1_pt[r.index],"h_jetpt1_pt_"+r.name,60,25,300);
  }

RecoEvent* eventvar = m_dataStore->getRecoEvent(m_eventHolderKey);
  if(!eventvar){
    std::cout << "WARNING " << m_name << " RecoEvent " << m_eventHolderKey << "  not found " << std::endl;
    return false;
  }
///.....
 m_jetOrderingTool.sortPt(*jets);

  vector<int> regions;
  getRegions(regions, eventvar->nJets, eventvar->nBJets, eventvar->leptonType, eventvar);

  for(auto region : regions){
    if(jets->size()>0)h_jetpt1_pt[region]->Fill((*jets)[0]->p4().Pt()/1000.,m_eventWeight);
  }

std::string localEventName = "cppmEvent"; /// the name of the event in the store
  std::string defaultBtagID = "isbtagged_85"; /// the default btag cut (e.g. nBTags in the event will be @85%,nBTags_85, nBTags_77 decorations ... are also available).

  //// create an event holder (TTHbb::RecoEvent class) that will hold all event based variables including choosed reco combinations
  EventHolderCreator* evtHolder = new EventHolderCreator("EventHolderCreator_"+localName);
  evtHolder->m_inputJets = "selected_jets";
  evtHolder->m_inputLeptons = "selected_leptons";
  evtHolder->m_inputNeutrinos = "basic_neutrinos";
  evtHolder->m_eventHolderKey=localEventName;
  evtHolder->m_btagIdName=defaultBtagID; 
  basic_seq->addAlgo(evtHolder);
  /// or use prepareEventHolder(basic_seq, localEventName, defaultBtagID); /// defined in util/JobOptions/EventHolderHelperFunc.h

  /// this will add few global variable to the event holder, needed to cut on later (like nBTags_77)
  addGlobalVars(basic_seq, localEventName); 
  addObjectVars(basic_seq, localEventName);
  addMVAVars(basic_seq, localEventName, config); /// do not forget to uncomment config in the funtion declaration

Now lets add regions to the algo. There a couple of ways to do it. We can create the region and add it explicitly like this:

///configure region with cuts. Can use decoration into the EventHolder or members like nJets and nBjets (@ default cut that the event holder was created with)
  RegionInterface::RegionDefinition r_6jupr;
  r_6jupr.name="6jUPR";
  r_6jupr.intLowerCut("nBTags_77",4);
  r_6jupr.intUpperCut("nBTags_77",10000);
  r_6jupr.nJetsMin=6;
  r_6jupr.nJetsMax=10000;

  RegionInterface::RegionDefinition r_6jsig =  getRegion("6jsig", "SemilepRegions_signalRegions", 601);

  tutoalgo->addRegion(r_6jupr);
  tutoalgo->addRegion(r_6jsig);
  tutoalgo->addRegion("6ji3be",6,10000,3,3,0); /// will create region with name, njetsmin, njetsmax, nbjetsmin(@default), nbjetsmax(@default), lepton flavor (0 for no cuts).
  tutoalgo->m_eventHolderKey=localEventName;

and now compile and run and see what you get.

-- GeorgesAad - 2017-04-23