Fork me on GitHub

source: git/modules/FastJetGridMedianEstimator.cc@ 0f8390e

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0f8390e was 1fa50c2, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

fix GPLv3 header

  • Property mode set to 100644
File size: 5.7 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class FastJetGridMedianEstimator
21 *
22 *
23 * Computes median energy density per event using a fixed grid.
24 *
25 * \author M. Selvaggi - UCL, Louvain-la-Neuve
26 *
27 */
28
29#include "modules/FastJetGridMedianEstimator.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesFormula.h"
34
35#include "ExRootAnalysis/ExRootResult.h"
36#include "ExRootAnalysis/ExRootFilter.h"
37#include "ExRootAnalysis/ExRootClassifier.h"
38
39#include "TMath.h"
40#include "TString.h"
41#include "TFormula.h"
42#include "TRandom3.h"
43#include "TObjArray.h"
44#include "TDatabasePDG.h"
45#include "TLorentzVector.h"
46
47#include <algorithm>
48#include <stdexcept>
49#include <iostream>
50#include <sstream>
51#include <vector>
52#include <utility>
53
54#include "fastjet/PseudoJet.hh"
55#include "fastjet/JetDefinition.hh"
56#include "fastjet/ClusterSequence.hh"
57#include "fastjet/Selector.hh"
58#include "fastjet/RectangularGrid.hh"
59#include "fastjet/ClusterSequenceArea.hh"
60#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
61
62#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
63
64#include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
65#include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
66#include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
67
68#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
69#include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
70#include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
71#include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
72
73using namespace std;
74using namespace fastjet;
75using namespace fastjet::contrib;
76
77
78//------------------------------------------------------------------------------
79
80FastJetGridMedianEstimator::FastJetGridMedianEstimator() :
81 fItInputArray(0)
82{
83
84}
85
86//------------------------------------------------------------------------------
87
88FastJetGridMedianEstimator::~FastJetGridMedianEstimator()
89{
90
91}
92
93//------------------------------------------------------------------------------
94
95void FastJetGridMedianEstimator::Init()
96{
97 // read eta ranges
98
99 ExRootConfParam param = GetParam("GridRange");
100 Long_t i, size;
101
102 fGrid.clear();
103 size = param.GetSize();
104 for(i = 0; i < size/4; ++i)
105 {
106 fGrid[make_pair(param[i*4].GetDouble(), param[i*4 + 1].GetDouble())] = make_pair(param[i*4 + 2].GetDouble(), param[i*4 + 3].GetDouble());
107 //cout<<param[i*4].GetDouble()<<","<< param[i*4 + 1].GetDouble()<<","<< param[i*4 + 2].GetDouble()<<","<< param[i*4 + 3].GetDouble()<<endl;
108
109 }
110
111
112 //cout<<fGrid[make_pair(0.0,2.5)].first<<","<<fGrid[make_pair(0.0,2.5)].second<<endl;
113
114 // import input array
115
116 fInputArray = ImportArray(GetString("InputArray", "Calorimeter/towers"));
117 fItInputArray = fInputArray->MakeIterator();
118
119 fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
120
121}
122
123//------------------------------------------------------------------------------
124
125void FastJetGridMedianEstimator::Finish()
126{
127 if(fItInputArray) delete fItInputArray;
128}
129
130//------------------------------------------------------------------------------
131
132void FastJetGridMedianEstimator::Process()
133{
134 Candidate *candidate;
135 TLorentzVector momentum;
136 Double_t deta, dphi, detaMin, detaMax;
137 Int_t number;
138 Double_t rho = 0;
139 PseudoJet jet;
140 vector<PseudoJet> inputList, outputList;
141
142 std::map< std::pair< Double_t , Double_t > , std::pair< Double_t , Double_t > >::iterator itGrid;
143
144 DelphesFactory *factory = GetFactory();
145
146 inputList.clear();
147
148 // loop over input objects
149 fItInputArray->Reset();
150 number = 0;
151 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
152 {
153 momentum = candidate->Momentum;
154 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
155 jet.set_user_index(number);
156 inputList.push_back(jet);
157 ++number;
158 }
159
160
161 // compute rho and store it
162
163 // cout<<"caio"<<endl;
164 for(itGrid = fGrid.begin(); itGrid != fGrid.end(); ++itGrid)
165 {
166 //Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
167 // JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
168
169 //cout<<itGrid->first.first<<endl;
170
171 detaMin = (itGrid->first).first;
172 detaMax = (itGrid->first).second;
173 deta = (itGrid->second).first;
174 dphi = (itGrid->second).second;
175
176 //cout<<detaMin<<","<<detaMax<<","<<deta<<","<<dphi<<endl;
177
178
179 RectangularGrid grid(detaMin, detaMax, deta, dphi);
180 //cout<<grid.is_initialised()<<endl;
181 //cout<<grid.rapmin()<<","<<grid.rapmax()<<","<<grid.drap()<<","<<grid.dphi()<<endl;
182
183 GridMedianBackgroundEstimator estimator(grid);
184
185 estimator.set_particles(inputList);
186
187 //cout<<estimator.description()<<endl;
188
189 rho = estimator.rho();
190 //cout<<rho<<endl;
191
192
193 candidate = factory->NewCandidate();
194 candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
195 candidate->Edges[0] = detaMin;
196 candidate->Edges[1] = detaMax;
197 fRhoOutputArray->Add(candidate);
198 }
199
200}
Note: See TracBrowser for help on using the repository browser.