AnalysisBase.h

Go to the documentation of this file.

#ifndef _ANALYSISBASE
#define _ANALYSISBASE

#include <iostream>
#include <fstream>
#include <stdio.h>
#include <map>
#include <math.h>

#include <TCanvas.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TH1.h>
#include <TObject.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TSystem.h>

#include "classes/DelphesClasses.h"
#include "external/ExRootAnalysis/ExRootTreeReader.h"
#include "external/ExRootAnalysis/ExRootResult.h"

#include "ETMiss.h"
#include "FinalStateObject.h"
#include "Units.h"

#include "mt2_bisect.h"
#include "mctlib.h"
#include "mt2bl_bisect.h"


class AnalysisBase {
 public:

    AnalysisBase(std::string inFile, std::string outFol, std::string outPre, double xs, double xserr, std::map<std::string, std::string> branches, std::map<std::string, std::vector<int> > flags);
    ~AnalysisBase(); 

    void loopOverEvents();

 protected:   


    virtual void initialize() {};
    virtual void analyze() {}; 
    virtual void finalize() {}; 

    std::vector<Electron*> electrons; 
    std::vector<Electron*> electronsLoose; 
    std::vector<Electron*> electronsMedium; 
    std::vector<Electron*> electronsTight; 
    std::vector<Muon*> muons; 
    std::vector<Muon*> muonsCombinedPlus; 
    std::vector<Muon*> muonsCombined; 
    std::vector<Jet*> jets; 
    std::vector<Jet*> jets2; 
    std::vector<Photon*> photons; 
    std::vector<Track*> tracks; 
    std::vector<Tower*> towers; 
    ETMiss* missingET; 



    template <class T>
    std::vector<T*> filterPhaseSpace(std::vector<T*> unfiltered, double pTmin = 0., double etamin = -100, double etamax = 100, bool exclude_overlap = false) {
    std::vector<T*> filtered;
    for (int i = 0; i < unfiltered.size(); i++) {
        T* cand = unfiltered[i];
        if((cand->PT > pTmin) && (cand->Eta > etamin) && (cand->Eta < etamax)) {
              if(!exclude_overlap)
        filtered.push_back(cand);
              else if( (fabs(cand->Eta) < 1.37) || (fabs(cand->Eta) > 1.52) )
                filtered.push_back(cand);
            }
    }
    return filtered;
    }       


    template <class X, class Y>
    std::vector<X*> overlapRemoval(std::vector<X*> candidates, std::vector<Y*> neighbours, double dR) {
      // If neighbours are empty, return candidates
      if(neighbours.size() == 0)
        return candidates;
      std::vector<X*> passed_candidates;
      // Loop over candidates
      for(int i = 0; i < candidates.size(); i++) {
        bool overlap = false;
        // If a neighbour is too close, declare overlap, break and don't save candidate
        for(int j = 0; j < neighbours.size(); j++) {
          if (candidates[i]->P4().DeltaR(neighbours[j]->P4()) < dR) {
            overlap = true;
            break;
          }
        }
        if (!overlap)
          passed_candidates.push_back(candidates[i]);
      }
      return passed_candidates;
    }


    template <class X>
    std::vector<X*> overlapRemoval(std::vector<X*> candidates, double dR) {
      // Same as above for the special case that candidates = neighbours. In that case, the removal 
      // can be formulated more effectively as the sum only has to run half as many times
      if(candidates.size() == 0)
        return candidates;
      std::vector<X*> passed_candidates;
      // Loop over candidates
      for(int i = 0; i < candidates.size(); i++) {
        bool overlap = false;
        // If one of the other, still untested, candidates is too close: remove
        // Since the list is order w.r.t pt, this will always remove the softer object
        for(int j = 0; j < i; j++) {
          if (candidates[i]->P4().DeltaR(candidates[j]->P4()) < dR) {
            overlap = true;
            break;
          }
        }
        if (!overlap)
          passed_candidates.push_back(candidates[i]);
      }
      return passed_candidates;
    }
    

    std::vector<Electron*> filterIsolation(std::vector<Electron*> unfiltered, std::vector<int> relative_flags = std::vector<int>()) {
      // Translate the relative isolation number of the analysis in the absolute number within all analyses
        std::vector<int> absolute_flags;
        for(int i = 0; i < relative_flags.size(); i++) {
          if (relative_flags[i] >= electronIsolationFlags.size()) {
            std::cerr << "Error: There is no electron isolation " << relative_flags[i] << std::endl;
            std::cerr << "Exiting... "<< std::endl;
            exit(1);
          }
          absolute_flags.push_back(electronIsolationFlags[relative_flags[i]]);
        }
        // if no flags are given, use all
        if (absolute_flags.size() == 0)          
          absolute_flags = electronIsolationFlags;
          
        return filterFlags(unfiltered, "isolation", absolute_flags);
    }
    
