source: HiSusy/trunk/Delphes-3.0.0/modules/Calorimeter.cc @ 1

/** \class Calorimeter
 *
 *  Fills calorimeter towers, performs calorimeter resolution smearing
 *  and preselects towers hit by electrons, photons and neutral particles.
 *
 *  $Date: 2012-11-19 14:05:19 +0100 (Mon, 19 Nov 2012) $
 *  $Revision: 829 $
 *
 *
 *  \author P. Demin - UCL, Louvain-la-Neuve
 *
 */

#include "modules/Calorimeter.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;

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

Calorimeter::Calorimeter() :
  fECalResolutionFormula(0), fHCalResolutionFormula(0), fItInputArray(0)
{
  fECalResolutionFormula = new DelphesFormula;
  fHCalResolutionFormula = new DelphesFormula;
}

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

Calorimeter::~Calorimeter()
{
  if(fECalResolutionFormula) delete fECalResolutionFormula;
  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
}

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

void Calorimeter::Init()
{
  ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
  Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
  Double_t ecalFraction, hcalFraction;
  TBinMap::iterator itEtaBin;
  set< Double_t >::iterator itPhiBin;
  vector< Double_t > *phiBins;

  // read eta and phi bins
  param = GetParam("EtaPhiBins");
  size = param.GetSize();
  fBinMap.clear();
  fEtaBins.clear();
  fPhiBins.clear();
  for(i = 0; i < size/2; ++i)
  {
    paramEtaBins = param[i*2];
    sizeEtaBins = paramEtaBins.GetSize();
    paramPhiBins = param[i*2 + 1];
    sizePhiBins = paramPhiBins.GetSize();

    for(j = 0; j < sizeEtaBins; ++j)
    {
      for(k = 0; k < sizePhiBins; ++k)
      {
        fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
      }
    }
  }

  // for better performance we transform map of sets to parallel vectors:
  // vector< double > and vector< vector< double >* >
  for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
  {
    fEtaBins.push_back(itEtaBin->first);
    phiBins = new vector< double >(itEtaBin->second.size());
    fPhiBins.push_back(phiBins);
    phiBins->clear();
    for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
    {
      phiBins->push_back(*itPhiBin);
    }
  }
  
  // read energy fractions for different particles
  param = GetParam("EnergyFraction");
  size = param.GetSize();
  
  // set default energy fractions values
  fFractionMap.clear();
  fFractionMap[0] = make_pair(0.0, 1.0);
  
  for(i = 0; i < size/2; ++i)
  {
    paramFractions = param[i*2 + 1];
    sizeFractions = paramFractions.GetSize();
  
    ecalFraction = paramFractions[0].GetDouble();
    hcalFraction = paramFractions[1].GetDouble();

    fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
  }
/*
  TFractionMap::iterator itFractionMap;
  for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
  {
    cout << itFractionMap->first << "   " << itFractionMap->second.first  << "   " << itFractionMap->second.second << endl;
  }
*/
  // read resolution formulas  
  fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
  fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));

  // import array with output from other modules
  fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/particles"));
  fItInputArray = fInputArray->MakeIterator();

  // create output arrays
  fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
  fNeutralOutputArray = ExportArray(GetString("NeutralOutputArray", "neutrals"));
  fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
}

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

void Calorimeter::Finish()
{
  vector< vector< Double_t>* >::iterator itPhiBin;
  if(fItInputArray) delete fItInputArray;
  for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
  {
    delete *itPhiBin;
  }

}

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

