ResSelector.cpp

Go to the documentation of this file.

#include "caf_util/ResSelector.hpp"
#include "kinem_util/AnglesUtil.hpp"
#include "cafe/Function.hpp"
#include "cafe/Config.hpp"
#include "cafe/Stat.hpp"
#include "TH1.h"
#include <iostream>
#include <string>
#include <sstream>

using namespace cafe; 
using namespace std;

namespace caf_util {
  ResSelector::ResSelector(const char *name) : Processor(name), 
                                               _vars1("fX", "fY", "fZ"), 
                                               _vars2("fX", "fY", "fZ"),
                                               _masshist(0) {
    
    Config config(name);
    branch1 = config.get("Branch1","Muon");
    branch2 =  config.get("Branch2",branch1);
    
    _target_mass = config.get("target_mass",-1.0);
    _low_window = config.get("lower_mass_bound", -1.0);
    _high_window = config.get("upper_mass_bound",-1.0);
    _hist_lower = config.get("hist_lower",0.);
    _hist_upper = config.get("hist_upper",120.);
    _hist_bins = config.get("hist_bins", 0);
    string type1 = config.get("Object1_type","TMBLorentzVector");
    string type2 = config.get("Object2_type","TMBLorentzVector");
    _pt1_method = config.get("Object1_method","best");
    _pt2_method = config.get("Object2_method","best");
    _useEOPetrack = config.get("UseEOPTrackMatchElectron", false) ;

    string charges = config.get("Charge","noselection");  // noselection or opposite or samesign
    
    if (charges.find("opposite") != string::npos) _charge_req = OPPOSITESIGN ;
    else if (charges.find("samesign") != string::npos) _charge_req = SAMESIGN ;
    else if (charges.find("noselection") != string::npos) _charge_req = NOSELECTION ;
    else {
      err() << "ResSelector WARNING: "
            << " option ChargeSign = \"" << charges << "\" is not reconized!"
            << endl ;
      exit(1) ;
    }
     
    //check the types and set up the appropriate case.  Muon-muon, mu-e, e-e, mu-jet, etc.
    
    if(type1.compare("TMBMuon")==0 || type1=="muon"|| type1=="Muon") {
      _obj1 = MUON ;
      //print an error message if someone tries to to the cal method with a muon
       if(_pt1_method=="cal" || _pt1_method == "calorimeter") 
         {
           err() << " ResSelector named \"" << name 
                 << "\" ERROR: "
                 <<"Calorimeter pT methods not available for muons."<<endl;
           exit(1) ;
         }       
    } else if(type1=="TMBEMCluster"||type1=="electron"|| type1=="Electron")
      _obj1 = ELECTRON ;
    else if(type1=="TMBJet"||type1=="jet"|| type1=="Jet") {
      _obj1 = JET;
       //error message if someone does the track method with a jet
      if(_pt1_method=="track" || _pt1_method == "Track") {
        err() << " ResSelector named \"" << name 
              << "\" ERROR: "
              << "Track pT determination methods not available for jets. " 
              << endl ;
         exit(1) ;
      }
    } else if (type1 != "TMBLorentzVector") {
      err() << " ResSelector named \"" << name 
            << "\" ERROR: Option \"Object1_type: " << type1 
            << "\" has wrong value \"" << type1 << "\""
            << endl ;
      exit(1) ;  
    } else _obj1 = LORENTZVECTOR;
    
    
    if(type2.compare("TMBMuon")==0 || type2=="muon"|| type2=="Muon") {
      _obj2 = MUON ;
      //print an error message if someone tries to to the cal method with a muon
      if(_pt1_method=="cal" || _pt1_method == "calorimeter") 
         {
           err() << " ResSelector named \"" << name 
                 << "\" ERROR: "
                 <<"Calorimeter pT methods not available for muons."<<endl;
           exit(1) ;
         }       
    } else if(type2=="TMBEMCluster"||type2=="electron"|| type2=="Electron")
      _obj2 = ELECTRON ;
     else if(type2=="TMBJet"||type2=="jet"|| type2=="Jet") {
       _obj2 = JET;
       //error message if someone does the track method with a jet
       if(_pt1_method=="track" || _pt1_method == "Track") {
         err() << " ResSelector named \"" << name 
             << "\" ERROR: "
               << "Track pT determination methods not available for jets. " 
               << endl ;
         exit(1) ;
       }
     } else if (type2 != "TMBLorentzVector") {
       err() << " ResSelector named \"" << name 
             << "\" ERROR: Option \"Object1_type: " << type2 
             << "\" has wrong value \"" << type2 << "\""
             << endl ;
       exit(1) ;         
     } else _obj2 = LORENTZVECTOR;


    if (_obj1==MUON) {
      _vars1.add("_chptr") ;
      _vars1.add("_hasCentral") ;
      if (charges != NOSELECTION) {
        _vars1.add("_best") ;
        _vars1.add("_localcorr._charge") ;
        _vars1.add("_centralcorr._charge") ;
        _vars1.add("_smearedMC._charge") ;
      }
    }

    if (_obj2==MUON) {
      _vars2.add("_chptr") ;
      _vars2.add("_hasCentral") ;
      if (charges != NOSELECTION) {
        _vars2.add("_best") ;
        _vars2.add("_localcorr._charge") ;
        _vars2.add("_centralcorr._charge") ;
        _vars2.add("_smearedMC._charge") ;
      }
    }

    if (_obj1==ELECTRON) {
      if (_useEOPetrack) {
        _vars1.add("_SpatialTrMatchChi2Prob") ;
        _vars1.add("_chptrBest") ;
      } else {
        _vars1.add("_SpatialTrMatchChi2ProbBest") ;
        _vars1.add("_chptr") ;
      }
      if (charges != NOSELECTION) _vars1.add("_q") ;
    }

    if (_obj2==ELECTRON) {
      if (_useEOPetrack) {
        _vars2.add("_SpatialTrMatchChi2Prob") ;
        _vars2.add("_chptrBest") ;
      } else {
        _vars2.add("_SpatialTrMatchChi2ProbBest") ;
        _vars2.add("_chptr") ;
      }
      if (charges != NOSELECTION) _vars2.add("_q") ;
    }
    
  }//end constructor
  