    std::vector<Electron*> filterIsolation(std::vector<Electron*> unfiltered, int relative_flag) {
      std::vector<int> absolute_flags;
      if (relative_flag >= electronIsolationFlags.size()) {
          std::cerr << "Error: There is no electron isolation " << relative_flag << std::endl;
          std::cerr << "Exiting... "<< std::endl;
          exit(1);
        }
        absolute_flags.push_back(electronIsolationFlags[relative_flag]);        
        return filterFlags(unfiltered, "isolation", absolute_flags);
    }    
    
    std::vector<Muon*> filterIsolation(std::vector<Muon*> unfiltered, std::vector<int> relative_flags = std::vector<int>()) {
      // Translate the relative isolation number of the analysis in the absolute number within all analyses
        std::vector<int> absolute_flags;
        for(int i = 0; i < relative_flags.size(); i++) {
          if (relative_flags[i] >= muonIsolationFlags.size()) {
            std::cerr << "Error: There is no muon isolation " << relative_flags[i] << std::endl;
            std::cerr << "Exiting... "<< std::endl;
            exit(1);
          }
          absolute_flags.push_back(muonIsolationFlags[relative_flags[i]]);
        }
        // if no flags are given, use all
        if (absolute_flags.size() == 0)          
          absolute_flags = muonIsolationFlags;
          
        return filterFlags(unfiltered, "isolation", absolute_flags);
    }
    
    std::vector<Muon*> filterIsolation(std::vector<Muon*> unfiltered, int relative_flag) {
        std::vector<int> absolute_flags;
        if (relative_flag >= muonIsolationFlags.size()) {
            std::cerr << "Error: There is no muon isolation " << relative_flag << std::endl;
            std::cerr << "Exiting... "<< std::endl;
            exit(1);
          }
          absolute_flags.push_back(muonIsolationFlags[relative_flag]);
        
        
        return filterFlags(unfiltered, "isolation", absolute_flags);
    }    
    
    std::vector<Photon*> filterIsolation(std::vector<Photon*> unfiltered, std::vector<int> relative_flags = std::vector<int>()) {
      // Translate the relative isolation number of the analysis in the absolute number within all analyses
        std::vector<int> absolute_flags;
        for(int i = 0; i < relative_flags.size(); i++) {
          if (relative_flags[i] >= photonIsolationFlags.size()) {
            std::cerr << "Error: There is no photon isolation " << relative_flags[i] << std::endl;
            std::cerr << "Exiting... "<< std::endl;
            exit(1);
          }
          absolute_flags.push_back(photonIsolationFlags[relative_flags[i]]);
        }
        // if no flags are given, use all
        if (absolute_flags.size() == 0)          
          absolute_flags = photonIsolationFlags;
          
        return filterFlags(unfiltered, "isolation", absolute_flags);
    }
    
    std::vector<Photon*> filterIsolation(std::vector<Photon*> unfiltered, int relative_flag) {
        std::vector<int> absolute_flags;
        if (relative_flag >= photonIsolationFlags.size()) {
            std::cerr << "Error: There is no photon isolation " << relative_flag << std::endl;
            std::cerr << "Exiting... "<< std::endl;
            exit(1);
          }
          absolute_flags.push_back(photonIsolationFlags[relative_flag]);
        
        
        return filterFlags(unfiltered, "isolation", absolute_flags);
    }    
    

    bool checkTauTag(Jet* candidate, std::string efficiency);
    

    bool checkBTag(Jet* candidate, int relative_flag = 0, std::string option = "");


    double mT(const TLorentzVector & vis, const TLorentzVector & invis);
    

    double mT2(const TLorentzVector & vis1, const TLorentzVector & vis2, double m_inv, const TLorentzVector & invis = TLorentzVector(0., 0., 0., 0.));


    double mCT(const TLorentzVector & v1, const TLorentzVector & v2);
    

    double mCTcorr(const TLorentzVector & v1, const TLorentzVector & v2, const TLorentzVector & vds, const TLorentzVector & invis, const double ecm = 8000.0, const double mxlo = 0.0);
    

    double mCTy(const TLorentzVector & v1, const TLorentzVector & v2, const TLorentzVector & vds, const TLorentzVector & invis);        
    

    double mT2_bl(const TLorentzVector & pl_in, const TLorentzVector & pb1_in, const TLorentzVector & pb2_in, const TLorentzVector & invis = TLorentzVector(0., 0., 0., 0.));    


    double alphaT(const std::vector<Jet*> & jets,  const double thresh_ET = 0.);
    

    std::vector<double> razor(const std::vector<TLorentzVector> & obj, const TLorentzVector & invis = TLorentzVector(0., 0., 0., 0.));

    void bookSignalRegions(std::string listOfRegions);  
    void bookControlRegions(std::string listOfRegions); 
    void bookCutflowRegions(std::string listOfRegions); 
    