void Calorimeter::Process()
{
  Candidate *candidate;
  TLorentzVector position, momentum;
  Short_t etaBin, phiBin, number, flags;
  Long64_t towerHit, towerEtaPhi, hitEtaPhi;
  Double_t ecalFraction, hcalFraction;
  Double_t ecalEnergy, hcalEnergy;
  Int_t pdgCode;

  TFractionMap::iterator itFractionMap;

  vector< Double_t >::iterator itEtaBin;
  vector< Double_t >::iterator itPhiBin;
  vector< Double_t > *phiBins;

  vector< Long64_t >::iterator itTowerHits;
  
  DelphesFactory *factory = GetFactory();
  fTowerHits.clear();
  fECalFractions.clear();
  fHCalFractions.clear();

  // loop over all particles
  fItInputArray->Reset();
  number = -1;
  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
  {
    const TLorentzVector &candidatePosition = candidate->Position;
    ++number;

    pdgCode = TMath::Abs(candidate->PID);

    itFractionMap = fFractionMap.find(pdgCode);
    if(itFractionMap == fFractionMap.end())
    {
      itFractionMap = fFractionMap.find(0);
    }

    ecalFraction = itFractionMap->second.first;
    hcalFraction = itFractionMap->second.second;

    fECalFractions.push_back(ecalFraction);
    fHCalFractions.push_back(hcalFraction);
    
    if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
    
    // find eta bin [1, fEtaBins.size - 1]
    itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), candidatePosition.Eta());
    if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
    etaBin = distance(fEtaBins.begin(), itEtaBin);
    
    // phi bins for given eta bin
    phiBins = fPhiBins[etaBin];

    // find phi bin [1, phiBins.size - 1]
    itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), candidatePosition.Phi());
    if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
    phiBin = distance(phiBins->begin(), itPhiBin);

    flags = (candidate->Charge == 0);
    flags |= (pdgCode == 22) << 1;

    // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 16-bits for candidate number, 16-bits for flags}
    towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | Long64_t(number) << 16 | Long64_t(flags);
    
    fTowerHits.push_back(towerHit);
  }

  // all hits are sorted first by eta bin number, then by phi bin number and then by candidate number
  sort(fTowerHits.begin(), fTowerHits.end());

  // loop over all hits
  towerEtaPhi = 0;
  fTower = 0;
  for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
  {
    towerHit = (*itTowerHits);
    flags = (towerHit) & 0x000000000000FFFFLL;
    number = (towerHit >> 16) & 0x000000000000FFFFLL;
    hitEtaPhi = towerHit >> 32;
    candidate = static_cast<Candidate*>(fInputArray->At(number));
    momentum = candidate->Momentum;

    if(towerEtaPhi != hitEtaPhi)
    {
      // switch to next tower
      towerEtaPhi = hitEtaPhi;

      // finalize previous tower
      FinalizeTower();

      // create new tower
      fTower = factory->NewCandidate();

      phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;    
      etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;

      // phi bins for given eta bin
      phiBins = fPhiBins[etaBin];

      // calculate eta and phi of the tower's center
      fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
      fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);

      fTowerECalEnergy = 0.0;
      fTowerHCalEnergy = 0.0;

      fTowerECalNeutralEnergy = 0.0;
      fTowerHCalNeutralEnergy = 0.0;
           
      fTowerHasNeutralHit = kFALSE;
      fTowerHasPhotonHit = kFALSE;
    }

    // fill current tower
    ecalEnergy = momentum.E() * fECalFractions[number];
    hcalEnergy = momentum.E() * fHCalFractions[number];
    
    fTowerECalEnergy += ecalEnergy;
    fTowerHCalEnergy += hcalEnergy;

    fTower->AddCandidate(candidate);   

    // check for neutral hits in current tower
    if(flags & 1)
    {
      fTowerHasNeutralHit = kTRUE;

      fTowerECalNeutralEnergy += ecalEnergy;
      fTowerHCalNeutralEnergy += hcalEnergy;
    }
    
    // check for photon hits in current tower
    if(flags & 2)
    {
      fTowerHasPhotonHit = kTRUE;
    }
  }

  // finalize last tower
  FinalizeTower();
}

//------------------------------------------------------------------------------
void Calorimeter::FinalizeTower()
{
  Candidate *tower;
  Double_t energy, pt;
  Double_t ecalNeutralFraction, hcalNeutralFraction;
  Double_t ecalEnergy, hcalEnergy;

  if(!fTower) return;

  ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
  if(ecalEnergy < 0.0) ecalEnergy = 0.0;

  hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
  if(hcalEnergy < 0.0) hcalEnergy = 0.0;

  energy = ecalEnergy + hcalEnergy;

  if(energy <= 0.0) return;

  pt = energy / TMath::CosH(fTowerEta);
  
  fTower->Position.SetPtEtaPhiE(1.0, fTowerEta, fTowerPhi, 0.0);
  fTower->Momentum.SetPtEtaPhiE(pt, fTowerEta, fTowerPhi, energy);
  fTower->Eem = ecalEnergy;
  fTower->Ehad = hcalEnergy;

  fTowerOutputArray->Add(fTower);
  if(fTowerHasPhotonHit) fPhotonOutputArray->Add(fTower);
  if(fTowerHasNeutralHit)
  {
    tower = static_cast<Candidate*>(fTower->Clone());

    ecalNeutralFraction = (fTowerECalEnergy > 0.0) ? fTowerECalNeutralEnergy / fTowerECalEnergy : 0.0;
    hcalNeutralFraction = (fTowerHCalEnergy > 0.0) ? fTowerHCalNeutralEnergy / fTowerHCalEnergy : 0.0;

    ecalEnergy *= ecalNeutralFraction;
    hcalEnergy *= hcalNeutralFraction;

    energy = ecalEnergy + hcalEnergy;
    pt = energy / TMath::CosH(fTowerEta);

    tower->Position.SetPtEtaPhiE(1.0, fTowerEta, fTowerPhi, 0.0);
    tower->Momentum.SetPtEtaPhiE(pt, fTowerEta, fTowerPhi, energy);
    tower->Eem = ecalEnergy;
    tower->Ehad = hcalEnergy;
    tower->Mother = fTower;

    fNeutralOutputArray->Add(tower);
  }
}

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