ReComputeMET.cpp

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include "caf_util/ReComputeMET.hpp"
#include "cafe/Config.hpp"

#include "TBranchElement.h"
#include "TDirectory.h"

using namespace cafe;
using namespace std;
using namespace metid;

namespace caf_util {

  ReComputeMET::ReComputeMET(const char *name) : cafe::Processor(name), 
                                                 _tree(0),_mcvars("_nchunks"),
                                                 _useDefaultElectron(false),
                                                 _useDefaultMuon(false)
  {
    cafe::Config config(name);

    _fromBranch = config.get("From", "");
    _toBranch   = config.get("To", "");
    _treeName   = config.get("Tree", name);
    _vars       = config.getVString("Variables", " ,");

    if(_fromBranch.empty()) {
      throw std::runtime_error(std::string("ReComputeMET[") + name + "] : From: branch in config file not set");
    }
    if(_toBranch.empty()) {
      throw std::runtime_error(std::string("ReComputeMET[") + name + "] : To: branch in config file not set");
    }
    out() << "ReComputeMET[" << name << "] From: " << _fromBranch << " To: " << _toBranch << std::endl ;
    out() << "ReComputeMET[" << name << "] Output Tree: " << _treeName << std::endl;

    _jetBranch  = config.get("jetBranch","JCCB"); 
    _muonBranch = config.get("muonBranch",""); 
    _electronBranch = config.get("electronBranch","");
    _useDefaultJet  = config.get("useDefaultGoodJetDef",true);

    if(_muonBranch.empty()) {
      _muonBranch = "Muon";
      _useDefaultMuon = true;
    }
    if(_electronBranch.empty()) {
      _electronBranch = "EMscone";
      _useDefaultElectron = true;
    }

    string jes_str  = config.get("JES","NotShifted");
    if (jes_str == "ShiftedPos") _jes = ShiftedPos ;
    else if (jes_str == "ShiftedNeg") _jes = ShiftedNeg ;
    else if (jes_str == "NotShifted") _jes = NotShifted ;
    else {
      err () << "ReComputeMET WARNING: JES config string \"" 
             << jes_str << "\" is not recongnized! Will use not shifted JES."
             << endl ;
      _jes = NotShifted ;
    }
      
    float jetConeSize = config.get("jetConeSize",0.5) ;
    if (jetConeSize != 0.5 && jetConeSize != 0.7) {
      err() << "RecomputeMET ERROR: jetConeSize should either 0.5 or 0.7, you value is \""
            << jetConeSize << "\". WILL USE 0.5 as default !" << endl ;
      jetConeSize = 0.5;
    }
    string jetAlgoType;
    if (jetConeSize==0.5) jetAlgoType = "JCCB";
    if (jetConeSize==0.7) jetAlgoType = "JCCA";

    out() << "ReComputeMET will use " << jetAlgoType << " jets. You're input jet branch should correspond to this jet algorithm" << endl;

    _valgo.push_back("corr"+jetAlgoType);
    _valgo.push_back("corrmu"+jetAlgoType);
    _valgo.push_back("smear_corr"+jetAlgoType);
    _valgo.push_back("smear_corrmu"+jetAlgoType);

    _debug = config.get("Debug", 0);

  }

  void ReComputeMET::inputFileOpened(TFile *file)
  {
    if(TTree *tmb_tree = (TTree *)file->Get("TMBTree")) {
      _tree = dynamic_cast<TTree*>(gROOT->FindObject(_treeName.c_str()));
      if (_tree == 0) {
        _tree = new TTree(_treeName.c_str(), fullName().c_str());
        tmb_tree->AddFriend(_tree);
      }

      _branch = _tree->GetBranch(_toBranch.c_str()) ;
      if (_branch == 0) {
        if (TBranchElement *br = (TBranchElement*) tmb_tree->GetBranch(_fromBranch.c_str())) {
          _result = new TMBMet(); 
          _branch = _tree->Branch(_toBranch.c_str(), &_result, 4096); 
        } else {
          err() << "ReComputeMET[" << name() << "] No such branch: " << _fromBranch << std::endl;
        } // TBranchElement
      } // _branch == 0
    } else {
      err() << "ReComputeMET[" << name() << "] No TMBTree" << std::endl;
    }

  }