  bool ResSelector::processEvent(cafe::Event& event)    
  {
    StatPointer stat ;
    event.get("StatPointer", stat) ;
    
    
     /* switch for the different cases.   
        The calculation is done exactly the same way each time; 
        only the types of collections change, depending on the case.  
     */     
    

    const TMBLorentzVector *first = 0, *second = 0 ; // objects forming the mass
    
    const TClonesArray* col1 = event.getClonesArray(branch1,_vars1) ;
    const TClonesArray* col2 = event.getClonesArray(branch2,_vars2) ;
    
    if (!col1) {
      err() << "caf_util::ResSelector named \"" << name()
            << "\" ERROR: branch \"" << branch1 << "\" could not be found!"
            << endl ;
      return false ; 
    }

    if (!col2) {
      err() << "caf_util::ResSelector named \"" << name()
            << "\" ERROR: branch \"" << branch2 << "\" could not be found!"
            << endl ;
      return false ; 
    }

    double rec_mass = -1.0 ;
    
    for (int i=0; i <= col1->GetLast(); i++) {
      
      const TMBLorentzVector* v1 = (const TMBLorentzVector*) col1 -> At(i);
      
      double pt1  = -1.0, eta1 = -10.0, phi1 = -10.0 ;
      int charge1 = 0;
      
      if (_obj1 == MUON) {
        if (!GetMuon(* dynamic_cast<const TMBMuon*> (v1), 
                     _pt1_method, pt1, eta1, phi1, charge1)) continue ;
      } else if (_obj1 == ELECTRON) {
        if (!GetElectron(* dynamic_cast<const TMBEMCluster*> (v1),  
                         _pt1_method, pt1, eta1, phi1, charge1)) continue ;
      } else {
        pt1 = v1->Pt() ;
        eta1 = v1->Eta() ;
        phi1 = v1->Phi() ;
        charge1 = 0 ;
      }
      
     int j = 0 ;
     if (branch1 == branch2) {
       if (i < col2->GetLast()) j = i+1 ;
       else continue ;
     }
     for (; j <= col2->GetLast(); j++) {
      
        const TMBLorentzVector* v2 = (const TMBLorentzVector*) col2 -> At(j);
        
        double pt2  = -1.0, eta2 = -10.0, phi2 = -10.0 ;
        int charge2 = 0 ;
        
        if (_obj2 == MUON) {
          if (!GetMuon(* dynamic_cast<const TMBMuon*> (v2), 
                       _pt2_method, pt2, eta2, phi2, charge2)) continue ;
        } else if (_obj2 == ELECTRON) {
          if (!GetElectron(* dynamic_cast<const TMBEMCluster*> (v2),
                           _pt2_method, pt2, eta2, phi2, charge2)) continue ;
        } else {
          pt2 = v2->Pt() ;
          eta2 = v2->Eta() ;
          phi2 = v2->Phi() ;
          charge2 = 0 ;
        }
        
        if (!mass(pt1, eta1, phi1, charge1, pt2, eta2, phi2, charge2, rec_mass)) continue ;
        first = v1 ;
        second = v2 ;
      }
    }
    
    if(rec_mass < 0) return false;// false will end event processing   

    //stat pointer stuff
    stat.EventSelected(name()+": Invarian mass selection") ; 

    //save the objects that form the resonance
    const TMBEMCluster *pe  = 0 ;
    const TMBMuon *pm  = 0 ;
    const TMBJet *pj  = 0 ;
    if (first) {
      if (pe = dynamic_cast<const TMBEMCluster*>(first))
        event.put(name()+":EM1", new TMBEMCluster(*pe)) ; 

      else if (pm = dynamic_cast<const TMBMuon*>(first))
        event.put(name()+":Muon1", new TMBMuon(*pm)) ; 

      else if (pj = dynamic_cast<const TMBJet*>(first))
        event.put(name()+":Jet1", new TMBJet(*pj)) ; 

      else 
        event.put(name()+":LorentzVector1", first) ; 
    }
    
    pe = 0 ; pm = 0 ; pj = 0 ;

    if (second) {
      if (pe = dynamic_cast<const TMBEMCluster*>(second))
        event.put(name()+":EM2", new TMBEMCluster(*pe)) ; 

      else if (pm = dynamic_cast<const TMBMuon*>(second))
        event.put(name()+":Muon2", new TMBMuon(*pm)) ; 

      else if (pj = dynamic_cast<const TMBJet*>(second))
        event.put(name()+":Jet2", new TMBJet(*pj)) ; 

      else 
        event.put(name()+":LorentzVector2", second) ; 
    }
    
    //output the  mass as well
    event.put(name()+":ResMass", (double) rec_mass);
    if(_masshist) _masshist->Fill(rec_mass); 


    if (debug() > 1) 
      out() << "caf_util::ResSelector named \"" << name() 
            << "\" invariant mass " << rec_mass << " GeV" << endl ;

    return true; 
    
  }    //processEvent
  
  
  void ResSelector::begin() {
    getDirectory()->cd();
    if(_hist_bins>0)
      _masshist = new TH1D(("masshist_"+name()).c_str(), "Invariant mass",
                           _hist_bins, _hist_lower, _hist_upper); 
  } 
  
  
  bool ResSelector::mass(double pt1, double eta1, double phi1, int charge1,
                         double pt2, double eta2, double phi2, int charge2, 
                         double& mass) {
    if (pt1 < 0 || pt2 < 0 || 
        eta1 < -9.9 || eta2 < -9.9 ||
        phi1 < -9.9 || phi2 < -9.9 ) {
      err() << "caf_util::ResSelector named \"" << name()
            << "\" ERROR: pt or eta or phi out of range!" << endl 
            << " pt1 = " << pt1 << " eta1 = " << eta1 << " phi1 = " << phi1 << endl
            << " pt2 = " << pt2 << " eta2 = " << eta2 << " phi2 = " << phi2 << endl ;
      return -1.0 ;
    } 
    
    if (_charge_req != NOSELECTION && (charge1 ==0 || charge2 ==0) ) {
      warn() << "ResSelector WARNING: "
             << " The charge selection is set to " ;
      if (_charge_req == OPPOSITESIGN) warn() << " \"opposite\" " ;
      if (_charge_req == SAMESIGN) warn() << " \"samesign\" " ;
      warn() << " but charge is not defined for the selected objects!"
             << endl ;       
      return false ;
    }
    if(_charge_req == OPPOSITESIGN && 
       ((charge1 >= 0 && charge2 >= 0) || (charge1 <= 0 && charge2 <= 0)) ) return false ;
    if(_charge_req == SAMESIGN &&
       ((charge1 >= 0 && charge2 <= 0) || (charge1 <= 0 && charge2 >=0)) ) return false ;
       
    //do the invariant mass calculation, assuming mass of the doughter particles is negligible compare to the parent
    double calc_mass = 2.*pt1*pt2*(cosh(eta1 - eta2)- cos(kinem::signed_delta_phi(phi1,phi2))) ;
     if (calc_mass < 0) {
       err() << "ResSelector ERROR: "
             << " error in the invariant mass calculation."
             << endl ;
       return false ;
     }
     
     calc_mass = sqrt(calc_mass) ;
     
     //make sure the mass is within the user-defined window
     if(_low_window >= 0 && calc_mass < _low_window) return false ;
     if (_high_window >= 0 &&  calc_mass > _high_window) return false ;
     
     if (mass < 0 ) {
       mass = calc_mass ;
       return true ;
     }
     
     if (_target_mass >= 0) {
       // if this latest resonance was is closer than the previous minimum, 
       // choose this pair as the new candidate            
       if(fabs(calc_mass - _target_mass)< fabs(mass - _target_mass)) {
         mass = calc_mass ;
         return true ;
       }
     } else {
       if(calc_mass > mass) {
         mass = calc_mass ;
         return true ;
       }
     }
     return false ;   
  }

//---------------------------------------------------------------------------------

