Fork me on GitHub

source: git/examples/ExternalFastJet/ExternalFastJetHepMC.cpp

Last change on this file was 283537c, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

update build instructions for examples/ExternalFastJet

  • Property mode set to 100644
File size: 8.9 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 in HepMC format are read via the DelphesHepMC reader.
26
27In order to run this example you first need to set the paths to your Delphes, FastJet
28and ROOT installations (DELPHES_DIR, FASTJET_DIR and ROOT_DIR):
29
30DELPHES_DIR=<path to Delphes installation>
31FASTJET_DIR=<path to FastJet installation>
32ROOT_DIR=<path to ROOT installation>
33
34Then run the following commands to build the executable:
35
36DELPHES_LIB="-Wl,-rpath,$DELPHES_DIR -L$DELPHES_DIR -lDelphesNoFastJet"
37
38FASTJET_INC=`$FASTJET_DIR/bin/fastjet-config --cxxflags`
39FASTJET_LIB=`$FASTJET_DIR/bin/fastjet-config --libs`
40
41ROOT_INC=`$ROOT_DIR/bin/root-config --incdir`
42ROOT_LIB=`$ROOT_DIR/bin/root-config --libs`
43
44CXXFLAGS="$FASTJET_INC -I$ROOT_INC -I$DELPHES_DIR -I$DELPHES_DIR/external"
45LDFLAGS="$FASTJET_LIB $ROOT_LIB -lEG $DELPHES_LIB"
46
47g++ $CXXFLAGS examples/ExternalFastJet/ExternalFastJetHepMC.cpp $LDFLAGS -o ExternalFastJetHepMC
48
49Then run (you need an event file in HepMC format):
50
51./ExternalFastJetHepMC cards/delphes_card_CMS_NoFastJet.tcl file.hepmc
52
53
54########################################################################
55*/
56
57#include <stdexcept>
58#include <iostream>
59#include <sstream>
60
61#include <vector>
62
63#include <signal.h>
64
65#include "TROOT.h"
66#include "TApplication.h"
67
68#include "TFile.h"
69#include "TObjArray.h"
70#include "TDatabasePDG.h"
71#include "TParticlePDG.h"
72#include "TLorentzVector.h"
73
74#include "modules/Delphes.h"
75#include "classes/DelphesClasses.h"
76#include "classes/DelphesFactory.h"
77#include "classes/DelphesHepMCReader.h"
78
79#include "ExRootAnalysis/ExRootTreeWriter.h"
80#include "ExRootAnalysis/ExRootTreeBranch.h"
81#include "ExRootAnalysis/ExRootProgressBar.h"
82
83#include "fastjet/PseudoJet.hh"
84#include "fastjet/JetDefinition.hh"
85#include "fastjet/ClusterSequence.hh"
86
87// #include "fastjet/contrib/Nsubjettiness.hh"
88// #include "fastjet/contrib/Njettiness.hh"
89// #include "fastjet/contrib/NjettinessPlugin.hh"
90// #include "fastjet/contrib/WinnerTakeAllRecombiner.hh"
91
92using namespace std;
93using namespace fastjet;
94// using namespace fastjet::contrib;
95
96//---------------------------------------------------------------------------
97
98static bool interrupted = false;
99
100void SignalHandler(int sig)
101{
102 interrupted = true;
103}
104
105//---------------------------------------------------------------------------
106
107int main(int argc, char *argv[])
108{
109 char appName[] = "ExternalFastJetHepMC";
110 stringstream message;
111 FILE *inputFile = 0;
112 ExRootConfReader *confReader = 0;
113 Delphes *modularDelphes = 0;
114 DelphesFactory *factory = 0;
115 TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
116 DelphesHepMCReader *reader = 0;
117 Int_t i, maxEvents, skipEvents;
118 Long64_t length, eventCounter;
119
120 TObjArray *inputArray = 0;
121 TIterator *inputIterator = 0;
122 Candidate *candidate = 0;
123 TLorentzVector momentum;
124
125 JetDefinition *definition = 0;
126// JetDefinition::Recombiner *recomb = 0;
127 vector<PseudoJet> inputList, outputList;
128 PseudoJet jet;
129
130 if(argc < 2)
131 {
132 cout << " Usage: " << appName << " config_file" << " [input_file(s)]" << endl;
133 cout << " config_file - configuration file in Tcl format," << endl;
134 cout << " input_file(s) - input file(s) in HepMC format," << endl;
135 cout << " with no input_file, or when input_file is -, read standard input." << endl;
136 return 1;
137 }
138
139 signal(SIGINT, SignalHandler);
140
141 gROOT->SetBatch();
142
143 int appargc = 1;
144 char *appargv[] = {appName};
145 TApplication app(appName, &appargc, appargv);
146
147 try
148 {
149 confReader = new ExRootConfReader;
150 confReader->ReadFile(argv[1]);
151
152 maxEvents = confReader->GetInt("::MaxEvents", 0);
153 skipEvents = confReader->GetInt("::SkipEvents", 0);
154
155 if(maxEvents < 0)
156 {
157 throw runtime_error("MaxEvents must be zero or positive");
158 }
159
160 if(skipEvents < 0)
161 {
162 throw runtime_error("SkipEvents must be zero or positive");
163 }
164
165 modularDelphes = new Delphes("Delphes");
166 modularDelphes->SetConfReader(confReader);
167
168 factory = modularDelphes->GetFactory();
169 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
170 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
171 partonOutputArray = modularDelphes->ExportArray("partons");
172
173 reader = new DelphesHepMCReader;
174
175 modularDelphes->InitTask();
176
177 ClusterSequence::print_banner();
178
179// recomb = new WinnerTakeAllRecombiner();
180// definition = new JetDefinition(antikt_algorithm, 0.5, recomb, Best);
181
182 definition = new JetDefinition(antikt_algorithm, 0.5);
183
184
185 // Define your input candidates to fastjet (by default particle-flow objects).
186 // If you want pure calorimeter towers change "EFlowMerger/eflow" into "Calorimeter/towers":
187
188 inputArray = modularDelphes->ImportArray("EFlowMerger/eflow");
189
190 inputIterator = inputArray->MakeIterator();
191
192
193 // start reading hepmc file
194
195 i = 2;
196 do
197 {
198 if(interrupted) break;
199
200 if(i == argc || strncmp(argv[i], "-", 2) == 0)
201 {
202 cout << "** Reading standard input" << endl;
203 inputFile = stdin;
204 length = -1;
205 }
206 else
207 {
208 cout << "** Reading " << argv[i] << endl;
209 inputFile = fopen(argv[i], "r");
210
211 if(inputFile == NULL)
212 {
213 message << "can't open " << argv[i];
214 throw runtime_error(message.str());
215 }
216
217 fseek(inputFile, 0L, SEEK_END);
218 length = ftello(inputFile);
219 fseek(inputFile, 0L, SEEK_SET);
220
221 if(length <= 0)
222 {
223 fclose(inputFile);
224 ++i;
225 continue;
226 }
227 }
228
229
230 reader->SetInputFile(inputFile);
231
232 // Loop over all objects
233 eventCounter = 0;
234 modularDelphes->Clear();
235 reader->Clear();
236 while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) &&
237 reader->ReadBlock(factory, allParticleOutputArray,
238 stableParticleOutputArray, partonOutputArray) && !interrupted)
239 {
240
241 // loop over events
242 if(reader->EventReady())
243 {
244 ++eventCounter;
245
246 if(eventCounter > skipEvents)
247 {
248
249 // run delphes reconstruction
250 modularDelphes->ProcessTask();
251
252 inputList.clear();
253 inputIterator->Reset();
254
255
256 // pass delphes candidates to fastjet clustering
257 while((candidate = static_cast<Candidate*>(inputIterator->Next())))
258 {
259 momentum = candidate->Momentum;
260 jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
261 inputList.push_back(jet);
262 }
263
264 // run fastjet clustering
265 ClusterSequence sequence(inputList, *definition);
266 outputList.clear();
267 outputList = sorted_by_pt(sequence.inclusive_jets(0.0));
268
269
270 // Prints for the first event:
271 // - the description of the algorithm used
272 // - show the inclusive jets as
273 // {index, rapidity, phi, pt}
274 //----------------------------------------------------------
275
276 if(eventCounter == skipEvents + 1)
277 {
278 cout << "Ran " << definition->description() << endl;
279
280 // label the columns
281 printf("%5s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt");
282
283 // print out the details for each jet
284 for (unsigned int i = 0; i < outputList.size(); i++) {
285 printf("%5u %15.8f %15.8f %15.8f\n",
286 i, outputList[i].rap(), outputList[i].phi(),
287 outputList[i].perp());
288 }
289 }
290 }
291
292 modularDelphes->Clear();
293 reader->Clear();
294 }
295 } // end of event loop
296
297 if(inputFile != stdin) fclose(inputFile);
298
299 ++i;
300 }
301 while(i < argc);
302
303 modularDelphes->FinishTask();
304
305 cout << "** Exiting..." << endl;
306
307 delete definition;
308// delete recomb;
309
310 delete reader;
311 delete modularDelphes;
312 delete confReader;
313
314 return 0;
315 }
316 catch(runtime_error &e)
317 {
318 cerr << "** ERROR: " << e.what() << endl;
319 return 1;
320 }
321}
Note: See TracBrowser for help on using the repository browser.