Fork me on GitHub

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

Last change on this file since e5ea42e 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: 9.1 KB
RevLine 
[8e602e5]1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
[1fa50c2]4 *
[8e602e5]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.
[1fa50c2]9 *
[8e602e5]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.
[1fa50c2]14 *
[8e602e5]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
[14ae668]19
[caba091]20/*
[14ae668]21########################################################################
22
23
[283537c]24This simple example shows how to use Delphes with an external FastJet installation.
25Events are specified via the multidimensional array "EVENTS" (for an example reading
26an HepMC file see ExternalFastJetHepMC.cpp).
[14ae668]27
[283537c]28In order to run this example you first need to set the paths to your Delphes, FastJet
[14ae668]29and ROOT installations (DELPHES_DIR, FASTJET_DIR and ROOT_DIR):
[a8782e8]30
31DELPHES_DIR=<path to Delphes installation>
32FASTJET_DIR=<path to FastJet installation>
33ROOT_DIR=<path to ROOT installation>
34
[14ae668]35Then run the following commands to build the executable:
[a8782e8]36
[283537c]37DELPHES_LIB="-Wl,-rpath,$DELPHES_DIR -L$DELPHES_DIR -lDelphesNoFastJet"
[a8782e8]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"
[283537c]46LDFLAGS="$FASTJET_LIB $ROOT_LIB -lEG $DELPHES_LIB"
[a8782e8]47
[283537c]48g++ $CXXFLAGS examples/ExternalFastJet/ExternalFastJetBasic.cpp $LDFLAGS -o ExternalFastJetBasic
[14ae668]49
[3d10d1f]50Then run:
[14ae668]51
[283537c]52./ExternalFastJetBasic cards/delphes_card_CMS_NoFastJet.tcl
[14ae668]53
54
55########################################################################
[a8782e8]56*/
57
[8e602e5]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
[14ae668]111
112// this function converts input event array into Delphes candidates (defined below)
113
[8e602e5]114void ConvertInput(Int_t event, DelphesFactory *factory, TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray);
115
[14ae668]116
117//----------------------------------------------------------------------------------------------------------------------------
118
[8e602e5]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
[14ae668]171
172 // fastjet definition
[8e602e5]173 ClusterSequence::print_banner();
174 definition = new JetDefinition(antikt_algorithm, 0.5);
[14ae668]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
[8e602e5]181 inputIterator = inputArray->MakeIterator();
182
183 // Event loop
184 for(event = 0; event < NEVENTS; ++event)
185 {
186 modularDelphes->Clear();
[14ae668]187
188 // convert EVENT input array into Delphes internal format
[8e602e5]189 ConvertInput(event, factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
[14ae668]190
191 // run Delphes reconstruction
[8e602e5]192 modularDelphes->ProcessTask();
193
194 inputList.clear();
195 inputIterator->Reset();
[14ae668]196
197
198 // pass delphes candidates to fastjet clustering
[8e602e5]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 }
[14ae668]205
206 // run clustering
[8e602e5]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
[14ae668]242
243// ------------------------------------------------------------------------------------------------------------------------------------
244
[8e602e5]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.