  void ReComputeMET::inputFileClosing(TFile *file)
  {
    if(TTree* oldtree = dynamic_cast<TTree *>(gROOT->Get(_treeName.c_str()))) { oldtree->Delete(); }
    _tree = 0;
    _branch = 0;
  }

  bool ReComputeMET::processEvent(cafe::Event &event)
  {
    if (!_result) return false;
    const TMBMet* oldmet = event.get<TMBMet>(_fromBranch.c_str(), _vars);

    // accessing input branches
    Collection<TMBJet> jetcol(event.getCollection<TMBJet>(_jetBranch.c_str()));

    // keep cell based MET untouched
    _result->setMET(oldmet->getMET());
    _result->setMETnoeta(oldmet->getMETnoeta());
    _result->setMETC(oldmet->getMETC());
    _result->setMETD(oldmet->getMETD());
    _result->setMETA(oldmet->getMETA());
    _result->setMETB(oldmet->getMETB());
    _result->setVertex(oldmet->getZvertex());

    // recompute MET corrections
    const vector<BMetQualInfo>* voldmetqual = oldmet->getMetQualInfo(); 
    vector<BMetQualInfo>::const_iterator it;
    vector<BMetQualInfo> vnewmetqual;

    for (int ialgo=0;ialgo<4;ialgo++) {

      string requestedalgo = _valgo[ialgo];

      for (it=voldmetqual->begin();it!=voldmetqual->end();it++) {
        string thisalgo = it->getJetAlgo();
        if (strcmp(requestedalgo.c_str(),thisalgo.c_str()) != 0) continue;

        bool isMuCorrected = (requestedalgo.find("mu")<requestedalgo.size());
        bool isSmeared = (requestedalgo.find("smear")<requestedalgo.size());

        if (_debug) cout << "ReComputeMET, requestedalgo=" << requestedalgo << " isMuCorrected=" << isMuCorrected
                         << " isSmeared=" << isSmeared << endl;

        BMetStruct oldMETcorr = it->getMETcorr();       
        BMetStruct oldCHcorr  = it->getCHcorr();
        BMetStruct oldJEScorr = it->getJEScorr();
        BMetStruct oldEMcorr  = it->getEMcorr();
        BMetStruct oldMUcorr  = it->getMUcorr();
        BMetStruct oldMUCalcorr = it->getMUCalcorr(); 
        BMetStruct oldBJcorr  = it->getBJcorr();
        BMetStruct oldbasemet = oldMETcorr - oldCHcorr - oldJEScorr - oldEMcorr - oldMUcorr - oldMUCalcorr;

        BMetStruct newCHcorr;
        BMetStruct newJEScorr;
        BMetStruct newEMcorr;
        BMetStruct newMUcorr;
        BMetStruct newMUCalcorr;
        BMetStruct newBJcorr;

        //--------------------------------
        // Compute EM correction
        //--------------------------------
        if (_useDefaultElectron) {
          newEMcorr = oldEMcorr;
        } else {
          Collection<TMBEMCluster> emcol(event.getCollection<TMBEMCluster>(_electronBranch.c_str()));
          for (Collection<TMBEMCluster>::const_iterator it_em = emcol.begin(); it_em != emcol.end(); ++it_em) { 
             // not exactly what is in d0correct:
             float emuncor = it_em->uncorrE();
             float emx = it_em->Vect().Px() - emuncor*cos(it_em->Phi())*sin(it_em->Theta());
             float emy = it_em->Vect().Py() - emuncor*sin(it_em->Phi())*sin(it_em->Theta());
             float emz = it_em->Vect().Pz() - emuncor*cos(it_em->Theta());
             float ems = sqrt(emx*emx + emy*emy);
             if (_debug) cout << "ReComputeMET : EM : "
                              << setw(10) << emx << " , "
                              << setw(10) << emy << " , "
                              << setw(10) << emz << " , "
                              << setw(10) << ems << endl;
             newEMcorr.AddVis(emx,emy,emz,ems);
          } // for em
          newEMcorr.ComputeMet();
          if (_debug) {
            cout << " old EM correction: " << endl;
            oldEMcorr.Print(cout);
            cout << " new EM correction: " << endl;
            newEMcorr.Print(cout);
          }
        } // else _useDefaultElectron

        //--------------------------------
        // Compute Muon correction
        //--------------------------------
        // using all muons you put into the input branch
        if (_useDefaultMuon) {
          newMUcorr = oldMUcorr;
          newMUCalcorr = oldMUCalcorr;
        } else {
          Collection<TMBMuon> mucol(event.getCollection<TMBMuon>(_muonBranch.c_str()));
          int imu = 0;
          for(Collection<TMBMuon>::const_iterator it_mu = mucol.begin(); it_mu != mucol.end(); ++it_mu) {

            float mux = it_mu->Px();
            float muy = it_mu->Py();
            float muz = it_mu->Pz();
            float mus = it_mu->Pt(); 

            newMUcorr.AddVis(mux,muy,muz,mus);
            float mutheta = it_mu->Theta();
            float muphi   = it_mu->Phi();
            float eloss   = it_mu->eloss();
            mux = eloss * sin(mutheta) * cos(muphi);
            muy = eloss * sin(mutheta) * sin(muphi);
            muz = eloss * cos(mutheta);
            mus = -1. * sqrt(mux*mux + muy*muy);

            if (_debug) cout << "ReComputeMET : MUC " << imu << " : "
                             << setw(10) << mux << " , "
                             << setw(10) << muy << " , "
                             << setw(10) << muz << " , "
                             << setw(10) << mus << endl;

            newMUCalcorr.Add(mux,muy,muz,mus);
            ++imu;
          } // for muon
          newMUcorr.ComputeMet();
          newMUCalcorr.ComputeMet();
          if (_debug) {
            cout << " old MU correction: " << endl;
            oldMUcorr.Print(cout);
            cout << " new MU correction: " << endl;
            newMUcorr.Print(cout);
          }
        } // else _useDefaultMuon

        //--------------------------------------------------------------------
        // Compute good/bad jet corrections : CH fraction and JES corrections
        //--------------------------------------------------------------------
        // using only good jets you put into the input branch
        int ijet = 0;
        for(Collection<TMBJet>::const_iterator it_jet = jetcol.begin(); it_jet != jetcol.end(); ++it_jet) {

          TMBJet jet = *it_jet;
          if (isMuCorrected) jet.ActAsJESMUCorrected(); else jet.ActAsJESCorrected();
          if (isSmeared) {
            if (isMuCorrected) jet.ActAsSmearedMU(); else jet.ActAsSmeared();
          } // isSmeared

          if (jet.isGood() && !jet.isEM()) {    

            if (!_useDefaultJet) {
            // not exactly what is in d0correct:
            float chf = jet.chETfraction();
            float chx = chf*cos(jet.Phi());
            float chy = chf*sin(jet.Phi());
            float chz = chf*cos(jet.Theta());
            float chs = fabs(jet.sET()) * jet.chETfraction();
            if (_debug) cout << "ReComputeMET : CH  " << ijet << " : "
              << setw(10) << chx << " , "
              << setw(10) << chy << " , "
              << setw(10) << chz << " , "
              << setw(10) << chs << endl;
            newCHcorr.AddVis(chx,chy,chz,chs);
            }

            // JES correction
            // Warning : the z component cannot be obtained from jetcorr function "forMET"
            // Take the correction from jet instead of MET
            float jesx, jesy, jesz, jess, metjes, jetjes, smearfactor;
            if (isMuCorrected) {
              metjes = jet.JESMU_metC();
              if(_jes == ShiftedPos) 
                metjes +=  sqrt(pow(jet.JESMU_metdC_stat(),2)+pow(jet.JESMU_metdC_sys(),2)) ;
              else if(_jes == ShiftedNeg) 
                metjes -=  sqrt(pow(jet.JESMU_metdC_stat(),2)+pow(jet.JESMU_metdC_sys(),2)) ;

              jetjes = jet.JESMU_C();
              if(_jes == ShiftedPos) 
                metjes +=  sqrt(pow(jet.JESMU_dC_stat_up(),2)+pow(jet.JESMU_dC_sys_up(),2)) ;
              else if(_jes == ShiftedNeg) 
                metjes -=  sqrt(pow(jet.JESMU_dC_stat_down(),2)+pow(jet.JESMU_dC_sys_down(),2)) ;

              smearfactor = jet.SmearingFactorMU();
            } else {
              metjes = jet.JES_metC();
              if(_jes == ShiftedPos) 
                metjes +=  sqrt(pow(jet.JES_metdC_stat(),2)+pow(jet.JES_metdC_sys(),2)) ;
              else if(_jes == ShiftedNeg) 
                metjes -=  sqrt(pow(jet.JES_metdC_stat(),2)+pow(jet.JES_metdC_sys(),2)) ;

              jetjes = jet.JES_C();
              if(_jes == ShiftedPos) 
                metjes +=  sqrt(pow(jet.JES_dC_stat_up(),2)+pow(jet.JES_dC_sys_up(),2)) ;
              else if(_jes == ShiftedNeg) 
                metjes -=  sqrt(pow(jet.JES_dC_stat_down(),2)+pow(jet.JES_dC_sys_down(),2)) ;

              smearfactor = jet.SmearingFactor();
            }
            jesx = metjes * cos(jet.Phi());
            jesy = metjes * sin(jet.Phi());
            jesz = jetjes != 0 ? jet.Pz() * (1./jetjes - 1.) : 0;
            jess = jet.sET() * (jetjes - 1.);
            if (isSmeared) {
              jesx *= smearfactor;
              jesy *= smearfactor;
              jesz  = jetjes != 0 && smearfactor != 0 ? jet.Pz()  * (1./jetjes/smearfactor - 1.) : 0;
              jess  = jet.sET() * (jetjes*smearfactor - 1.);
            } // if isSmeared

            if (_debug) cout << "ReComputeMET : JES " << ijet << " : "
                             << setw(10) << jesx << " , "
                             << setw(10) << jesy << " , "
                             << setw(10) << jesz << " , "
                             << setw(10) << jess << endl;
            newJEScorr.Add(jesx,jesy,jesz,jess);

          } // isGood

          // @@ pxCH, pyCH are not stored in CAF for the moment !!
          /*
            if (jet.isBad()) {
            float bjx = jet.pxCH();
            float bjy = jet.pyCH();
            float bjz = jet.pzCH();
            float bjs = fabs(jet.sET()) * jet.chETfraction();
            if (_debug) cout << "ReComputeMET : BJ  " << ijet << " : "
            << setw(10) << bjx << " , "
            << setw(10) << bjy << " , "
            << setw(10) << bjz << " , "
            << setw(10) << bjs << endl;
            newBJcorr.Add(bjx,bjy,bjz,bjs);

            } // isBad
          */

          ++ijet;
        } // for jets
        if (_useDefaultJet) { newCHcorr = oldCHcorr; newBJcorr = oldBJcorr; }
        if (_debug) {
           cout << " old CH correction: " << endl;
           oldCHcorr.Print(cout);
           cout << " new CH correction: " << endl;
           newCHcorr.Print(cout);
        }

        newCHcorr.ComputeMet();
        newJEScorr.ComputeMet();
        newBJcorr.ComputeMet();

        //--------------------------------
        // Compute final MissingET
        //--------------------------------
        BMetStruct newMETcorr = oldbasemet + newCHcorr + newJEScorr + newEMcorr + newMUcorr + newMUCalcorr;

        BMetQualInfo result;
        result.setJetAlgo(thisalgo);

        result.setMETcorr(newMETcorr);
        result.setCHcorr(newCHcorr);
        result.setJEScorr(newJEScorr);
        result.setEMcorr(newEMcorr);
        result.setMUcorr(newMUcorr);
        result.setMUCalcorr(newMUCalcorr);
        result.setBJcorr(newBJcorr);

        // no recomputing of isolation variables for the moment
        result.setIsoPhiGoodJet(it->getIsoPhiGoodJet());
        result.setIsoPhiMu(it->getIsoPhiMu());
        result.setIsoPhiEM(it->getIsoPhiEM());
        result.setIsoPhiBadJet(it->getIsoPhiBadJet());
        result.setIsoPhiUE(it->getIsoPhiUE());
        result.setIsoPhiBadTower(it->getIsoPhiBadTower());

        vnewmetqual.push_back(result);
        break;
      } // for voldmetqual

    } // ialgo

    _result->setMetQualInfo(vnewmetqual);

    return true  ;
  }

} // namespace

ClassImp(caf_util::ReComputeMET);

---