source: git/modules/JetFakeParticle.cc@ 68dd0df

/*
 *  Delphes: a framework for fast simulation of a generic collider experiment
 *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */


/** \class JetFakeParticle
 *
 *  Converts jet into particle with some PID,
 *  according to parametrized probability.
 *
 *  \author M. Selvaggi - UCL, Louvain-la-Neuve
 *
 */

#include "modules/JetFakeParticle.h"

#include "classes/DelphesClasses.h"
#include "classes/DelphesFactory.h"
#include "classes/DelphesFormula.h"

#include "ExRootAnalysis/ExRootResult.h"
#include "ExRootAnalysis/ExRootFilter.h"
#include "ExRootAnalysis/ExRootClassifier.h"

#include "TMath.h"
#include "TString.h"
#include "TFormula.h"
#include "TRandom3.h"
#include "TObjArray.h"
#include "TDatabasePDG.h"
#include "TLorentzVector.h"

#include <algorithm>
#include <stdexcept>
#include <iostream>
#include <sstream>


using namespace std;

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

JetFakeParticle::JetFakeParticle() :
  fItInputArray(0)
{
}

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

JetFakeParticle::~JetFakeParticle()
{
}

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

void JetFakeParticle::Init()
{
  TFakeMap::iterator itEfficiencyMap;
  ExRootConfParam param;
  DelphesFormula *formula;
  Int_t i, size, pdg;

  // read efficiency formulas
  param = GetParam("EfficiencyFormula");
  size = param.GetSize();

  fEfficiencyMap.clear();
 
  for(i = 0; i < size/2; ++i)
  {
    formula = new DelphesFormula;
    formula->Compile(param[i*2 + 1].GetString());
    pdg = param[i*2].GetInt();
    
    if(TMath::Abs(pdg) != 11 && TMath::Abs(pdg) != 13 && TMath::Abs(pdg) != 22)
    {
      throw runtime_error("Jets can only fake into electrons, muons or photons. Other particles are not authorized.");
    }
    
    fEfficiencyMap[param[i*2].GetInt()] = formula;
  }

  // set default efficiency formula
  itEfficiencyMap = fEfficiencyMap.find(0);
  if(itEfficiencyMap == fEfficiencyMap.end())
  {
    formula = new DelphesFormula;
    formula->Compile("0.0");

    fEfficiencyMap[0] = formula;
  }
 
 
  // import input array

  fInputArray = ImportArray(GetString("InputArray", "FastJetFinder/jets"));
  fItInputArray = fInputArray->MakeIterator();

  // create output array

  fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "electrons_fake"));
  fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "muons_fake"));
  fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons_fake"));
  fJetOutputArray = ExportArray(GetString("JetOutputArray", "jets"));

}

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

void JetFakeParticle::Finish()
{
  if(fItInputArray) delete fItInputArray;

  TFakeMap::iterator itEfficiencyMap;
  DelphesFormula *formula;
  for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
  {
    formula = itEfficiencyMap->second;
    if(formula) delete formula;
  }
}

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

void JetFakeParticle::Process()
{
  Candidate *candidate, *fake;
  Double_t pt, eta, phi, e;
  TFakeMap::iterator itEfficiencyMap;
  DelphesFormula *formula;
  Int_t pdgCodeIn, pdgCodeOut, charge;
  Bool_t faked;
  
  Double_t p, r, rs, total;

  fItInputArray->Reset();
  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
  {
    const TLorentzVector &candidatePosition = candidate->Position;
    const TLorentzVector &candidateMomentum = candidate->Momentum;
    eta = candidatePosition.Eta();
    phi = candidatePosition.Phi();
    pt = candidateMomentum.Pt();
    e = candidateMomentum.E();
   
    r = gRandom->Uniform();
    total = 0.0;
    faked = 0;
    
    // loop over map for this jet
    for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
    {
    
  
    
      formula = itEfficiencyMap->second;
      pdgCodeOut = itEfficiencyMap->first;

      p = formula->Eval(pt, eta, phi, e);

      if(total <= r && r < total + p)
      {
        
        fake = static_cast<Candidate*>(candidate->Clone());
        faked = 1;
        
        // convert jet
        
        if(TMath::Abs(pdgCodeOut) == 11 || TMath::Abs(pdgCodeOut) == 13)
        {
          
        // for electrons and muons fake use the jet charge (if non-zero, otherwise randomly assign sign) 
          
          if(candidate->Charge != 0)fake->PID = -(candidate->Charge)*TMath::Abs(pdgCodeOut);
          else
          {
            rs = gRandom->Uniform();
            fake->PID = (rs < 0.5) ? -TMath::Abs(pdgCodeOut) : TMath::Abs(pdgCodeOut);
          } 
        
        }
        
        if(TMath::Abs(pdgCodeOut) == 22) fake->PID = 22;
        
        if(TMath::Abs(pdgCodeOut) == 11) fElectronOutputArray->Add(fake);
        if(TMath::Abs(pdgCodeOut) == 13) fMuonOutputArray->Add(fake);
        if(TMath::Abs(pdgCodeOut) == 22) fPhotonOutputArray->Add(fake);

        break;

      }

      total += p;
    }
 
  if(!faked) fJetOutputArray->Add(candidate);
 
 }
}

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