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 | This simple example shows how to use Delphes with an external FastJet installation.
|
---|
25 | Events are specified via the multidimensional array "EVENTS" (for an example reading
|
---|
26 | an HepMC file see ExternalFastJetHepMC.cpp).
|
---|
27 |
|
---|
28 | In order to run this example you first need to set the paths to your Delphes, FastJet
|
---|
29 | and ROOT installations (DELPHES_DIR, FASTJET_DIR and ROOT_DIR):
|
---|
30 |
|
---|
31 | DELPHES_DIR=<path to Delphes installation>
|
---|
32 | FASTJET_DIR=<path to FastJet installation>
|
---|
33 | ROOT_DIR=<path to ROOT installation>
|
---|
34 |
|
---|
35 | Then run the following commands to build the executable:
|
---|
36 |
|
---|
37 | DELPHES_LIB="-Wl,-rpath,$DELPHES_DIR -L$DELPHES_DIR -lDelphesNoFastJet"
|
---|
38 |
|
---|
39 | FASTJET_INC=`$FASTJET_DIR/bin/fastjet-config --cxxflags`
|
---|
40 | FASTJET_LIB=`$FASTJET_DIR/bin/fastjet-config --libs`
|
---|
41 |
|
---|
42 | ROOT_INC=`$ROOT_DIR/bin/root-config --incdir`
|
---|
43 | ROOT_LIB=`$ROOT_DIR/bin/root-config --libs`
|
---|
44 |
|
---|
45 | CXXFLAGS="$FASTJET_INC -I$ROOT_INC -I$DELPHES_DIR -I$DELPHES_DIR/external"
|
---|
46 | LDFLAGS="$FASTJET_LIB $ROOT_LIB -lEG $DELPHES_LIB"
|
---|
47 |
|
---|
48 | g++ $CXXFLAGS examples/ExternalFastJet/ExternalFastJetBasic.cpp $LDFLAGS -o ExternalFastJetBasic
|
---|
49 |
|
---|
50 | Then run:
|
---|
51 |
|
---|
52 | ./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 |
|
---|
85 | using namespace std;
|
---|
86 | using namespace fastjet;
|
---|
87 |
|
---|
88 | //---------------------------------------------------------------------------
|
---|
89 |
|
---|
90 | const int NEVENTS = 2;
|
---|
91 | const int NPARTICLES = 4;
|
---|
92 |
|
---|
93 | float 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 |
|
---|
114 | void ConvertInput(Int_t event, DelphesFactory *factory, TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray);
|
---|
115 |
|
---|
116 |
|
---|
117 | //----------------------------------------------------------------------------------------------------------------------------
|
---|
118 |
|
---|
119 | int 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 |
|
---|
245 | void 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 |
|
---|