Fork me on GitHub

source: git/examples/ExternalFastJetBasic.cpp@ 3d10d1f

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 3d10d1f was 3d10d1f, checked in by Michele Selvaggi <michele.selvaggi@…>, 10 years ago

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