TMBMuon.cpp

Go to the documentation of this file.

//                                                                      //
//      TMBTree Muon class                                              //
//                                                                      //

#include "tmb_tree/TMBMuon.hpp"
#include "TError.h"
#include "TMath.h"
#include "tmb_tree/TMBTrack.hpp"
#include "kinem_util/AnglesUtil.hpp"
#include <iostream>

using namespace std;

namespace {
  Double_t MuonMass = 0.1057;
}

ClassImp(TMBMuon)

//_____________________________________________________________________________
TMBMuon::TMBMuon()
    :  TMBLorentzVector(0,0,0,MuonMass, kPtEtaPhiM),
       _best(kCentralCorr)
{
}

TMBMuon::TMBMuon(Double_t pT, Double_t eta, Double_t phi, Type best,
                 const TMBMuonType& local,   const TMBMuonType& localcorr,
                 const TMBMuonType& central, const TMBMuonType& centralcorr,
                 const TMBMuonType& global,  const TMBMuonType& smearedMC)
    : TMBLorentzVector(pT,eta,phi,MuonMass,kPtEtaPhiM),
      _local(local),
      _localcorr(localcorr),
      _central(central),
      _centralcorr(centralcorr),
      _global(global),
      _smearedMC(smearedMC),
      _best(best)
{
}

void TMBMuon::SetLocalInfo(Int_t nhit, Int_t ndeck, Int_t region, Int_t octant,
                Float_t chisqloc, Float_t sctimeA, Float_t sctimeB, Float_t sctimeC,
                Float_t xA, Float_t yA, Float_t zA, Float_t bdl) 
{
  _nhit     = nhit;
  _ndeck    = ndeck;
  _region   = region;
  _octant   = octant;
  _chisqloc = chisqloc;
  _sctimeA  = sctimeA;
  _sctimeB  = sctimeB;
  _sctimeC  = sctimeC;
  _xA       = xA;
  _yA       = yA;
  _zA       = zA;
  _bdl      = bdl;
}

void TMBMuon::SetMatchingInfo(Int_t nseg, Int_t ndof,
                Float_t deltaPhi, Float_t deltaEta, Float_t deltaDrift,
                Float_t chisq, Float_t zAtPca, Float_t impPar, Float_t impParSig,
                Float_t err_zAtPca, Float_t err_impPar, Float_t dca,
                TRef chptr, TRef vtxref)
{
  _nseg       = nseg;
  _ndof       = ndof;
  _deltaPhi   = deltaPhi;
  _deltaEta   = deltaEta;
  _deltaDrift = deltaDrift;
  _chisq      = chisq;
  _zAtPca     = zAtPca;
  _impPar     = impPar;
  _impParSig  = impParSig;
  _err_zAtPca = err_zAtPca;
  _err_impPar = err_impPar;
  _dca        = dca;
  _chptr      = chptr;
  _vtxref     = vtxref;
}

void TMBMuon::SetDetEta(Float_t deteta)
{
  _deteta     = deteta;
}

void TMBMuon::SetMTCInfo(Int_t nmtc, Int_t calnLayer, Float_t etrack_best,
                Float_t caleSig, Float_t calEta, Float_t calPhi, Float_t eloss)
{
  _nmtc        = nmtc;
  _calnLayer   = calnLayer;
  _etrack_best = etrack_best;
  _caleSig     = caleSig;
  _calEta      = calEta;
  _calPhi      = calPhi;
  _eloss       = eloss;
}

void TMBMuon::SetIsolationInfo(Int_t nTrk5, Float_t EInCone1, Float_t EInCone15, Float_t EInCone2,
                Float_t EInCone4, Float_t EInCone6,
                Float_t drJet5, Float_t etTrkCone5, Float_t etHalo)
{
  _nTrk5      = nTrk5;
  _EInCone1   = EInCone1;
  _EInCone15  = EInCone15;
  _EInCone2   = EInCone2;
  _EInCone4   = EInCone4;
  _EInCone6   = EInCone6;
  _drJet5     = drJet5;
  _etTrkCone5 = etTrkCone5;
  _etHalo     = etHalo;
}
 
void TMBMuon::SetFlags(Int_t isLoose, Int_t isMedium, Int_t isTight,
                Int_t hasLocal, Int_t hasCentral, Int_t hasCal,
                Int_t isMuonEventOK, Int_t isCosmic, Int_t isCosmicT)
{
  _isLoose       = isLoose;
  _isMedium      = isMedium;
  _isTight       = isTight;
  _hasLocal      = hasLocal;
  _hasCentral    = hasCentral;
  _hasCal        = hasCal;
  _isMuonEventOK = isMuonEventOK;
  _isCosmic      = isCosmic;
  _isCosmicT     = isCosmicT;
}