  bool ResSelector::GetMuon(const TMBMuon& muon, const string& pt_method,
                             double& pt, double& eta, double& phi, int& charge) {
     
     
     //determine where to get the pt from 
     if(pt_method == "track" || pt_method == "Track")
       {
         if (!muon.hasCentral()) {
           warn() << "ResSelector WARNING! "
                  << "You should use only muon with central track macthed when"
                  << " invarinat mass is calculated with track pT."
                  << endl ;
           return false ;
         }
         
         const TMBTrack* trk = muon.GetChargedTrack() ;
         if(!trk) {
           err() << "ResSelector ERROR: "
                 << " Track pointer is 0 for the track macthed muon!"
                 << endl ;
           return false ;
         }
         
         pt = trk->Pt();
         phi = trk->Phi();
         eta = trk->Eta();
         charge = trk->charge() ;
         
       } else {       
       eta = muon.Eta();
       pt = muon.Pt(); 
       phi = muon.Phi();
       charge = muon.charge() ;              
     }

     return true ;
   }


  bool ResSelector::GetElectron(const TMBEMCluster& electron, const string& pt_method,
                                double& pt, double& eta, double& phi, int& charge) {
    
    //determine where to get the pt from 
    if(pt_method == "track" || pt_method == "Track") {
        const TMBTrack* trk  = 0 ;
        if ( (_useEOPetrack && !electron.has_track_match(0.0) ) ||
             (!_useEOPetrack && !electron.has_spatial_track_match(0.0))) {
          warn() << "ResSelector WARNING! "
                 << "You should use only electron with track macthed when"
                 << " invarinat mass is calculated with track pT."
                 << endl ;
          return false ;
        }

        if (!_useEOPetrack) 
          trk = electron.GetChargedTrack() ;
        else 
          trk = electron.getPtrChp() ;

        if(!trk) {
          err() << "ResSelector ERROR: "
                << " Track pointer is 0 for the track macthed electron!"
                << endl ;
          return false ;
        }
        
        pt = trk->Pt();
        phi = trk->Phi();
        eta = trk->Eta();
        charge = trk->charge() ;
        
      } else {       
      
      eta = electron.Eta();
      pt = electron.Pt(); 
      phi = electron.Phi();
      if (!_useEOPetrack) charge = (int) electron.charge() ;
      else {
        const TMBTrack* trk  = 0 ;
        if ( !electron.has_track_match(0.0)) {
          warn() << "ResSelector WARNING! "
                  << "You should use only electron with track macthed when"
                 << " determine charge with EOP track matched electron."
                 << endl ;
         return false ;
        }
        
        trk = electron.getPtrChp() ;
        if(!trk) {
          err() << "ResSelector ERROR: "
                << " Track pointer is 0 for the track macthed electron!"
                << endl ;
          return false ;
        }
        charge = (int) trk->charge() ;
      }
    }
    return true ;
  }

}//namespace caf_util

  
ClassImp(caf_util::ResSelector);

---