mstein/triggerStudy/triggerStudy.C

/////////////////////////////////////////////////////////
//                  Trigger Study
/////////////////////////////////////////////////////////
//to run, execute: root run.C
//
//In this file, the implementation of the T&P method and 
//the evaluation of the trigger efficiency is performed. 
//some documentation can be found here:
//https://twiki.cern.ch/twiki/bin/viewauth/CMSsandbox/MatthiasSteinTriggerStudy


#include <iostream>
#include <fstream>
#include <map>
#include <vector>

#include <Math/LorentzVector.h>

#include "TString.h"
#include "TH1.h"
#include "TH1D.h"
#include "TTree.h"
#include "TFile.h"

#include "plotMaker.h"
#include "tools.h"
//--------------------------------------------------- declare some special types
#ifdef __MAKECINT__
#pragma link C++ class pair<string,bool>+;
#pragma link C++ class pair<string,string>+;
#pragma link C++ class map<string,bool>+;
#pragma link C++ class map<string,string>+;
#pragma link C++ class ROOT::Math::PtEtaPhiM4D<float>+;
#pragma link C++ class ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<float> >+;
#pragma link C++ class vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<float> > >+;
#endif

typedef ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<float> > LorentzM;

//--------------------------------------------------- Structure to store the whole event information in order to pass it to functions
struct eventInfo{
  vector<string> *DESYtriggerElMatchedTriggerFilter;
  vector<string> *DESYtriggerMuMatchedTriggerFilter;  
  map<string, string> *DESYtriggerNameMap;
  vector<int> *electronCharge;
  vector<LorentzM> *electronP4;
  int Event;
  double HT;
  vector<LorentzM> *jetP4;
  double MET;
  vector<int> *muonCharge;
  vector<LorentzM> *muonP4;
  map<string, int> *prescaled;
  int Run;
  int PUInter;
  double PUWeight;
  map<string, bool> *triggered;
  vector<vector<TString> > triggerFilterEl;
  vector<vector<TString> > triggerFilterMu;
  double Weight;
};

//--------------------------------------------------- functions
// void init(TString inputFileName);
void handleTriggerMatchingInfo(eventInfo& evt);
void initBranches(TTree *tree, eventInfo &event);
void investigateTriggerFilterMu(eventInfo& evt, map<TString, plotMaker> &plot, TString filterName);
void Loop(TString inputFileName, TString outputFileName);
void makeEfficiencyPlots(map<TString, plotMaker> &plot);
void plotHistograms(eventInfo& evt, map<TString, plotMaker> &plot);
void printElectronInfo(eventInfo& evt, ostream& os=cout);
void printEventInfo(eventInfo& evt, ostream& os=cout);
void printMuonInfo(eventInfo& evt, ostream& os=cout);
void setEfficiencyErrors(TH1D* &eff, TH1D* &h_passed, TH1D* &h_failed);

//--------------------------------------------------- global Variables
vector<string> triggerFilterToCheckMu;
vector<string> triggerFilterToCheckEl;
//----------------------------------------
const Int_t    binPt = 25;
const Double_t lowPt = 0.;
const Double_t upPt  = 250.;
//----------------------------------------
const Int_t    binEta = 50;
const Double_t lowEta = -2.5;
const Double_t upEta  = 2.5;
//----------------------------------------
const Int_t    binPU = 64;
const Double_t lowPU = 0.;
const Double_t upPU  = 64.;


