[b443089] | 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 | *
|
---|
[b443089] | 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 | *
|
---|
[b443089] | 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 | *
|
---|
[b443089] | 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 |
|
---|
[975405a] | 19 | #include <algorithm>
|
---|
[d7d2da3] | 20 | #include <stdexcept>
|
---|
| 21 | #include <iostream>
|
---|
| 22 | #include <sstream>
|
---|
| 23 | #include <memory>
|
---|
| 24 |
|
---|
| 25 | #include <map>
|
---|
[975405a] | 26 | #include <vector>
|
---|
[d7d2da3] | 27 |
|
---|
| 28 | #include <stdlib.h>
|
---|
| 29 | #include <signal.h>
|
---|
| 30 | #include <stdio.h>
|
---|
| 31 |
|
---|
| 32 | #include "TROOT.h"
|
---|
| 33 | #include "TApplication.h"
|
---|
| 34 |
|
---|
| 35 | #include "TFile.h"
|
---|
| 36 | #include "TObjArray.h"
|
---|
| 37 | #include "TStopwatch.h"
|
---|
| 38 | #include "TDatabasePDG.h"
|
---|
| 39 | #include "TParticlePDG.h"
|
---|
| 40 | #include "TLorentzVector.h"
|
---|
| 41 |
|
---|
| 42 | #include "modules/Delphes.h"
|
---|
| 43 | #include "classes/DelphesStream.h"
|
---|
| 44 | #include "classes/DelphesClasses.h"
|
---|
| 45 | #include "classes/DelphesFactory.h"
|
---|
| 46 |
|
---|
| 47 | #include "ExRootAnalysis/ExRootTreeWriter.h"
|
---|
| 48 | #include "ExRootAnalysis/ExRootTreeBranch.h"
|
---|
| 49 | #include "ExRootAnalysis/ExRootProgressBar.h"
|
---|
| 50 |
|
---|
| 51 | #include "FWCore/FWLite/interface/AutoLibraryLoader.h"
|
---|
| 52 |
|
---|
| 53 | #include "DataFormats/FWLite/interface/Event.h"
|
---|
| 54 | #include "DataFormats/FWLite/interface/Handle.h"
|
---|
| 55 | #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
|
---|
[f29758e] | 56 | #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
|
---|
| 57 | #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
|
---|
[3241a0e] | 58 | #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
|
---|
| 59 | #include "SimDataFormats/GeneratorProducts/interface/WeightsInfo.h"
|
---|
[d7d2da3] | 60 |
|
---|
| 61 | using namespace std;
|
---|
| 62 |
|
---|
| 63 | //---------------------------------------------------------------------------
|
---|
| 64 |
|
---|
[f29758e] | 65 | void ConvertInput(fwlite::Event &event, Long64_t eventCounter,
|
---|
[3241a0e] | 66 | ExRootTreeBranch *branchEvent, ExRootTreeBranch *branchRwgt,
|
---|
| 67 | DelphesFactory *factory, TObjArray *allParticleOutputArray,
|
---|
| 68 | TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray)
|
---|
[d7d2da3] | 69 | {
|
---|
[f29758e] | 70 | fwlite::Handle< GenEventInfoProduct > handleGenEventInfo;
|
---|
| 71 |
|
---|
[3241a0e] | 72 | fwlite::Handle< LHEEventProduct > handleLHEEvent;
|
---|
| 73 |
|
---|
[d7d2da3] | 74 | fwlite::Handle< vector< reco::GenParticle > > handleParticle;
|
---|
[975405a] | 75 | vector< reco::GenParticle >::const_iterator itParticle;
|
---|
| 76 |
|
---|
| 77 | vector< const reco::Candidate * > vectorCandidate;
|
---|
| 78 | vector< const reco::Candidate * >::iterator itCandidate;
|
---|
| 79 |
|
---|
[f29758e] | 80 | handleGenEventInfo.getByLabel(event, "generator");
|
---|
[ebe36ed] | 81 | handleLHEEvent.getByLabel(event, "externalLHEProducer");
|
---|
[d7d2da3] | 82 | handleParticle.getByLabel(event, "genParticles");
|
---|
| 83 |
|
---|
[f29758e] | 84 | HepMCEvent *element;
|
---|
[3241a0e] | 85 | Weight *weight;
|
---|
[d7d2da3] | 86 | Candidate *candidate;
|
---|
| 87 | TDatabasePDG *pdg;
|
---|
| 88 | TParticlePDG *pdgParticle;
|
---|
| 89 | Int_t pdgCode;
|
---|
| 90 |
|
---|
| 91 | Int_t pid, status;
|
---|
[80d4a34] | 92 | Double_t px, py, pz, e, mass;
|
---|
[d7d2da3] | 93 | Double_t x, y, z;
|
---|
| 94 |
|
---|
[3241a0e] | 95 | const vector< gen::WeightsInfo > &vectorWeightsInfo = handleLHEEvent->weights();
|
---|
| 96 | vector< gen::WeightsInfo >::const_iterator itWeightsInfo;
|
---|
| 97 |
|
---|
[f29758e] | 98 | element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
|
---|
| 99 |
|
---|
| 100 | element->Number = eventCounter;
|
---|
| 101 |
|
---|
| 102 | element->ProcessID = handleGenEventInfo->signalProcessID();
|
---|
| 103 | element->MPI = 1;
|
---|
[529fe78] | 104 | element->Weight = handleGenEventInfo->weight();
|
---|
[f29758e] | 105 | element->Scale = handleGenEventInfo->qScale();
|
---|
| 106 | element->AlphaQED = handleGenEventInfo->alphaQED();
|
---|
| 107 | element->AlphaQCD = handleGenEventInfo->alphaQCD();
|
---|
| 108 |
|
---|
| 109 | element->ID1 = 0;
|
---|
| 110 | element->ID2 = 0;
|
---|
| 111 | element->X1 = 0.0;
|
---|
| 112 | element->X2 = 0.0;
|
---|
| 113 | element->ScalePDF = 0.0;
|
---|
| 114 | element->PDF1 = 0.0;
|
---|
| 115 | element->PDF2 = 0.0;
|
---|
| 116 |
|
---|
| 117 | element->ReadTime = 0.0;
|
---|
| 118 | element->ProcTime = 0.0;
|
---|
| 119 |
|
---|
[3241a0e] | 120 | for(itWeightsInfo = vectorWeightsInfo.begin(); itWeightsInfo != vectorWeightsInfo.end(); ++itWeightsInfo)
|
---|
| 121 | {
|
---|
| 122 | weight = static_cast<Weight *>(branchRwgt->NewEntry());
|
---|
| 123 | weight->Weight = itWeightsInfo->wgt;
|
---|
| 124 | }
|
---|
| 125 |
|
---|
[d7d2da3] | 126 | pdg = TDatabasePDG::Instance();
|
---|
| 127 |
|
---|
[975405a] | 128 | for(itParticle = handleParticle->begin(); itParticle != handleParticle->end(); ++itParticle)
|
---|
[d7d2da3] | 129 | {
|
---|
[975405a] | 130 | vectorCandidate.push_back(&*itParticle);
|
---|
| 131 | }
|
---|
| 132 |
|
---|
| 133 | for(itParticle = handleParticle->begin(); itParticle != handleParticle->end(); ++itParticle)
|
---|
| 134 | {
|
---|
| 135 | const reco::GenParticle &particle = *itParticle;
|
---|
[d7d2da3] | 136 |
|
---|
[1e1f73f] | 137 | pid = particle.pdgId();
|
---|
| 138 | status = particle.status();
|
---|
[80d4a34] | 139 | px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.energy(); mass = particle.mass();
|
---|
[1e1f73f] | 140 | x = particle.vx(); y = particle.vy(); z = particle.vz();
|
---|
[d7d2da3] | 141 |
|
---|
[1e1f73f] | 142 | candidate = factory->NewCandidate();
|
---|
[d7d2da3] | 143 |
|
---|
[1e1f73f] | 144 | candidate->PID = pid;
|
---|
| 145 | pdgCode = TMath::Abs(candidate->PID);
|
---|
[d7d2da3] | 146 |
|
---|
[1e1f73f] | 147 | candidate->Status = status;
|
---|
[d7d2da3] | 148 |
|
---|
[f29758e] | 149 | if(particle.mother())
|
---|
| 150 | {
|
---|
| 151 | itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.mother());
|
---|
| 152 | if(itCandidate != vectorCandidate.end()) candidate->M1 = distance(vectorCandidate.begin(), itCandidate);
|
---|
| 153 | }
|
---|
| 154 |
|
---|
[975405a] | 155 | itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(0));
|
---|
| 156 | if(itCandidate != vectorCandidate.end()) candidate->D1 = distance(vectorCandidate.begin(), itCandidate);
|
---|
| 157 |
|
---|
| 158 | itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(particle.numberOfDaughters() - 1));
|
---|
| 159 | if(itCandidate != vectorCandidate.end()) candidate->D2 = distance(vectorCandidate.begin(), itCandidate);
|
---|
| 160 |
|
---|
[1e1f73f] | 161 | pdgParticle = pdg->GetParticle(pid);
|
---|
| 162 | candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
|
---|
[bba4646] | 163 | candidate->Mass = mass;
|
---|
[d7d2da3] | 164 |
|
---|
[1e1f73f] | 165 | candidate->Momentum.SetPxPyPzE(px, py, pz, e);
|
---|
[d7d2da3] | 166 |
|
---|
[1e1f73f] | 167 | candidate->Position.SetXYZT(x, y, z, 0.0);
|
---|
[d7d2da3] | 168 |
|
---|
[1e1f73f] | 169 | allParticleOutputArray->Add(candidate);
|
---|
[d7d2da3] | 170 |
|
---|
[2f82259] | 171 | if(!pdgParticle) continue;
|
---|
[d7d2da3] | 172 |
|
---|
[1e1f73f] | 173 | if(status == 1)
|
---|
| 174 | {
|
---|
| 175 | stableParticleOutputArray->Add(candidate);
|
---|
| 176 | }
|
---|
| 177 | else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
|
---|
| 178 | {
|
---|
| 179 | partonOutputArray->Add(candidate);
|
---|
[d7d2da3] | 180 | }
|
---|
| 181 | }
|
---|
| 182 | }
|
---|
| 183 |
|
---|
| 184 | //---------------------------------------------------------------------------
|
---|
| 185 |
|
---|
| 186 | static bool interrupted = false;
|
---|
| 187 |
|
---|
| 188 | void SignalHandler(int sig)
|
---|
| 189 | {
|
---|
| 190 | interrupted = true;
|
---|
| 191 | }
|
---|
| 192 |
|
---|
| 193 | //---------------------------------------------------------------------------
|
---|
| 194 |
|
---|
| 195 | int main(int argc, char *argv[])
|
---|
| 196 | {
|
---|
| 197 | char appName[] = "DelphesCMSFWLite";
|
---|
| 198 | stringstream message;
|
---|
| 199 | TFile *inputFile = 0;
|
---|
| 200 | TFile *outputFile = 0;
|
---|
| 201 | TStopwatch eventStopWatch;
|
---|
| 202 | ExRootTreeWriter *treeWriter = 0;
|
---|
[3241a0e] | 203 | ExRootTreeBranch *branchEvent = 0, *branchRwgt = 0;
|
---|
[d7d2da3] | 204 | ExRootConfReader *confReader = 0;
|
---|
| 205 | Delphes *modularDelphes = 0;
|
---|
| 206 | DelphesFactory *factory = 0;
|
---|
| 207 | TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
|
---|
| 208 | Int_t i;
|
---|
[5ca3d52] | 209 | Long64_t eventCounter, numberOfEvents;
|
---|
[d7d2da3] | 210 |
|
---|
| 211 | if(argc < 4)
|
---|
| 212 | {
|
---|
| 213 | cout << " Usage: " << appName << " config_file" << " output_file" << " input_file(s)" << endl;
|
---|
| 214 | cout << " config_file - configuration file in Tcl format," << endl;
|
---|
| 215 | cout << " output_file - output file in ROOT format," << endl;
|
---|
| 216 | cout << " input_file(s) - input file(s) in ROOT format." << endl;
|
---|
| 217 | return 1;
|
---|
| 218 | }
|
---|
| 219 |
|
---|
| 220 | signal(SIGINT, SignalHandler);
|
---|
| 221 |
|
---|
| 222 | gROOT->SetBatch();
|
---|
| 223 |
|
---|
| 224 | int appargc = 1;
|
---|
| 225 | char *appargv[] = {appName};
|
---|
| 226 | TApplication app(appName, &appargc, appargv);
|
---|
| 227 |
|
---|
| 228 | AutoLibraryLoader::enable();
|
---|
| 229 |
|
---|
| 230 | try
|
---|
| 231 | {
|
---|
| 232 | outputFile = TFile::Open(argv[2], "CREATE");
|
---|
| 233 |
|
---|
| 234 | if(outputFile == NULL)
|
---|
| 235 | {
|
---|
| 236 | message << "can't open " << argv[2] << endl;
|
---|
| 237 | throw runtime_error(message.str());
|
---|
| 238 | }
|
---|
| 239 |
|
---|
| 240 | treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
|
---|
| 241 |
|
---|
[f29758e] | 242 | branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
|
---|
[3241a0e] | 243 | branchRwgt = treeWriter->NewBranch("Rwgt", Weight::Class());
|
---|
[f29758e] | 244 |
|
---|
[d7d2da3] | 245 | confReader = new ExRootConfReader;
|
---|
| 246 | confReader->ReadFile(argv[1]);
|
---|
| 247 |
|
---|
| 248 | modularDelphes = new Delphes("Delphes");
|
---|
| 249 | modularDelphes->SetConfReader(confReader);
|
---|
| 250 | modularDelphes->SetTreeWriter(treeWriter);
|
---|
| 251 |
|
---|
| 252 | factory = modularDelphes->GetFactory();
|
---|
| 253 | allParticleOutputArray = modularDelphes->ExportArray("allParticles");
|
---|
| 254 | stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
|
---|
| 255 | partonOutputArray = modularDelphes->ExportArray("partons");
|
---|
| 256 |
|
---|
| 257 | modularDelphes->InitTask();
|
---|
| 258 |
|
---|
| 259 | for(i = 3; i < argc && !interrupted; ++i)
|
---|
| 260 | {
|
---|
| 261 | cout << "** Reading " << argv[i] << endl;
|
---|
| 262 |
|
---|
| 263 | inputFile = TFile::Open(argv[i]);
|
---|
| 264 |
|
---|
| 265 | if(inputFile == NULL)
|
---|
| 266 | {
|
---|
| 267 | message << "can't open " << argv[i] << endl;
|
---|
| 268 | throw runtime_error(message.str());
|
---|
| 269 | }
|
---|
| 270 |
|
---|
| 271 | fwlite::Event event(inputFile);
|
---|
| 272 |
|
---|
[5ca3d52] | 273 | numberOfEvents = event.size();
|
---|
[d7d2da3] | 274 |
|
---|
[5ca3d52] | 275 | if(numberOfEvents <= 0) continue;
|
---|
[d7d2da3] | 276 |
|
---|
[a0538b9] | 277 | // ExRootProgressBar progressBar(numberOfEvents - 1);
|
---|
| 278 | ExRootProgressBar progressBar(-1);
|
---|
[d7d2da3] | 279 |
|
---|
| 280 | // Loop over all objects
|
---|
[5ca3d52] | 281 | eventCounter = 0;
|
---|
[d7d2da3] | 282 | modularDelphes->Clear();
|
---|
| 283 | treeWriter->Clear();
|
---|
| 284 | for(event.toBegin(); !event.atEnd() && !interrupted; ++event)
|
---|
| 285 | {
|
---|
[3241a0e] | 286 | ConvertInput(event, eventCounter, branchEvent, branchRwgt, factory,
|
---|
[f29758e] | 287 | allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
|
---|
[d7d2da3] | 288 | modularDelphes->ProcessTask();
|
---|
| 289 |
|
---|
| 290 | treeWriter->Fill();
|
---|
| 291 |
|
---|
| 292 | modularDelphes->Clear();
|
---|
| 293 | treeWriter->Clear();
|
---|
| 294 |
|
---|
[a0538b9] | 295 | progressBar.Update(eventCounter, eventCounter);
|
---|
[5ca3d52] | 296 | ++eventCounter;
|
---|
[d7d2da3] | 297 | }
|
---|
[a0538b9] | 298 |
|
---|
| 299 | progressBar.Update(eventCounter, eventCounter, kTRUE);
|
---|
[d7d2da3] | 300 | progressBar.Finish();
|
---|
| 301 |
|
---|
| 302 | inputFile->Close();
|
---|
| 303 | }
|
---|
| 304 |
|
---|
| 305 | modularDelphes->FinishTask();
|
---|
| 306 | treeWriter->Write();
|
---|
| 307 |
|
---|
| 308 | cout << "** Exiting..." << endl;
|
---|
| 309 |
|
---|
| 310 | delete modularDelphes;
|
---|
| 311 | delete confReader;
|
---|
| 312 | delete treeWriter;
|
---|
| 313 | delete outputFile;
|
---|
| 314 |
|
---|
| 315 | return 0;
|
---|
| 316 | }
|
---|
| 317 | catch(runtime_error &e)
|
---|
| 318 | {
|
---|
| 319 | if(treeWriter) delete treeWriter;
|
---|
| 320 | if(outputFile) delete outputFile;
|
---|
| 321 | cerr << "** ERROR: " << e.what() << endl;
|
---|
| 322 | return 1;
|
---|
| 323 | }
|
---|
| 324 | }
|
---|