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

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 3.7 KB
Line 
1
2/** \class EnergySmearing
3 *
4 *  Performs energy resolution smearing.
5 *
6 *  $Date: 2012-11-19 14:05:19 +0100 (Mon, 19 Nov 2012) $
7 *  $Revision: 829 $
8 *
9 *
10 *  \author P. Demin - UCL, Louvain-la-Neuve
11 *
12 */
13
14#include "modules/EnergySmearing.h"
15
16#include "classes/DelphesClasses.h"
17#include "classes/DelphesFactory.h"
18#include "classes/DelphesFormula.h"
19
20#include "ExRootAnalysis/ExRootResult.h"
21#include "ExRootAnalysis/ExRootFilter.h"
22#include "ExRootAnalysis/ExRootClassifier.h"
23
24#include "TMath.h"
25#include "TString.h"
26#include "TFormula.h"
27#include "TRandom3.h"
28#include "TObjArray.h"
29#include "TDatabasePDG.h"
30#include "TLorentzVector.h"
31
32#include <algorithm> 
33#include <stdexcept>
34#include <iostream>
35#include <sstream>
36
37using namespace std;
38
39//------------------------------------------------------------------------------
40
41EnergySmearing::EnergySmearing() :
42  fItInputArray(0)
43{
44}
45
46//------------------------------------------------------------------------------
47
48EnergySmearing::~EnergySmearing()
49{
50}
51
52//------------------------------------------------------------------------------
53
54void EnergySmearing::Init()
55{
56  map< Int_t, DelphesFormula * >::iterator itFormulaMap;
57  ExRootConfParam param;
58  DelphesFormula *formula;
59  Int_t i, size;
60
61  // read resolution formulas
62  param = GetParam("ResolutionFormula");
63  size = param.GetSize();
64 
65  fFormulaMap.clear();
66  for(i = 0; i < size/2; ++i)
67  {
68    formula = new DelphesFormula;
69    formula->Compile(param[i*2 + 1].GetString());
70
71    fFormulaMap[param[i*2].GetInt()] = formula;
72  }
73
74  // set default Resolution formula
75  itFormulaMap = fFormulaMap.find(0);
76  if(itFormulaMap == fFormulaMap.end())
77  {
78    formula = new DelphesFormula;
79    formula->Compile("0.0");
80
81    fFormulaMap[0] = formula;
82  }
83
84  // import array with output from filter/classifier module
85
86  fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/candidates"));
87  fItInputArray = fInputArray->MakeIterator();
88
89  // create output arrays
90
91  fOutputArray = ExportArray(GetString("OutputArray", "candidates"));
92}
93
94//------------------------------------------------------------------------------
95
96void EnergySmearing::Finish()
97{
98  map< Int_t, DelphesFormula * >::iterator itFormulaMap;
99  DelphesFormula *formula;
100 
101  if(fItInputArray) delete fItInputArray;
102 
103  for(itFormulaMap = fFormulaMap.begin(); itFormulaMap != fFormulaMap.end(); ++itFormulaMap)
104  {
105    formula = itFormulaMap->second;
106    if(formula) delete formula;
107  }
108}
109
110//------------------------------------------------------------------------------
111
112void EnergySmearing::Process()
113{
114  Candidate *candidate, *mother;
115  Double_t energy, eta, phi;
116  Int_t pdgCode; 
117  map< Int_t, DelphesFormula * >::iterator itEfficiencyMap, itFormulaMap;
118  DelphesFormula *formula;
119
120  fItInputArray->Reset();
121  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
122  {
123    const TLorentzVector &candidatePosition = candidate->Position;
124    const TLorentzVector &candidateMomentum = candidate->Momentum;
125    pdgCode = TMath::Abs(candidate->PID);
126    eta = candidatePosition.Eta();
127    phi = candidatePosition.Phi();
128    energy = candidateMomentum.E();
129
130    // find smearing formula
131    itFormulaMap = fFormulaMap.find(pdgCode);
132    if(itFormulaMap == fFormulaMap.end())
133    {
134      itFormulaMap = fFormulaMap.find(0);
135    }
136    formula = itFormulaMap->second;
137
138    // apply smearing formula
139    energy = gRandom->Gaus(energy, formula->Eval(0.0, eta, 0.0, energy));
140   
141    if(energy <= 0.0) continue;
142
143    mother = candidate;
144    candidate = static_cast<Candidate*>(candidate->Clone());
145    candidate->Momentum.SetE(energy);
146    candidate->Mother = mother;
147   
148    fOutputArray->Add(candidate);
149  }
150}
151
152//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.