//================================================================================= triggerStudy
//main function which will be called
//here, it is possible to run over multiple files and to estimate trigger efficiencies for different samples
void triggerStudy(){
  
  //--------------------------------------------------- define the muon TriggerFilters to be checked
  triggerFilterToCheckMu.push_back("hltL1sMu16Eta2p1");
  triggerFilterToCheckMu.push_back("hltL2fL1sMu16Eta2p1L1f0L2Filtered16Q");
  triggerFilterToCheckMu.push_back("hltL3crIsoL1sMu16Eta2p1L1f0L2f16QL3f20L3crIsoFiltered10");
  triggerFilterToCheckMu.push_back("hltIsoMu202p1TriCentralPFJet30MuCleaned");
  triggerFilterToCheckMu.push_back("hltIsoMu202p1TriCentralPFJet303020MuCleaned");
  triggerFilterToCheckMu.push_back("hltDiMuonGlb22Trk8DzFiltered0p2");  
  triggerFilterToCheckMu.push_back("hltL3crIsoL1sMu14Eta2p1L1f0L2f16QL3f20L3crIsoRhoFiltered0p15");  
  //--------------------------------------------------- define the electron TriggerFilters to be checked
//   triggerFilterToCheckEl.push_back("hltL1sL1SingleEG20orL1SingleEG22");  
//   triggerFilterToCheckEl.push_back("hltEle25CaloIdVTTrkIdTCaloIsoTTrkIsoTTrackIsoFilter");  
//   triggerFilterToCheckEl.push_back("hltEle25CaloIdVTCaloIsoTTrkIdTTrkIsoTTriCentralPFJet30EleCleaned");  
//   triggerFilterToCheckEl.push_back("hltEle25CaloIdVTCaloIsoTTrkIdTTrkIsoTTriCentralPFNoPUJet30EleCleaned");  
//   triggerFilterToCheckEl.push_back("hltEle25CaloIdVTCaloIsoTTrkIdTTrkIsoTTriCentralPFNoPUJet303020EleCleaned");  
  
  //--------------------------------------------------- define the sample(s) to loop over and the output file
  //these are small test files of ONE single file af a sample
//   Loop("/scratch/hh/current/cms/user/mstein/output/TTJets_mu_tree.root", "out_TTJets_mu_tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/TTJets_el_tree.root", "out_TTJets_el_tree.root");
//  Loop("/scratch/hh/current/cms/user/mstein/output/DYJets_mu_tree.root", "out_DYJets_mu_tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/DYJets_el_tree.root", "out_DYJets_el_tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/SingleMu_mu_tree.root", "out_SingleMu_mu_tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/SingleElectron_el_tree.root", "out_SingleElectron_el_tree.root");
  //these are the whole samples
//   Loop("/scratch/hh/current/cms/user/mstein/output/TTJets/SUMMER12/TTJets_SUMMER12_TrigStudy-mu_noTail_Tree.root", "out_TTJets_SUMMER12_TrigStudy-mu_noTail_Tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/TTJets/SUMMER12/TTJets_SUMMER12_TrigStudy-el_noTail_Tree.root", "out_TTJets_SUMMER12_TrigStudy-el_noTail_Tree.root");
   Loop("/scratch/hh/current/cms/user/mstein/output/DYJets/M-50/DYJets_M-50_TrigStudy-mu_noTail_Tree.root", "out_DYJets_M-50_TrigStudy-mu_noTail_Tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/DYJets/M-50/DYJets_M-50_TrigStudy-el_noTail_Tree.root", "out_DYJets_M-50_TrigStudy-el_noTail_Tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/SingleElectron/Run2012B-13Jul2012-v1/SingleElectron_Run2012B-13Jul2012-v1_TrigStudy-el_noTail_Tree.root", "out_SingleElectron_Run2012B-13Jul2012-v1_TrigStudy-el_noTail_Tree.root");
//   Loop("/scratch/hh/current/cms/user/mstein/output/SingleMu/Run2012B-13Jul2012-v1/SingleMu_Run2012B-13Jul2012-v1_TrigStudy-mu_noTail_Tree.root", "out_SingleMu_Run2012B-13Jul2012-v1_TrigStudy-mu_noTail_Tree.root");

  
}


//================================================================================= plotHistograms 
void plotHistograms(eventInfo& evt, map<TString, plotMaker> &plot){
  plot["event"].addh1d("muMultiplicity", "muMultiplicity", "muon multiplicity", "events", 5, -0.5, 4.5, evt.muonP4->size(), evt.Weight);
  plot["event"].addh1d("elMultiplicity", "elMultiplicity", "electron multiplicity", "events", 5, -0.5, 4.5, evt.electronP4->size(), evt.Weight);
  plot["event"].addh1d("jetMultiplicity", "jetMultiplicity", "jet multiplicity", "events", 10, -0.5, 9.5, evt.jetP4->size(), evt.Weight);

  for(Int_t i=0,N=triggerFilterToCheckMu.size(); i<N; ++i){
    investigateTriggerFilterMu(evt, plot, triggerFilterToCheckMu.at(i));
  }
//   for(Int_t i=0,N=triggerFilterToCheckEl.size(); i<N; ++i){
//     investigateTriggerFilterEl(evt, plot, triggerFilterToCheckEl.at(i));
//   }
}