void TMBMuon::SetExpectedHits(Int_t expWhitsA, Int_t expWhitsBC, 
                              Int_t expShitsA, Int_t expShitsBC)
{
  _expWhitsA      = expWhitsA;
  _expWhitsBC     = expWhitsBC;
  _expShitsA      = expShitsA;
  _expShitsBC     = expShitsBC;
}

void TMBMuon::SetMCsmearingRand(Float_t rand)
{
  _rand = rand;
}

void TMBMuon::CorrectPt(double dca)
{
  //If muo does not have track do not try to correct
  if(!_hasCentral) return;
  const TMBTrack* trk = GetChargedTrack();
  if(trk) {
    double err_rqpt = trk->trerrs(4, 0);
    double err_rr = trk->trerrs(0, 0);    
    float qopt = trk->qpt();
    qopt -= dca * err_rqpt / err_rr;
    double pTcorr = 1 / qopt;
    int q = 1;
    if(pTcorr<0) {
      pTcorr *= -1;
      q = -1;
    }
    double scale = pTcorr / trk->Pt();
    double corrpx = scale * trk->Px();
    double corrpy = scale * trk->Py();
    double corrpz = scale * trk->Pz();
    double corrE = scale * trk->E();
    
    //change TMBMuon 4-vector
    SetPxPyPzE(corrpx,corrpy,corrpz,corrE);
    //change _centralcorr 4-vector & charge
    TMBMuonType new_centralcorr(pTcorr,Eta(),Phi(),err_pT(),err_eta(),err_phi(),q);
    _centralcorr = new_centralcorr;
    //Set best muon info to CentralCorr
    _best = TMBMuon::kCentralCorr; 
  }

}

Int_t TMBMuon::charge() const
{
  switch(_best) {
    case kLocalCorr: return _localcorr.charge();
    case kCentralCorr: return _centralcorr.charge();
    case kSmearedMC: return _smearedMC.charge();
    default: Error("TMBMuon::charge","bad best muon type");
  }
  return 0;
}

Double_t TMBMuon::tlm() const
{
  switch(_best) {
    case kLocalCorr: return _localcorr.tlm();
    case kCentralCorr: return _centralcorr.tlm();
    case kSmearedMC: return _smearedMC.tlm();
    default: Error("TMBMuon::tlm","bad best muon type");
  }
  return 0;
}

Double_t TMBMuon::err_phi() const
{
  switch(_best) {
    case kLocalCorr: return _localcorr.err_phi();
    case kCentralCorr: return _centralcorr.err_phi();
    case kSmearedMC: return _smearedMC.err_phi();
    default: Error("TMBMuon::err_phi","bad best muon type");
  }
  return 0;
}

Double_t TMBMuon::err_eta() const
{
  switch(_best) {
    case kLocalCorr: return _localcorr.err_eta();
    case kCentralCorr: return _centralcorr.err_eta();
    case kSmearedMC: return _smearedMC.err_eta();
    default: Error("TMBMuon::err_eta","bad best muon type");
  }
  return 0;
}

Double_t TMBMuon::err_pT() const
{
  switch(_best) {
    case kLocalCorr: return _localcorr.err_pT();
    case kCentralCorr: return _centralcorr.err_pT();
    case kSmearedMC: return _smearedMC.err_pT();
    default: Error("TMBMuon::err_pT","bad best muon type");
  }
  return 0;
}

Int_t TMBMuon::wireHitsA() const
{
  return (_nhit % 10);
}

Int_t TMBMuon::wireHitsB() const
{
  return (_nhit / 10)%10;
}

Int_t TMBMuon::wireHitsC() const
{
  return (_nhit / 100)%10;
}

Int_t TMBMuon::wireHitsBC() const
{
  return (wireHitsB() + wireHitsC());
}

int TMBMuon::scintHitsA() const {
  return (_nhit / 1000)%10;
}

int TMBMuon::scintHitsBC() const {
  return scintHitsB() + scintHitsC();
}

int TMBMuon::scintHitsB() const {
  return (_nhit/10000)%10;
}

int TMBMuon::scintHitsC() const {
  return (_nhit/100000)%10;
}

Double_t TMBMuon::GetEtaDetector(Double_t eta, Double_t z, Double_t rdet) const
{
  float z0 = z + rdet*sinh(eta);
  float temp = z0/rdet + sqrt(z0*z0/rdet/rdet + 1.);
  return (temp>0?log(temp):eta);
}

bool TMBMuon::isInHole(Int_t solpol) const
{
 double deta = GetEtaDetector(solpol);
 return ( fabs(deta)<1.25 && this->Phi()>4.25 && this->Phi()<5.15);
}

