My Analysis Tools

This page contains the documentation for my analysis code in UserCode/musella/Tools. The library can be downloaded and compiled as:

cd /src
export CVSROOT=:gserver:cmscvs.cern.ch:/cvs_server/repositories/CMSSW
cvs co -d musella/Tools UserCode/musella/Tools
scramv1 b -j 4 

here

from ROOT import gSystem 
gSystem.Load('${CMSSW_BASE}/lib/${SCRAM_ARCH}/libmusellaTools.so')
from ROOT import Event, Photons, Electrons, NtupleAnalyzer 

Generic NTuple reading

The most efficient way to save complex objects into a ROOT is to use dictionaries. In this way the serialization/deserialization process is left to the framework and the efforts on the user side are reduced to the creation of classes and dictionaries. In this case, however, all users need to share the dictionaries and care has to be taken to ensure backward compatibility of the produced trees.

For these and other reasons physicists often end up producing plain ntuples. Also, in big collaborations, ntuples are often produced by analysis groups and different groups usually use different ntuples. This is why, I ended up developing a set of classes that would allow to read different type of ntuples and present them through the same interface, leaving the downstream analysis code unchanged.

The solution I came up with is probably not very elegant nor extremely efficient, but it is very flexible and easy to use.

• The Event class manages the booking of branches and the looping over the tree. It class contains the Collection objects and allows to read only a subset of the events.
  Usage example:

