/* * Delphes: a framework for fast simulation of a generic collider experiment * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ #include #include #include #include #include #include #include #include #include "TROOT.h" #include "TApplication.h" #include "TFile.h" #include "TObjArray.h" #include "TStopwatch.h" #include "TDatabasePDG.h" #include "TParticlePDG.h" #include "TLorentzVector.h" #include "classes/DelphesStream.h" #include "classes/DelphesClasses.h" #include "classes/DelphesFactory.h" #include "classes/DelphesPileUpReader.h" #include "ExRootAnalysis/ExRootTreeWriter.h" #include "ExRootAnalysis/ExRootTreeBranch.h" #include "ExRootAnalysis/ExRootProgressBar.h" using namespace std; //------------------------------------------------------------------------------ void ProcessEvent(DelphesPileUpReader *reader, ExRootTreeBranch *branch) { GenParticle *particle; Int_t pid; Float_t x, y, z, t; Float_t px, py, pz, e; TDatabasePDG *pdg = TDatabasePDG::Instance(); TParticlePDG *pdgParticle; TLorentzVector momentum; Double_t pt, signPz, cosTheta, eta, rapidity; while(reader->ReadParticle(pid, x, y, z, t, px, py, pz, e)) { particle = static_cast(branch->NewEntry()); particle->PID = pid; particle->X = x; particle->Y = y; particle->Z = z; particle->T = t; particle->Px = px; particle->Py = py; particle->Pz = pz; particle->E = e; particle->Status = 1; particle->IsPU = 1; particle->M1 = -1; particle->M2 = -1; particle->D1 = -1; particle->D2 = -1; pdgParticle = pdg->GetParticle(pid); particle->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999; particle->Mass = pdgParticle ? pdgParticle->Mass() : -999.9; momentum.SetPxPyPzE(px, py, pz, e); pt = momentum.Pt(); cosTheta = TMath::Abs(momentum.CosTheta()); signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0; eta = (cosTheta == 1.0 ? signPz*999.9 : momentum.Eta()); rapidity = (cosTheta == 1.0 ? signPz*999.9 : momentum.Rapidity()); particle->Eta = eta; particle->Phi = momentum.Phi(); particle->PT = pt; particle->Rapidity = rapidity; } } //--------------------------------------------------------------------------- static bool interrupted = false; void SignalHandler(int sig) { interrupted = true; } //--------------------------------------------------------------------------- int main(int argc, char *argv[]) { char appName[] = "pileup2root"; stringstream message; TFile *outputFile = 0; ExRootTreeWriter *treeWriter = 0; ExRootTreeBranch *branchParticle = 0; DelphesPileUpReader *reader = 0; Long64_t entry, allEntries; if(argc != 3) { cout << " Usage: " << appName << " output_file" << " input_file" << endl; cout << " output_file - output file in ROOT format," << endl; cout << " input_file - input binary pile-up file." << endl; return 1; } signal(SIGINT, SignalHandler); gROOT->SetBatch(); int appargc = 1; char *appargv[] = {appName}; TApplication app(appName, &appargc, appargv); try { outputFile = TFile::Open(argv[1], "CREATE"); if(outputFile == NULL) { message << "can't open " << argv[1]; throw runtime_error(message.str()); } cout << "** Reading " << argv[2] << endl; reader = new DelphesPileUpReader(argv[2]); allEntries = reader->GetEntries(); cout << "** Input file contains " << allEntries << " events" << endl; if(allEntries > 0) { treeWriter = new ExRootTreeWriter(outputFile, "Delphes"); branchParticle = treeWriter->NewBranch("Particle", GenParticle::Class()); ExRootProgressBar progressBar(allEntries - 1); // Loop over all events for(entry = 0; entry < allEntries && !interrupted; ++entry) { if(!reader->ReadEntry(entry)) { cerr << "** ERROR: cannot read event " << entry << endl; break; } treeWriter->Clear(); ProcessEvent(reader, branchParticle); treeWriter->Fill(); progressBar.Update(entry); } treeWriter->Write(); progressBar.Finish(); delete treeWriter; } delete reader; cout << "** Exiting..." << endl; return 0; } catch(runtime_error &e) { if(treeWriter) delete treeWriter; if(reader) delete reader; if(outputFile) delete outputFile; cerr << "** ERROR: " << e.what() << endl; return 1; } }