Double_t TMBMuon::GetEtaDetector(Int_t solpol) const
{
 double theta = this->Theta();
 double phi   = this->Phi();
 double qpt   = (this->Pt()==0)?0:this->charge()/this->Pt();
 phi+= solpol *qpt *0.36; // correction for bending in soleneoid
 double z0    = _zAtPca;
 double deta = -99.;
 TMBTrack* track = this->GetChargedTrack();
 if(track){
   theta = track->Theta();
   phi = track->Phi();
   z0  = track->z();
   deta = det_eta(theta, phi, z0);
 } else {
   deta = kinem::eta(_xA,_yA,_zA);
 }

 return deta;
}

Double_t TMBMuon::det_eta(double theta, double phi, double z0) const
{
   float x_A = 292.;
   float y_A = 294.;
   float z_A = 442.;
   float pi = kinem::PI;

   if(phi<0) phi += 2*pi;
   if(phi>2*pi) phi -= 2*pi;

   double x = 0;
   double y = 0;
   double z = 0;
   int forward = 0;

   if(theta < pi/2) z = z_A;
   else z = -z_A;

   //see if we hit the forward or backward system:
   double r = fabs( (z - z0)*tan(theta) );
   x = r*cos(phi);
   y = r*sin(phi);
   if (fabs(x)<x_A && fabs(y)<y_A) forward = 1;

   //otherwise calculate x, y and z for central:
   if(!forward){
     if(phi > 7*pi/4 || phi < pi/4){
       x = x_A;
       y = x* tan(phi);
     }
     else if(phi > pi/4 && phi < 3*pi/4){
       y = y_A;
       x = y/tan(phi);
     }
     else if(phi > 3*pi/4 && phi < 5*pi/4){
       x = -x_A;
       y = x* tan(phi);
     }
     else if(phi > 5*pi/4 && phi < 7*pi/4){
       y = -y_A;
       x = y/tan(phi);
     }

     r = sqrt(x*x + y*y);
     z = z0 + r/tan(theta);
   }

   float eta = -log( tan( atan2(r,z)/2));
   return eta;
}

Double_t TMBMuon::chisqProb() const
{
  return TMath::Prob(_chisq,_ndof);
}


// Interface to apply Muon Smearing added 8 Feb 2006
void TMBMuon::ActAsSmeared(Float_t * par, bool zjpsiparams) {
  float rndnum=MCsmearingRand();
 if ( rndnum>98 ) return; // random number not set
//   if ( rndnum>98 ){
//       rndnum=RandGauss::shoot();   // random number used for smearing
//       SetMCsmearingRand(rndnum);
//   }
  if ( !hasCentral() || GetChargedTrack() == NULL  ) return;
  int numSMTHits = GetChargedTrack()->nsmt();      // number of hits in the SMT
  int numCFTHits = GetChargedTrack()->ncft();       // number of hits in the CFT
  int ncft_cut=0;
  float eta_cut=100;
  float eta_cft_cut=100;
  int ipar_off=0;       // parameter offset +1 for large eta or few cft hit
                        //                  +2 for no smt hits
  TMBMuonType Muon;
  float q_over_pt;
  if (numSMTHits == 0 ) {                      // smear vertex corrected pt if no SMT hits
    Muon=CentralCorr();
    ipar_off+=2;
    ncft_cut=int(par[I_nosmt_ncft_cut_TMBMuon]+0.5);
    eta_cut =par[I_nosmt_eta_cut_TMBMuon];
    eta_cft_cut =par[I_nosmt_eta_cft_cut_TMBMuon];
  } else {
    Muon=Central();
    ncft_cut=int(par[I_smt_ncft_cut_TMBMuon]+0.5);
    eta_cut =par[I_smt_eta_cut_TMBMuon];
    eta_cft_cut =par[I_smt_eta_cft_cut_TMBMuon];
  }
  float det_etaCFT=0;
  if ( GetChargedTrack() != NULL )  det_etaCFT=GetChargedTrack()->det_etaCFT();
  if ( Muon.Pt()< 1e-4) return;

  if (!zjpsiparams) {
    if ( fabs(det_etaCFT) >  eta_cft_cut || fabs(Muon.Eta() ) >  eta_cut || numCFTHits <  ncft_cut )
      ipar_off+=1;
    float par_A=par[ipar_off+I_par_A_TMBMuon];
    float par_B=par[ipar_off+I_par_B_TMBMuon];
    float par_Bpt=par[ipar_off+I_par_Bpt_TMBMuon];
    float scale=par[ipar_off+I_scale_TMBMuon];
    float width=par_A+par_B/Muon.Pt()+par_Bpt*Muon.Pt();
    q_over_pt = scale*Muon.charge()/Muon.Pt()+width*rndnum;
  } else {
    if ( fabs(det_etaCFT) > eta_cft_cut) ipar_off+=1;
    float par_A=par[ipar_off+I_par_A_TMBMuon];
    float par_B=par[ipar_off+I_par_B_TMBMuon];
    float coshpt = cosh(Muon.Eta())/Muon.Pt()/Muon.Pt();
    float width  =  sqrt(par_A*par_A + par_B*par_B*coshpt);
    q_over_pt    = Muon.charge()/Muon.Pt()+width*rndnum;  
    // cout << " numSMTHits=" << numSMTHits << ", det_etaCFT=" << det_etaCFT << ", par_A=" << par_A << ", par_B=" << par_B << endl;
  } // else name.find("zjpsi")

  float pT=Muon.Pt();
  if (q_over_pt != 0 ) pT=fabs(1/q_over_pt);
  int charge_sm=1;
  if (q_over_pt<0)  charge_sm=-1;
  TMBMuonType smear_muon(pT, Muon.Eta(),Muon.Phi(), 
        Muon.err_pT(), Muon.err_eta(), Muon.err_phi(),
        charge_sm);
  _smearedMC = smear_muon; 
  _best=  kSmearedMC;
  this->SetPtEtaPhiM(pT,this->Eta(),this->Phi(),this->M());
  return;

}

 
void TMBMuon::ActAsSmeared_pre()
{
  float par[npar_TMBMuon];
// smaring parameter for pre shutown data 
  par[I_par_A_TMBMuon+0]=0.0011; 
  par[I_par_A_TMBMuon+1]=0.0000; 
  par[I_par_A_TMBMuon+2]=0.0009; 
  par[I_par_A_TMBMuon+3]=0.0009; 
  par[I_par_B_TMBMuon+0]=0.013; 
  par[I_par_B_TMBMuon+1]=0.000; 
  par[I_par_B_TMBMuon+2]=0.018; 
  par[I_par_B_TMBMuon+3]=0.018; 
  par[I_par_Bpt_TMBMuon+0]=0.0;
  par[I_par_Bpt_TMBMuon+1]=0.0;
  par[I_par_Bpt_TMBMuon+2]=0.0;
  par[I_par_Bpt_TMBMuon+3]=0.0;
  par[I_scale_TMBMuon+0]=1; 
  par[I_scale_TMBMuon+1]=1; 
  par[I_scale_TMBMuon+2]=1; 
  par[I_scale_TMBMuon+3]=1; 
  par[I_smt_ncft_cut_TMBMuon]=0;
  par[I_nosmt_ncft_cut_TMBMuon]=0;
  par[I_smt_eta_cut_TMBMuon]=1.5;
  par[I_nosmt_eta_cut_TMBMuon]=1.5;
  par[I_smt_eta_cft_cut_TMBMuon]=0.;
  par[I_nosmt_eta_cft_cut_TMBMuon]=0.;
   ActAsSmeared(par,"dummy");
  return;
}
 