// event loop 
TFile * file = TFile::Open("<some file>"); 
TTree * tree = file-&gt;Get("<some tree>"); 
Event ev(tree); 
for( ev.toBegin(); !ev.atEnd(); ev.next() ){ 
      Photons & photons = ev.photons(); 
      if( photons.size() > 0 ) { std::cout << photons[0].pt() << endl; 
} 

// Loop over 1/3 of the events. 
std::vector<int> keep_first(1,0); 
ev.prescale(3,keep_first); 
for( ev.toBegin(); !ev.atEnd(); ev.next() ){ 
      std::cout << (int) ev << std::endl; 
}

// Read the other 2/3 of the events. 
std::vector<int> keep_others; keep_others.push_back(1); keep_others.push_back(2); 
ev.prescale(3,keep_others); 
for( ev.toBegin(); !ev.atEnd(); ev.next() ){ 
      std::cout <<; (int) ev << std::endl; 
} 

• The Collection and Particle classes are used to access the content of the tree branches.
  The Collection class a aggregates several tree branches into a coherent structure and keeps track of selected candidates.
  The Particle class indexes the Collection class and presents the value of all branches as class methods.

Particle p = collection[0]; 
cout << p.energy() - c.energy(0) << std::endl; 

class Collection { 
public: 
      //! Define branch names 
      void branchName(const std::string & name, const std::string & val) { branchNames_[name] = val; } 
      //! Common prefix for branch names 
      void prefix(const std::string & pfx) { prefix_ = pfx; }; 
}

ev.photons().branchName("energy", "energy"); 
ev.photons().branchName("eta", "eta"); 
ev.photons().branchName("phi", "phi"); 

// Define buffers and getters 
class Collection { 
public: 
      COLL_GETTER_MULT(energy, Float_t, escale_ ); 
      COLL_GETTER(eta , Float_t ); 
      COLL_GETTER(phi , Float_t ); 
}

// Define particle getters 
class Particle { 
public: 
      PART_GETTER(energy, Float_t, c_ );
      PART_GETTER(eta, Float_t, c_ ); 
      PART_GETTER(phi, Float_t, c_ ); 
}

// Branch binding 
void Collection::getBranches(Event & ev) { 
      GET_BRANCH(energy); 
      GET_BRANCH(eta); 
      GET_BRANCH(phi); 
} 

• Physics objects classes inherit from the Collection and Particle classes:
  • Photons
  • Electrons
  • Jets
  • MET
  • Vertexes
  • Tracks

NTuple analyser

The analysis code is supposed to be contained in classes denominated NTuple analyzers. A base class called NTupleAnalyzer implements all the basic blocks needed to loop over events.

The class Implements the loop over datasets and events, histograms bookkeeping and basic event selection. It provides hook for the analysis code defining virtual methods to be overridden by inheriting classes.

• Analyzer definition. Show the code.... Hide the code

// Analyzer
class MyAnalyzer : public NtupleAnalyzer {
begin:
      void begin() {
            openOutput( "myHistoFile.root" );
            bookHistogram("bx_id","Bunch crossing id",0xdec,1,0xdec+1);
      }

      virtual void analyze(Event &e) {
            histogram("bx_id")->Fill(e.bx());
      }

      void end() {
           closeOutput();
      }
}

• Loop over events.

// n.events to analyze per dataset (-1 for all)
int maxEvent = 100;
// Event handling object
Event ev;
// Fill list of datasets
std::vector<Dataset *> myDatasets;
myDatasets.push_back( <dataset 1> );
myDatasets.push_back( <dataset 2> );
...

// Instantiate analyzer
MyAnalyzer myAnalyzer;

// Loop over datatsers
myAnalyzer.analyze(myDatasets,&ev,maxEvent);

Helper classes

A few helper classes can be used in association with the NtupleAnalyzer class, to build analyses.

• The AutoPlotter uses the ROOT TTreeFormula class to provide automatic plot generation. Show usage example.... Hide the code

class MyAnalyzer : public NtupleAnalyzer {
begin:
      AutoPlot myPlotter_;

      begin() {
             openOutput("myHistos.root");
             // define plots
             vector<string> plots;
             plots.push_back("bx:(0xdec:1:0xded):Bunch crossing; bx; events");
             plots.push_back("metEt:(200:0:200):MET; MET [GeV]; Events/GeV");
             myPlotter_.setPlots(plots);
             // allocate functors
             myPlotter_.makeFormulas();
             // add to the list of plotters
             plotters_.push_back(&myPlotter_);
      }

      void beginDataset() {
            // Book histograms for all categories
            myPlotter_.bookHistograms("");
            myPlotter_.bookHistograms("cat1");
            myPlotter_.bookHistograms("cat2");
      }

      virtual void analyze(Event &e) {
           float weight = unweigh_ ? 1. : e.getTree()->GetWeight();
           if( eventSelected ) {
                   // define which categories should be filled for this event
                   vector<string> categories(1,"");
                   if( cat1 ) categories.push_back("cat1");
                   if( cat2 ) categories.push_back("cat2");
                   // fill all plots
                   myPlotter_.fillPlots(weight, categories);
           }
      }

      void end() {
           closeOutput();
      }
}

• The PhotonSelector class can be used to define a set of selection criteria for photon objects. The class keeps track of which selection criteria are satisfied and allows to define groups of selection criteria and study N-M groups of cuts.

Sample weighing

The Dataset class is used to manage different dataset with different weights. The weight of a simulated sample is calculated through the formula where is the cross section of the simulated process, is the number of simulated events and the luminosity being simulated. The class inherits from the TChain class and provides weights trough the standard GetWeight() method from the TTree class, thus in a completely transparent way for the analysis code. The relevant methods of the class are shown below. Show the code.... Hide the code.

class Dataset : public TChain {
public:
        // Add file to dataset
	void AddFile(const char* fname, double xsec, TString processId="", double minScale=-1., double maxScale=1.e+99, 
		     Long64_t nentries = kBigNumber, const char* tname = "");

        // Set equiv. lumi for weight calculation
	void setRefLumi(double lumi) { refLumi_ = lumi; };
	
        // Define event scale expression (uses TTreeFormula)
	void eventScale(const char * formula);
	
        // Recalculate weights
	void reweight();
	
        // Get weight for current event
	virtual Double_t GetWeight() const;
}

Static web pages

Static web pages can be generated using the HtmlHelper class.

HtmlHelper helper("<output director>");

helper.body().add(new HtmlTag("h1")) << "Usage example of the HtmlHelper class";
helper.body().add(new HtmlTag("p")) << "This will generate a static page with one plot.";

TCanvas * canv  = new TCanvas();
// generate the plot
....
helper.body().add(new HtmlPlot(canv)).caption("This is the plot");

helper.dump();

Cut optimization

http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/musella/Tools/interface/Optimizer.h?view=markup

Template fitting

http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/musella/Tools/interface/Purity.h?view=markup

http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/UserCode/musella/Tools/interface/purity_fit.h?view=markup