//================================================================================= investigateTriggerFilterMu 
//check, if the positive lepton can be matched to the trigger filter
//if yes, take this lepton as Tag
//then probe the other lepton
void investigateTriggerFilterMu(eventInfo& evt, map<TString, plotMaker> &plot, TString filterName, TString label){
  //--------------------------------------------------- do some consistency check - just to be sure
  if(evt.muonP4->size() != 2){
    cout<<"WARNING: The size of the muon vector is NOT two!"<<endl;
    cout<<"evt.muonP4->size()          = " << evt.muonP4->size() << endl;
    return;
  }  
  if(evt.muonP4->size() != evt.triggerFilterMu.size()){
    cout<<"WARNING: The size of the muon vector and Matching-info vector is different!"<<endl;
    cout<<"evt.muonP4->size()         = " << evt.muonP4->size() << endl;
    cout<<"evt.triggerFilterMu.size() = " << evt.triggerFilterMu.size() << endl;
    return;
  }
  if(evt.muonP4->size() != evt.muonCharge->size()){
    cout<<"WARNING: The size of the muon vector and muon-charge vector is different!"<<endl;
    cout<<"evt.muonP4->size()     = " << evt.muonP4->size() << endl;
    cout<<"evt.muonCharge->size() = " << evt.muonCharge->size() << endl;
    return;
  }
//   if(evt.electronP4->size() != evt.triggerFilterEl.size()){
//     cout<<"WARNING: The size of the electron vector and Matching-info vector is different!"
//     cout<<"evt.electronP4->size()      = " << evt.electronP4->size() << endl;
//     cout<<"evt.triggerFilterEl->size() = " << evt.triggerFilterEl->size() << endl;
//     return;
//   }
//   if(evt.electronP4->size() != evt.electronP4.size()){
//     cout<<"WARNING: The size of the electron vector and muon-charge vector is different!"
//     cout<<"evt.electronP4->size() = " << evt.electronP4->size() << endl;
//     cout<<"evt.electronP4->size() = " << evt.electronP4->size() << endl;
//     return;
//   }
  //--------------------------------------------------- check which lepton has positive and which negative charge
  Int_t index_pos = 0;
  Int_t index_neg = 1;  
  if(evt.muonCharge->at(0)==-1){
    index_pos = 0;
    index_pos = 1;
  }
  //--------------------------------------------------- check which lepton has been matched to the trigger filter 
  vector<bool> matched;
  for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
    matched.push_back(false);
    for(Int_t j=0,M=evt.triggerFilterMu.at(i).size(); j<M; ++j){
      if(filterName == evt.triggerFilterMu.at(i).at(j)){
        matched.back() = true;
        break;
      }
    }
  }
  //--------------------------------------------------- do the plotting according to the outcome, which lepton could be matched
  
  //if the positive lepton IS NOT matched => no Tag found => return
  if(!matched.at(index_pos)) return;
  
  //if the positive lepton IS matched => plot Tag and Probe, and passed and failed
  //plot the Tag
  plot[filterName].addh1d(filterName+"_pt_tag", filterName+"_pt_tag", "p_{T}(tag) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(index_pos).Pt(), evt.Weight);
  plot[filterName].addh1d(filterName+"_eta_tag", filterName+"_eta_tag", "#eta(tag)", "events", binEta, lowEta, upEta, evt.muonP4->at(index_pos).eta(), evt.Weight);
  //plot the Probe
  plot[filterName].addh1d(filterName+"_pt_probe", filterName+"_pt_probe", "p_{T}(probe) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(index_neg).pt(), evt.Weight);    
  plot[filterName].addh1d(filterName+"_eta_probe", filterName+"_eta_probe", "#eta(probe)", "events", binEta, lowEta, upEta, evt.muonP4->at(index_neg).eta(), evt.Weight);
  //plot the passed
  if(matched.at(index_neg)){
    plot[filterName].addh1d(filterName+"_pt_passed", filterName+"_pt_passed", "p_{T}(passed) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(index_neg).pt(), evt.Weight);    
    plot[filterName].addh1d(filterName+"_eta_passed", filterName+"_eta_passed", "#eta(passed)", "events", binEta, lowEta, upEta, evt.muonP4->at(index_neg).eta(), evt.Weight);
    plot[filterName].addh1d(filterName+"_PUVertices_passed", filterName+"_PUVertices_passed", "number of primary vertices (passed)", "events", binPU, lowPU, upPU, evt.PUInter, evt.Weight);
  }
  else{//plot the failed
    plot[filterName].addh1d(filterName+"_pt_failed", filterName+"_pt_failed", "p_{T}(failed) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(index_neg).pt(), evt.Weight);    
    plot[filterName].addh1d(filterName+"_eta_failed", filterName+"_eta_failed", "#eta(failed)", "events", binEta, lowEta, upEta, evt.muonP4->at(index_neg).eta(), evt.Weight);
    plot[filterName].addh1d(filterName+"_PUVertices_failed", filterName+"_PUVertices_failed", "number of primary vertices (failed)", "events", binPU, lowPU, upPU, evt.PUInter, evt.Weight);
  }
}


//================================================================================= makeEfficiencyPlots
//eff = passed/(passed+failed)
void makeEfficiencyPlots(map<TString, plotMaker> &plot){
  cout<<"----------> Create efficiency histograms"<<endl;

  vector<string> allTriggerFiltersToCheck = triggerFilterToCheckMu;
  for(Int_t i=0,N=triggerFilterToCheckEl.size(); i<N; ++i){
    allTriggerFiltersToCheck.push_back(triggerFilterToCheckEl.at(i));
  }
  vector<string> histNames;
  histNames.push_back("_pt_passed");
  histNames.push_back("_eta_passed");
  histNames.push_back("_PUVertices_passed");
  
  for(Int_t i=0,N=allTriggerFiltersToCheck.size(); i<N; ++i){
    TString triggerFilter = allTriggerFiltersToCheck.at(i);
    for(Int_t j=0,M=histNames.size(); j<M; ++j){
      TString passedName = triggerFilter;
      TString failedName = triggerFilter;
      TString effName    = triggerFilter;
      passedName += histNames.at(j);
      failedName += histNames.at(j);
      effName += histNames.at(j);
      failedName.ReplaceAll("passed", "failed");
      effName.ReplaceAll("passed", "eff");

      if(plot[triggerFilter].isNewh1d(passedName)) continue;
      if(plot[triggerFilter].isNewh1d(failedName)) continue;

      TH1D *h_passed = (TH1D*)(plot[triggerFilter].h1d[passedName])->Clone();
      TH1D *h_failed = (TH1D*)(plot[triggerFilter].h1d[failedName])->Clone();
      TH1D *den = (TH1D*)(h_passed)->Clone();
      den->Add((TH1D*)(h_failed)->Clone());

      TH1D *eff = (TH1D*)(h_passed)->Clone();
      eff->Divide(den);
      setEfficiencyErrors(eff, h_passed, h_failed);
      
      eff->SetName(effName);
      eff->SetTitle(effName);
      
      TString xTitle_eff = eff->GetXaxis()->GetTitle(); 
      xTitle_eff.ReplaceAll("(passed)",""); 
      eff->GetXaxis()->SetTitle(xTitle_eff);
      eff->GetYaxis()->SetTitle("efficiency");
      //eff->SetAxisRange(0.,1.,"Y");
      plot[triggerFilter].addh1d(eff);
    }
  }
}


//================================================================================= setEfficiencyErrors
//p = passed; f = failed
// f(eff) = p/(p+f)
// df/dt = 1/(p+f) - p/(p+f)^2
// df/dp = - p/(p+f)^2
// eff_sig = sqrt(p_sig^2*(df/dp)^2 + f_sig^2*(df/df)^2)
void setEfficiencyErrors(TH1D* &eff, TH1D* &h_passed, TH1D* &h_failed){
  if(h_passed->GetNbinsX() !=  h_failed->GetNbinsX()){
    cout<<"WARNING: the number of bins of the passed and failed histograms differ! The efficiency calculation might be depricated!"<<endl;
    cout<<"Number of bins (passed | failed): " << h_passed->GetNbinsX() << " | " << h_failed->GetNbinsX() << endl;
  }
  if(h_passed->GetXaxis()->GetXmax() !=  h_failed->GetXaxis()->GetXmax()){
    cout<<"WARNING: the upper x-range of the passed and failed histograms differ! The efficiency calculation might be depricated!"<<endl;
    cout<<"Number of bins (passed | failed): " << h_passed->GetXaxis()->GetXmax() << " | " << h_failed->GetXaxis()->GetXmax() << endl;
  }
  if(h_passed->GetXaxis()->GetXmin() !=  h_failed->GetXaxis()->GetXmin()){
    cout<<"WARNING: the lower x-range of the passed and failed histograms differ! The efficiency calculation might be depricated!"<<endl;
    cout<<"Number of bins (passed | failed): " << h_passed->GetXaxis()->GetXmin() << " | " << h_failed->GetXaxis()->GetXmin() << endl;
  }
  //loop over single bins of histograms an set the efficiency error
  for(Int_t i=1,N=h_passed->GetNbinsX(); i<=N; ++i){
    Double_t p = h_passed  ->GetBinContent(i);
    Double_t f = h_failed->GetBinContent(i);
    Double_t p_sig = h_passed  ->GetBinError(i);
    Double_t f_sig = h_failed->GetBinError(i);
    
    if(eff->GetBinContent(i)==0.) continue;
    
    Double_t dedt = 1/(p+f)-p/((p+f)*(p+f));
    Double_t dedp = -p/((p+f)*(p+f));
    
    Double_t eff_sig = sqrt(p_sig*p_sig*dedt*dedt + f_sig*f_sig*dedp*dedp);
    eff->SetBinError(i, eff_sig);
  }
}


//================================================================================= Loop
void Loop(TString inputFileName, TString outputFileName){
  //timer.Start("Loop");
  TString treeName = "trigStudyTree";
  map<TString, plotMaker> myPlots;      
  
  //--------------------------------------------------- open file and read in leafs
  TFile* inputFile = new TFile(inputFileName,"READ");
  if (!inputFile->IsOpen()) std::cout<<"Could not open file "<<inputFileName<<std::endl;
  TTree* inputTree= (TTree*)inputFile->Get(treeName);
  
  eventInfo event;
  initBranches(inputTree, event);
  //--------------------------------------------------- Event Loop
  int nEvents=inputTree->GetEntries();
  cout<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"////////////////////////////// Starting Event Loop ////////////////////////////"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Investigate File: " << inputFileName << endl;
  cout<<"Read out TTree:   " << treeName << endl;
  cout<<"Number of events: " << nEvents << endl;
  cout<<"-------------------------------------------------------------------------------"<<endl;
  cout<<endl;
  cout<<"Looping over events..."<<endl;
  nEvents=50000;
  for (int iEvent=0; iEvent<nEvents; ++iEvent){
    tools::progress(iEvent);
    inputTree->GetEntry(iEvent);
    handleTriggerMatchingInfo(event);
//     printEventInfo(event);
//      printMuonInfo(event);
     //printElectronInfo(event);
    plotHistograms(event, myPlots);
    //tools::progressReset();
  }
  
  makeEfficiencyPlots(myPlots);
  //--------------------------------------------------- Save all histograms in the appropriate folders
  TString option="RECREATE";
  for(typename map<TString, plotMaker>::iterator it = myPlots.begin(); it != myPlots.end(); ++it){
    it->second.savePlots(outputFileName, option, it->first);
    option="UPDATE";
  }
  //timer.Stop("Loop");
}


//================================================================================= handleTriggerMatchingInfo
// converting the trigger Info to a more practical format:  vector<string> --> vector<vector<string> >
void handleTriggerMatchingInfo(eventInfo& evt){
  evt.triggerFilterEl.clear();
  evt.triggerFilterMu.clear();
  //make simple check
  
  if(evt.electronP4->size() != evt.DESYtriggerElMatchedTriggerFilter->size()){
    cout<<"WARNING: Size of electron vector does not correspond to size of matching info vector!"<<endl;
    cout<<"evt.electronP4->size()                        = " << evt.electronP4->size() << endl;
    cout<<"evt.DESYtriggerElMatchedTriggerFilter->size() = " << evt.DESYtriggerElMatchedTriggerFilter->size() << endl;
    
    tools::printVector("evt.electronP4", *(evt.electronP4));
    tools::printVector("evt.DESYtriggerElMatchedTriggerFilter", *(evt.DESYtriggerElMatchedTriggerFilter));

    return;
  }
  if(evt.muonP4->size() != evt.DESYtriggerMuMatchedTriggerFilter->size()){
    cout<<"WARNING: Size of muon vector does not correspond to size of matching info vector!"<<endl;
    cout<<"evt.muonP4->size()                            = " << evt.muonP4->size() << endl;
    cout<<"evt.DESYtriggerMuMatchedTriggerFilter->size() = " << evt.DESYtriggerMuMatchedTriggerFilter->size() << endl;
    return;
  }
  
  for(Int_t i=0,N=evt.electronP4->size(); i<N; ++i){
    evt.triggerFilterEl.push_back(tools::splitToWords(evt.DESYtriggerElMatchedTriggerFilter->at(i), ";"));
    //cout<<"evt.DESYtriggerElMatchedTriggerFilter->at("<<i<<") = " << evt.DESYtriggerElMatchedTriggerFilter->at(i) << endl;
    //tools::printVector<TString>("triggerFilterEl", evt.triggerFilterEl.at(i));
  }
//     printMuonInfo(evt);
//     printElectronInfo(evt);
    
    //tools::printVector("evt.triggerFilterMu", evt.triggerFilterMu);
//     tools::printVector("evt.DESYtriggerMuMatchedTriggerFilter", *(evt.DESYtriggerMuMatchedTriggerFilter));
//     tools::printVector("evt.muonP4", *(evt.muonP4));
    
  for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
    evt.triggerFilterMu.push_back(tools::splitToWords(evt.DESYtriggerMuMatchedTriggerFilter->at(i), ";"));
    //cout<<"evt.DESYtriggerMuMatchedTriggerFilter->at("<<i<<") = " << evt.DESYtriggerMuMatchedTriggerFilter->at(i) << endl;
    //tools::printVector<TString>("triggerFilterMu", evt.triggerFilterMu.at(i));
  }  
  //cout<<"size: muonP4 | filterinfo = "
  //<<evt.muonP4->size() << " | "
  //<<evt.DESYtriggerMuMatchedTriggerFilter->size()
  //<<endl;
}


//================================================================================= initBranches
void initBranches(TTree *tree, eventInfo& evt){
  //--------------------------------------------------- initialize event info
  evt.DESYtriggerElMatchedTriggerFilter = 0;
  evt.DESYtriggerMuMatchedTriggerFilter = 0;
  evt.DESYtriggerNameMap                = 0;
  evt.electronCharge                    = 0;
  evt.electronP4                        = 0;
  evt.Event                             = 0;
  evt.HT                                = 0.;
  evt.jetP4                             = 0;
  evt.MET                               = 0.;
  evt.muonCharge                        = 0;
  evt.muonP4                            = 0;
  evt.prescaled                         = 0;
  evt.Run                               = 0;
  evt.PUInter                           = 0.;
  evt.PUWeight                          = 0.;
  evt.triggered                         = 0;
  evt.Weight                            = 0.;
  //--------------------------------------------------- assign branches to event structure
  tree->SetBranchAddress("DESYtriggerElMatchedTriggerFilter", &(evt.DESYtriggerElMatchedTriggerFilter));
  tree->SetBranchAddress("DESYtriggerMuMatchedTriggerFilter",&(evt.DESYtriggerMuMatchedTriggerFilter));
  tree->SetBranchAddress("DESYtriggerNameMap",&(evt.DESYtriggerNameMap));
  tree->SetBranchAddress("electronCharge",&(evt.electronCharge));
  tree->SetBranchAddress("electronP4",&(evt.electronP4));
  tree->SetBranchAddress("Event",&(evt.Event));
  tree->SetBranchAddress("HT",&(evt.HT));
  tree->SetBranchAddress("jetP4",&(evt.jetP4));
  tree->SetBranchAddress("MET",&(evt.MET));
  tree->SetBranchAddress("muonCharge",&(evt.muonCharge));
  tree->SetBranchAddress("muonP4",&(evt.muonP4));
  tree->SetBranchAddress("prescaled",&(evt.prescaled));
  tree->SetBranchAddress("Run",&(evt.Run));
  tree->SetBranchAddress("PUInter",&(evt.PUInter));
  tree->SetBranchAddress("PUWeight",&(evt.PUWeight));
  tree->SetBranchAddress("triggered",&(evt.triggered));
  tree->SetBranchAddress("Weight",&(evt.Weight));
}


//================================================================================= printEventInfo
void printEventInfo(eventInfo& evt, ostream& os){
  os<<endl;
  os<<"///////////////////////////////////////////////////////"<<endl;
  os<<"///////////////////// Event Info //////////////////////"<<endl;
  os<<"///////////////////////////////////////////////////////"<<endl;
  os<<endl;
  os<<"-----------------------------------------"<<endl;
  os<<"Run      = "<<evt.Run<<endl;
  os<<"Event    = "<<evt.Event<<endl;
  os<<"Weight   = "<<evt.Weight<<endl;
  os<<"PUWeight = "<<evt.PUWeight<<endl;
  os<<"PUInter  = "<<evt.PUInter<<endl;
  os<<"-----------------------------------------"<<endl;
  os<<"HT               = "<<evt.HT<<endl;
  os<<"MET              = "<<evt.MET<<endl;
  os<<"-----------------------------------------"<<endl;
  os<<"No. of Muons     = "<<evt.muonP4->size()<<endl;
  os<<"No. of Electrons = "<<evt.electronP4->size()<<endl;
  os<<"No. of Jets      = "<<evt.jetP4->size()<<endl;
  os<<"-----------------------------------------"<<endl;
  os<<"muonCharge.size()     = "<<evt.muonCharge->size()<<endl;
  os<<"electronCharge.size() = "<<evt.electronCharge->size()<<endl;
  os<<"-----------------------------------------"<<endl;
  os<<"MuMatchedTriggerFilter.size() = "<<evt.DESYtriggerMuMatchedTriggerFilter->size()<<endl;
  os<<"ElMatchedTriggerFilter.size() = "<<evt.DESYtriggerElMatchedTriggerFilter->size()<<endl;
  os<<"-----------------------------------------"<<endl;
  os<<"triggered.size()      = "<<evt.triggered->size()<<endl;
  os<<"prescaled.size()      = "<<evt.prescaled->size()<<endl;
  os<<"DESYTriggerMap.size() = "<<evt.DESYtriggerNameMap->size()<<endl;
  os<<"-----------------------------------------"<<endl;
  os<<endl;
  
  //tools::printVector<string>("jetP4", *(evt.jetP4));  
  //tools::printVector<string>("DESYtriggerMuMatchedTriggerFilter", *(evt.DESYtriggerMuMatchedTriggerFilter));
  //tools::printVector<string>("DESYtriggerElMatchedTriggerFilter", *(evt.DESYtriggerElMatchedTriggerFilter));
  //tools::printVector<string>("DESYtriggerNameMap", *(evt.DESYtriggerNameMap));
  //tools::printVector<string>("prescaled", *(evt.prescaled));
  //tools::printVector<string>("triggered", *(evt.triggered));
}


//================================================================================= printMuonInfo
void printMuonInfo(eventInfo& evt, ostream& os){
  os<<endl;
  os<<"///////////////////////////////////////////////////////"<<endl;
  os<<"////////////////////// Muon Info //////////////////////"<<endl;
  os<<"///////////////////////////////////////////////////////"<<endl;
  os<<endl;
  tools::printVector<LorentzM>("muonP4", *(evt.muonP4), os);
  tools::printVector<int>("muonCharge", *(evt.muonCharge), os);
  //------------------------ print invariant mass from all possible muon pairs
  if(evt.muonP4->size()>=2){
    os<<"-------------> Invariant masses of all possible muon pairs:"<<endl;
    for(Int_t i=0,N=evt.muonP4->size(); i<N-1; ++i){
      for(Int_t j=i+1; j<N; ++j){
        Double_t minv = (evt.muonP4->at(i) + evt.muonP4->at(j)).mass();
        os<<"Minv("<<i<<" + " <<j<<") = "<<minv<<endl;
      }
    }
    os<<endl;
  }
  //------------------------ print matching info
  for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
    os<<"----------------------------------------------------------------------------> Matches for muon "<<i<<":"<<
    endl;
    tools::printVector<TString>("matched TriggerFilter", evt.triggerFilterMu.at(i), os);
  }
  os<<endl;
}


//================================================================================= printElectronInfo
void printElectronInfo(eventInfo& evt, ostream& os){
  os<<endl;
  os<<"///////////////////////////////////////////////////////"<<endl;
  os<<"///////////////////// Electron Info ///////////////////"<<endl;
  os<<"///////////////////////////////////////////////////////"<<endl;
  os<<endl;
  tools::printVector<LorentzM>("electronP4", *(evt.electronP4), os);
  tools::printVector<int>("electronCharge", *(evt.electronCharge), os);
  //------------------------ print invariant mass from all possible electron pairs
  if(evt.electronP4->size()>=2){
    os<<"-------------> Invariant masses of all possible electron pairs:"<<endl;
    for(Int_t i=0,N=evt.electronP4->size(); i<N-1; ++i){
      for(Int_t j=i+1; j<N; ++j){
        Double_t minv = (evt.electronP4->at(i) + evt.electronP4->at(j)).mass();
        os<<"Minv("<<i<<" + " <<j<<") = "<<minv<<endl;
      }
    }
    os<<endl;
  }
  //------------------------ print matching info
  for(Int_t i=0,N=evt.electronP4->size(); i<N; ++i){
    os<<"----------------------------------------------------------------------------> Matches for electron "<<i<<":"<<
    endl;
    tools::printVector<TString>("matched TriggerFilter", evt.triggerFilterEl.at(i), os);
  }
  os<<endl;
}


//Thats another version of implementing the T&P method
//This version uses ALL Objects within an event to extract the trigger efficiency. 
//However, this method leads to a more complex treatment of the counting and the error propagation.
//The advantage is higher statistict. But the statistis is anayhow large enough...
// //================================================================================= investigateTriggerFilterMu 
// // eff = 2*NTT / (2*NTT + NTP)
// void investigateTriggerFilterMu(eventInfo& evt, map<TString, plotMaker> &plot, TString filterName){
//   //create vector which stores, if a muon is matched to the trigger filter with the name "filterName"
//   vector<bool> matched;
//   vector<int> probeIndex;
//   for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
//     matched.push_back(false);
//     for(Int_t j=0,M=evt.triggerFilterMu.at(i).size(); j<M; ++j){
//       if(filterName == evt.triggerFilterMu.at(i).at(j)){
//      matched.back() = true;
//      break;
//       }
//     }
//     if(!matched.at(i)) probeIndex.push_back(i);
//   }
//   
//   if(probeIndex.size()>2){
//     cout<<"WARNING: Something is awkward with the number of Probe leptons: " << probeIndex.size() << endl;
//     cout<<"The number of probe leptons is expected to be 1 or 2, because here a T&P method is applied!" << endl;
//   }
//   else if(probeIndex.size()==2) return; //make only plots if at least one Tag exists in an event
//   
//   //Fill the Tag histograms
//   //here, we have TWO Tag leptons -> use BOTH for the efficiency
//   if(probeIndex.size()==0){
//     for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
//       plot[filterName].addh1d("pt_tag", "pt_tag", "p_{T}(tag) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(i).Pt(), evt.Weight);
//       plot[filterName].addh1d("eta_tag", "eta_tag", "#eta(tag)", "events", binEta, lowEta, upEta, evt.muonP4->at(i).eta(), evt.Weight);
//     }
//     //the factor of '2' is due to the fact, that there are two matched tag leptons
//     plot[filterName].addh1d("PUVertices_tag", "PUVertices_tag", "number of primary vertices (tag)", "events", binPU, lowPU, upPU, evt.PUInter, 2*evt.Weight);
//   }
// 
//   //Fill the Probe histograms
//   //here, we have ONE Tag lepton -> the probe is NOT matched to the trigger filter  
//   if(probeIndex.size()==1){
//     plot[filterName].addh1d("pt_probe", "pt_probe", "p_{T}(probe) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(probeIndex.back()).pt(), evt.Weight);    
//     plot[filterName].addh1d("eta_probe", "eta_probe", "#eta(probe)", "events", binEta, lowEta, upEta, evt.muonP4->at(probeIndex.back()).eta(), evt.Weight);
//     plot[filterName].addh1d("PUVertices_probe", "PUVertices_probe", "number of primary vertices (probe)", "events", binPU, lowPU, upPU, evt.PUInter, evt.Weight);
//   }
// }
// 
// 
// //================================================================================= makeEfficiencyPlots
// //eff = tag/(tag+probe)
// void makeEfficiencyPlots(map<TString, plotMaker> &plot){
//   cout<<"----------> Create efficiency histograms"<<endl;
// 
//   vector<string> allTriggerFiltersToCheck = triggerFilterToCheckMu;
//   for(Int_t i=0,N=triggerFilterToCheckEl.size(); i<N; ++i){
//     allTriggerFiltersToCheck.push_back(triggerFilterToCheckEl.at(i));
//   }
//   vector<string> histNames;
//   histNames.push_back("pt_tag");
//   histNames.push_back("eta_tag");
//   histNames.push_back("PUVertices_tag");
//   
//   for(Int_t i=0,N=allTriggerFiltersToCheck.size(); i<N; ++i){
//     TString triggerFilter = allTriggerFiltersToCheck.at(i);
//     for(Int_t j=0,M=histNames.size(); j<M; ++j){
//       TString tagName   = histNames.at(j);
//       TString probeName = histNames.at(j);
//       TString effName   = histNames.at(j);
//       probeName.ReplaceAll("tag", "probe");
//       effName.ReplaceAll("tag", "eff");
//       if(plot[triggerFilter].isNewh1d(tagName)) continue;
//       if(plot[triggerFilter].isNewh1d(probeName)) continue;
// 
//       TH1D *h_tag   = (TH1D*)(plot[triggerFilter].h1d[tagName])->Clone();
//       TH1D *h_probe = (TH1D*)(plot[triggerFilter].h1d[probeName])->Clone();
//       TH1D *den = (TH1D*)(h_tag)->Clone();
//       den->Add((TH1D*)(h_probe)->Clone());
// 
//       TH1D *eff = (TH1D*)(h_tag)->Clone();
//       eff->Divide(den);
//       setEfficiencyErrors(eff, h_tag, h_probe);
//       
//       eff->SetName(effName);
//       eff->SetTitle(effName);
//       
//       TString xTitle_eff = eff->GetXaxis()->GetTitle(); 
//       xTitle_eff.ReplaceAll("(tag)",""); 
//       eff->GetXaxis()->SetTitle(xTitle_eff);
//       eff->GetYaxis()->SetTitle("efficiency");
//       //eff->SetAxisRange(0.,1.,"Y");
//       plot[triggerFilter].addh1d(eff);
//     }
//   }
// }
// 
// 
// //================================================================================= setEfficiencyErrors
// // f(eff) = t/(t+p)
// // df/dt = 1/(t+p) - t/(t+p)^2
// // df/dp = - t/(t+p)^2
// // eff_sig = sqrt(t_sig^2*(df/dt)^2 + p_sig^2*(df/dp)^2)
// void setEfficiencyErrors(TH1D* &eff, TH1D* &h_tag, TH1D* &h_probe){
//   if(h_tag->GetNbinsX() !=  h_probe->GetNbinsX()){
//     cout<<"WARNING: the number of bins of the tag and probe histograms differ! The efficiency calculation might be depricated!"<<endl;
//     cout<<"Number of bins (Tag | Probe): " << h_tag->GetNbinsX() << " | " << h_probe->GetNbinsX() << endl;
//   }
//   if(h_tag->GetXaxis()->GetXmax() !=  h_probe->GetXaxis()->GetXmax()){
//     cout<<"WARNING: the upper x-range of the tag and probe histograms differ! The efficiency calculation might be depricated!"<<endl;
//     cout<<"Number of bins (Tag | Probe): " << h_tag->GetXaxis()->GetXmax() << " | " << h_probe->GetXaxis()->GetXmax() << endl;
//   }
//   if(h_tag->GetXaxis()->GetXmin() !=  h_probe->GetXaxis()->GetXmin()){
//     cout<<"WARNING: the lower x-range of the tag and probe histograms differ! The efficiency calculation might be depricated!"<<endl;
//     cout<<"Number of bins (Tag | Probe): " << h_tag->GetXaxis()->GetXmin() << " | " << h_probe->GetXaxis()->GetXmin() << endl;
//   }
//   //loop over single bins of histograms an set the efficiency error
//   for(Int_t i=1,N=h_tag->GetNbinsX(); i<=N; ++i){
//     Double_t t = h_tag  ->GetBinContent(i);
//     Double_t p = h_probe->GetBinContent(i);
//     Double_t t_sig = h_tag  ->GetBinError(i);
//     Double_t p_sig = h_probe->GetBinError(i);
//     
//     if(eff->GetBinContent(i)==0.) continue;
//     
//     Double_t dedt = 1/(t+p)-t/((t+p)*(t+p));
//     Double_t dedp = -t/((t+p)*(t+p));
//     
//     Double_t eff_sig = sqrt(t_sig*t_sig*dedt*dedt + p_sig*p_sig*dedp*dedp);
//     eff->SetBinError(i, eff_sig);
//   }
// }


Revision:	1.5
Committed:	Thu Dec 6 12:27:27 2012 UTC (12 years, 5 months ago) by mstein
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+269 -110 lines
Log Message:	simple trigger efficiency estimation