void TMBMuon::ActAsSmeared_post()
{
  float par[npar_TMBMuon];
// smaring parameter for pre shutown data 
  par[I_par_A_TMBMuon+0]=0.0015; 
  par[I_par_A_TMBMuon+1]=0.0016; 
  par[I_par_A_TMBMuon+2]=0.0017; 
  par[I_par_A_TMBMuon+3]=0.0017; 
  par[I_par_B_TMBMuon+0]=0.018; 
  par[I_par_B_TMBMuon+1]=0.019; 
  par[I_par_B_TMBMuon+2]=0.034; 
  par[I_par_B_TMBMuon+3]=0.034; 
  par[I_par_Bpt_TMBMuon+0]=0.0;
  par[I_par_Bpt_TMBMuon+1]=0.0;
  par[I_par_Bpt_TMBMuon+2]=0.0;
  par[I_par_Bpt_TMBMuon+3]=0.0;
  par[I_scale_TMBMuon+0]=1; 
  par[I_scale_TMBMuon+1]=1; 
  par[I_scale_TMBMuon+2]=1; 
  par[I_scale_TMBMuon+3]=1; 
  par[I_smt_ncft_cut_TMBMuon]=0;
  par[I_nosmt_ncft_cut_TMBMuon]=0;
  par[I_smt_eta_cut_TMBMuon]=1.5;
  par[I_nosmt_eta_cut_TMBMuon]=1.5;
  par[I_smt_eta_cft_cut_TMBMuon]=0.;
  par[I_nosmt_eta_cft_cut_TMBMuon]=0.;
  ActAsSmeared(par,"dummy");
  return;
}




void TMBMuon:: FixMuonPt(){
  TMBMuonType* best;
  switch(_best) {
  case kLocalCorr:    best= &_localcorr; break;
  case kCentralCorr:  best=&_centralcorr; break;
  case kSmearedMC  :  best=&_smearedMC; break;
  default:  Error("TMBMuon::FixMuonPt","bad best muon type");
 }


  SetPxPyPzE(best->X(),best->Y(),best->Z(),best->E());

}


//______________________________________________________________________________

---