    inline void countSignalEvent(std::string region) {      
      signalRegions[region] += weight;
      signalRegions2[region] += weight*weight;
    }
    inline void countControlEvent(std::string region) {      
      controlRegions[region] += weight;
      controlRegions2[region] += weight*weight;
    }
    inline void countCutflowEvent(std::string region) {
      cutflowRegions[region] += weight;
      cutflowRegions2[region] += weight*weight;
    }    

    // File streams for generalized output files
    std::vector<std::ofstream*> fStreams;  
    std::vector<std::string> fNames; 


    int bookFile(std::string name, bool noheader = false); 
        
    


    inline void setInformation(std::string s) {
    information = s;
    };
    

    inline void setLuminosity(double l) {
      luminosity = l;
    };
    

    inline void setAnalysisName(std::string name) {
      analysis = name;
    };
        

    inline double normalize(double x) {
    return x*xsect*luminosity/sumOfWeights;
    };
    
    ExRootResult *result;
    

    void ignore(std::string ignore_what); 
    double weight; 
    
 private:
    // Overall reader of the ROOT tree
    ExRootTreeReader *treeReader;   
    TChain *chain;
    
    // Objects for all the branches in the ROOT file
    TClonesArray *branchElectron;
    TClonesArray *branchMuon;
    TClonesArray *branchJet;
    TClonesArray *branchJet2;
    TClonesArray *branchPhoton;
    TClonesArray *branchMissingET;
    TClonesArray *branchEvent;
    TClonesArray *branchTrack;
    TClonesArray *branchTower;    
        
    // Information on files
    std::string outputFolder;
    std::string outputPrefix;
    std::string inputFile;

    // Information about the analysis to be printed at the start of each output file
    std::string analysis;
    std::string information;
    
    // Global parameters
    Long64_t nEvents;
    double sumOfWeights;
    double sumOfWeights2;
    double xsect;
    double xsecterr;
    double luminosity;
        
    // Sums up weights (and weights^2) that fall into control, signal or cutflow regions 
    std::map<std::string, double> controlRegions;
    std::map<std::string, double> signalRegions;    
    std::map<std::string, double> cutflowRegions;
    std::map<std::string, double> controlRegions2;
    std::map<std::string, double> signalRegions2;    
    std::map<std::string, double> cutflowRegions2;
    
    // Absolute flag values for the isolation criteria of the individual analysis
    // [e.g. if there are 2 analyses with 2 isolation criteria each, analysis 1 gets (1, 2)
    //  whereas analysis 2 gets (3, 4) ]
    std::vector<int> electronIsolationFlags;
    std::vector<int> muonIsolationFlags;
    std::vector<int> photonIsolationFlags;
    std::vector<int> jetBTagFlags;
    std::vector<int> jet2BTagFlags;    
    
    // If these flags are set, the corresponding containers won't be set up (to save computing time)
    bool ignoreElectrons;
    bool ignoreElectrons_LooseIsolation;
    bool ignoreElectrons_MediumEfficiency;
    bool ignoreElectrons_TightEfficiency;
    bool ignoreMuons;
    bool ignoreMuons_LooseIsolation;
    bool ignoreMuons_CombinedPlusEfficiency;
    bool ignoreMuons_CombinedEfficiency;
    bool ignorePhotons;
    bool ignorePhotons_LooseIsolation;
    bool ignoreTracks;
    bool ignoreTowers;

    // A given set of objects (Electrons, Jets, ...) can be filtered w.r.t given flag values
    template <class T>
    std::vector<T*> filterFlags(std::vector<T*> unfiltered, std::string whichFlag, std::vector<int> flags) {
      std::vector<T*> filtered;
      for (int i = 0; i < unfiltered.size(); i++) {
        T* cand = unfiltered[i];
        // First, decode the member's isoflag into a vector of valid flags
        int code = 0;
        if(whichFlag == "isolation") {
          code = cand->IsolationFlags;
        }        
        else if(whichFlag == "efficiency")
          code = cand->EfficiencyFlags;
        
        std::vector<int> candidates_flags = code_to_flags(code);
        
        // Then, check if iso_flags are all within the member's flags        
        bool passes = true;
        for (int iso = 0; iso < flags.size(); iso++) {
          if(std::find(candidates_flags.begin(), candidates_flags.end(), flags[iso]) == candidates_flags.end()) {
            passes = false;
            break;
          }
        }
        
        // Only if memebr passed it should be saved in filtered list
        if(passes)
          filtered.push_back(cand);
      }
      return filtered;
    }
    
    // The flagnumber, read in binary, tells which flags have been set to true
    std::vector<int> code_to_flags(int code) {
      std::vector<int> flags;
      int flag = 0;
      while(code > 0) {
        if(code % 2 == 1)
          flags.push_back(flag);
          code -= code % 2;
        code /= 2;
        flag++;
      }
      return flags;
    }
    
    // Used by alphaT code
    struct fabs_less { 
      bool operator()(const double x, const double y) const { 
        return fabs(x) < fabs(y); 
      } 
    };

};


#endif