Fork me on GitHub

source: git/examples/ExternalFastJet/ExternalFastJetBasic.cpp@ aef818f

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

move ExternalFastJet examples to examples/ExternalFastJet
  • Property mode set to 100644
File size: 9.1 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/*
21########################################################################
22
23
24This simple example shows how to use Delphes with an external fastjet installation.
25Events are specified via the multidimentionnal array "EVENTS" (for an example reading
26an hepmc file see ExternalFastJetHepMC.cpp).
27
28In order to run this example you first, you need to set the paths to your Delphes, FastJet
29and ROOT installations (DELPHES_DIR, FASTJET_DIR and ROOT_DIR):
30
31DELPHES_DIR=<path to Delphes installation>
32FASTJET_DIR=<path to FastJet installation>
33ROOT_DIR=<path to ROOT installation>
34
35Then run the following commands to build the executable:
36
37DELPHES_LIB="-Wl,-rpath $DELPHES_DIR -L$DELPHES_DIR -lDelphesNoFastJet"
38
39FASTJET_INC=`$FASTJET_DIR/bin/fastjet-config --cxxflags`
40FASTJET_LIB=`$FASTJET_DIR/bin/fastjet-config --libs`
41
42ROOT_INC=`$ROOT_DIR/bin/root-config --incdir`
43ROOT_LIB=`$ROOT_DIR/bin/root-config --libs`
44
45CXXFLAGS="$FASTJET_INC -I$ROOT_INC -I$DELPHES_DIR -I$DELPHES_DIR/external"
46LDFLAGS="$FASTJET_LIB $ROOT_LIB $DELPHES_LIB"
47
48g++ $CXXFLAGS $LDFLAGS examples/ExternalFastJetBasic.cpp -o examples/ExternalFastJetBasic
49
50Then run:
51
52./examples/ExternalFastJetBasic cards/delphes_card_CMS_NoFastJet.tcl
53
54
55########################################################################
56*/
57
58#include <stdexcept>
59#include <iostream>
60#include <sstream>
61
62#include <stdio.h>
63
64#include "TROOT.h"
65#include "TApplication.h"
66
67#include "TFile.h"
68#include "TObjArray.h"
69#include "TDatabasePDG.h"
70#include "TParticlePDG.h"
71#include "TLorentzVector.h"
72
73#include "modules/Delphes.h"
74#include "classes/DelphesClasses.h"
75#include "classes/DelphesFactory.h"
76
77#include "ExRootAnalysis/ExRootTreeWriter.h"
78#include "ExRootAnalysis/ExRootTreeBranch.h"
79#include "ExRootAnalysis/ExRootProgressBar.h"
80
81#include "fastjet/PseudoJet.hh"
82#include "fastjet/JetDefinition.hh"
83#include "fastjet/ClusterSequence.hh"
84
85using namespace std;
86using namespace fastjet;
87
88//---------------------------------------------------------------------------
89
90const int NEVENTS = 2;
91const int NPARTICLES = 4;
92
93float EVENTS[NEVENTS][NPARTICLES][11] =
94{
95 {
96 {-211, 1, -2.91e-01, 1.99e-01, -1.30e+00, 1.35e+00, 1.3957e-01, 0.0, 0.0, 0.0, 0.0},
97 { 211, 1, -2.62e-02, 1.90e-01, 1.42e+00, 1.44e+00, 1.3957e-01, 0.0, 0.0, 0.0, 0.0},
98 {-211, 1, -3.08e-02, -2.30e-01, 5.13e+00, 5.14e+00, 1.3957e-01, 0.0, 0.0, 0.0, 0.0},
99 { 211, 1, 1.21e-01, -3.52e-01, 2.12e+01, 2.12e+01, 1.3957e-01, 0.0, 0.0, 0.0, 0.0}
100 },
101 {
102 {22, 1, -7.80e-02, -9.70e-02, 2.03e+01, 2.03e+01, 0.0, 0.0, 0.0, 0.0, 0.0},
103 {22, 1, -1.16e-02, -2.08e-02, 4.89e-01, 4.90e-01, 0.0, 0.0, 0.0, 0.0, 0.0},
104 {22, 1, 5.64e-02, -1.34e-02, 5.81e+00, 5.82e+00, 0.0, 0.0, 0.0, 0.0, 0.0},
105 {22, 1, 2.35e-01, -1.04e-01, 9.12e+00, 9.12e+00, 0.0, 0.0, 0.0, 0.0, 0.0}
106 }
107};
108
109//---------------------------------------------------------------------------
110
111
112// this function converts input event array into Delphes candidates (defined below)
113
114void ConvertInput(Int_t event, DelphesFactory *factory, TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray);
115
116
117//----------------------------------------------------------------------------------------------------------------------------
118
119int main(int argc, char *argv[])
120{
121 // Declaration of variables
122 ExRootConfReader *confReader = 0;
123 Delphes *modularDelphes = 0;
124 DelphesFactory *factory = 0;
125 TObjArray *allParticleOutputArray = 0;
126 TObjArray *stableParticleOutputArray = 0;
127 TObjArray *partonOutputArray = 0;
128
129 Int_t event;
130
131 TObjArray *inputArray = 0;
132 TIterator *inputIterator = 0;
133 Candidate *candidate = 0;
134 TLorentzVector momentum;
135
136 JetDefinition *definition = 0;
137 vector<PseudoJet> inputList, outputList;
138 PseudoJet jet;
139
140 gROOT->SetBatch();
141
142 int appargc = 1;
143 char appName[] = "ExternalFastJetBasic";
144 char *appargv[] = {appName};
145 TApplication app(appName, &appargc, appargv);
146
147 if(argc != 2)
148 {
149 cout << " Usage: " << appName << " config_file" << endl;
150 cout << " config_file - configuration file in Tcl format." << endl;
151 return 1;
152 }
153
154 try
155 {
156 // Initialization
157 confReader = new ExRootConfReader;
158 confReader->ReadFile(argv[1]);
159
160 modularDelphes = new Delphes("Delphes");
161 modularDelphes->SetConfReader(confReader);
162
163 factory = modularDelphes->GetFactory();
164
165 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
166 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
167 partonOutputArray = modularDelphes->ExportArray("partons");
168
169 modularDelphes->InitTask();
170
171
172 // fastjet definition
173 ClusterSequence::print_banner();
174 definition = new JetDefinition(antikt_algorithm, 0.5);
175
176 // Define your input candidates to fastjet (by default particle-flow objects).
177 // If you want pure calorimeter towers change "EFlowMerger/eflow" into "Calorimeter/towers":
178
179 inputArray = modularDelphes->ImportArray("EFlowMerger/eflow");
180
181 inputIterator = inputArray->MakeIterator();
182
183 // Event loop
184 for(event = 0; event < NEVENTS; ++event)
185 {
186 modularDelphes->Clear();
187
188 // convert EVENT input array into Delphes internal format
189 ConvertInput(event, factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
190
191 // run Delphes reconstruction
192 modularDelphes->ProcessTask();
193
194 inputList.clear();
195 inputIterator->Reset();
196
197
198 // pass delphes candidates to fastjet clustering
199 while((candidate = static_cast<Candidate*>(inputIterator->Next())))
200 {
201 momentum = candidate->Momentum;
202 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
203 inputList.push_back(jet);
204 }
205
206 // run clustering
207 ClusterSequence sequence(inputList, *definition);
208 outputList.clear();
209 outputList = sorted_by_pt(sequence.inclusive_jets(0.0));
210
211 // tell the user what was done
212 // - the description of the algorithm used
213 // - show the inclusive jets as
214 // {index, rapidity, phi, pt}
215 //----------------------------------------------------------
216 cout << "Ran " << definition->description() << endl;
217
218 // label the columns
219 printf("%5s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt");
220
221 // print out the details for each jet
222 for (unsigned int i = 0; i < outputList.size(); i++) {
223 printf("%5u %15.8f %15.8f %15.8f\n",
224 i, outputList[i].rap(), outputList[i].phi(),
225 outputList[i].perp());
226 }
227 }
228
229 // Finalization
230 modularDelphes->FinishTask();
231 delete modularDelphes;
232 delete confReader;
233 return 0;
234 }
235 catch(runtime_error &e)
236 {
237 cerr << "** ERROR: " << e.what() << endl;
238 return 1;
239 }
240}
241
242
243// ------------------------------------------------------------------------------------------------------------------------------------
244
245void ConvertInput(Int_t event, DelphesFactory *factory, TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray)
246{
247 Int_t particle;
248
249 Candidate *candidate;
250 TDatabasePDG *pdg;
251 TParticlePDG *pdgParticle;
252 Int_t pdgCode;
253
254 Int_t pid, status;
255 Double_t px, py, pz, e, mass;
256 Double_t x, y, z, t;
257
258 pdg = TDatabasePDG::Instance();
259
260 for(particle = 0; particle < NPARTICLES; ++particle)
261 {
262 pid = EVENTS[event][particle][0];
263 status = EVENTS[event][particle][1];
264 px = EVENTS[event][particle][2];
265 py = EVENTS[event][particle][3];
266 pz = EVENTS[event][particle][4];
267 e = EVENTS[event][particle][5];
268 mass = EVENTS[event][particle][6];
269 x = EVENTS[event][particle][7];
270 y = EVENTS[event][particle][8];
271 z = EVENTS[event][particle][9];
272 t = EVENTS[event][particle][10];
273
274 candidate = factory->NewCandidate();
275
276 candidate->PID = pid;
277 pdgCode = TMath::Abs(candidate->PID);
278
279 candidate->Status = status;
280
281 pdgParticle = pdg->GetParticle(pid);
282 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
283 candidate->Mass = mass;
284
285 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
286
287 candidate->Position.SetXYZT(x, y, z, t);
288
289 allParticleOutputArray->Add(candidate);
290
291 if(!pdgParticle) return;
292
293 if(status == 1)
294 {
295 stableParticleOutputArray->Add(candidate);
296 }
297 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
298 {
299 partonOutputArray->Add(candidate);
300 }
301 }
302}
303
Note: See TracBrowser for help on using the repository browser.