Page MenuHomeHEPForge
Files F8310511

No OneTemporary
Actions

View Options

Index: trunk/examples/sartreMain.cpp
===================================================================
--- trunk/examples/sartreMain.cpp (revision 446)
+++ trunk/examples/sartreMain.cpp (revision 447)
@@ -1,297 +1,364 @@
//==============================================================================
// sartreMain.cpp
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2021 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// 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.
// 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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Example main program. Useful to get started and easy to modify to fit
// the user need. Data are stored in a simple but complete Root TTree.
//
//==============================================================================
#include <iostream>
#include <cmath>
#include "TTree.h"
#include "TFile.h"
#include "TGenPhaseSpace.h"
#include "TH1D.h"
#include "TROOT.h"
#include "TLorentzVector.h"
#include "TClonesArray.h"
#include "Sartre.h"
#include "Settings.h"
#include "DipoleModelParameters.h"
#include "TwoBodyVectorMesonDecay.h"
+#include "EicSmearFormatWriter.h"
#define PR(x) cout << #x << " = " << (x) << endl;
+#define EIC_SMEAR_OUTPUT 1
+
using namespace std;
struct rootSartreEvent {
double t;
double Q2;
double x;
double s;
double y;
double W;
double xpom;
int iEvent;
int pol; // 0=transverse or 1=longitudinal
int dmode; // 0=coherent, 1=Incoherent
};
rootSartreEvent myRootSartreEvent;
void randomlyReverseBeams(Event* ); // for UPC mode
int main(int argc, char *argv[])
{
//
// Check command line arguments
//
if (! (argc == 2 || argc == 3) ) {
cout << "Usage: " << argv[0] << " runcard [rootfile]" << endl;
return 2;
}
string runcard = argv[1];
//
// Create the generator and initialize it.
// Once initialized you cannot (should not) change
// the settings w/o re-initialing sartre.
//
Sartre sartre;
bool ok = sartre.init(runcard);
if (!ok) {
cerr << "Error: Initialization of sartre failed." << endl;
return 1;
}
EventGeneratorSettings* settings = sartre.runSettings();
//
// ROOT file
// Use the one from command line unless not given
// in which case the one in the runcard is used.
//
string rootfile;
if (argc == 3) {
rootfile = argv[2];
settings->setRootfile(argv[2]);
}
else
rootfile = settings->rootfile();
+
//
// Print out all Sartre settings
//
settings->list();
TFile *hfile = 0;
if (rootfile.size()) {
hfile = new TFile(rootfile.c_str(),"RECREATE");
cout << "ROOT file is '" << rootfile.c_str() << "'." << endl;
}
-
+
+ //
+ // Derive filename for eic smear compatible output
+ // from root file
+ //
+#if defined(EIC_SMEAR_OUTPUT)
+ string eicSmearFilename;
+ eicSmearFilename = rootfile.substr(0, rootfile.find(".root"));
+ eicSmearFilename += ".out";
+ EicSmearFormatWriter eicSmearWriter;
+ eicSmearWriter.open(eicSmearFilename, settings->enableNuclearBreakup());
+ cout << "eic-smear compatible output file is '" << eicSmearFilename.c_str() << "'." << endl;
+#endif
+
//
// Setup ROOT tree
//
TLorentzVector *eIn = new TLorentzVector;
TLorentzVector *pIn = new TLorentzVector;
TLorentzVector *vm = new TLorentzVector;
TLorentzVector *eOut = new TLorentzVector;
TLorentzVector *pOut = new TLorentzVector;
TLorentzVector *gamma = new TLorentzVector;
TLorentzVector *vmDaughter1 = new TLorentzVector;
TLorentzVector *vmDaughter2 = new TLorentzVector;
TClonesArray protons("TLorentzVector");
TClonesArray neutrons("TLorentzVector");
TClonesArray remnants("TLorentzVector");
TTree tree("tree","sartre");
tree.Branch("event", &myRootSartreEvent.t,
"t/D:Q2/D:x/D:s/D:y/D:W/D:xpom/D:iEvent/I:pol/I:dmode/I");
tree.Branch("eIn", "TLorentzVector", &eIn, 32000, 0);
tree.Branch("pIn", "TLorentzVector", &pIn, 32000, 0);
tree.Branch("vm", "TLorentzVector", &vm, 32000, 0);
tree.Branch("eOut", "TLorentzVector", &eOut, 32000, 0);
tree.Branch("pOut", "TLorentzVector", &pOut, 32000, 0);
tree.Branch("gamma","TLorentzVector", &gamma, 32000, 0);
tree.Branch("vmDaughter1", "TLorentzVector", &vmDaughter1, 32000, 0);
tree.Branch("vmDaughter2", "TLorentzVector", &vmDaughter2, 32000, 0);
-
- if(settings->enableNuclearBreakup()){
+
+ if(settings->enableNuclearBreakup()) {
tree.Branch("protons", &protons);
tree.Branch("neutrons", &neutrons);
tree.Branch("nuclearRemnants", &remnants);
}
//
// Prepare event generation
//
TwoBodyVectorMesonDecay decayEngine;
pair<TLorentzVector, TLorentzVector> daughters;
int daughterID = settings->userInt();
bool doPerformDecay = false;
if (daughterID && settings->vectorMesonId() != 22) {
if (abs(daughterID) == 11 || abs(daughterID) == 13 ||
abs(daughterID) == 211 || abs(daughterID) == 321) {
doPerformDecay = true;
cout << "Will decay vector meson: ";
cout << settings->lookupPDG(settings->vectorMesonId())->GetName();
cout << " -> ";
cout << settings->lookupPDG(daughterID)->GetName();
cout << " ";
cout << settings->lookupPDG(-daughterID)->GetName();
cout << endl;
}
else {
cerr << "Error: Cannot decay vector meson to daughters with ID=" << daughterID << "." << endl;
return 1;
}
}
//
// Events and how often to show status
//
int nPrint;
if (settings->timesToShow())
nPrint = settings->numberOfEvents()/settings->timesToShow();
else
nPrint = 0;
unsigned long maxEvents = settings->numberOfEvents();
cout << "Generating " << maxEvents << " events." << endl << endl;
//
// Event generation
//
for (unsigned long iEvent = 0; iEvent < maxEvents; iEvent++) {
//
// Generate one event
//
Event *event = sartre.generateEvent();
if (nPrint && (iEvent+1)%nPrint == 0 && iEvent != 0) {
cout << "processed " << iEvent+1 << " events" << endl;
}
//
// If Sartre is run in UPC mode, half of the events needs to be
// rotated around and axis perpendicular to z:
//
if(settings->UPC() && settings->A() == settings->UPCA()){
randomlyReverseBeams(event);
}
//
- // Print out (here only for the first few events)
- //
- if (iEvent < 10) event->list();
-
- //
// Fill ROOT tree
//
myRootSartreEvent.iEvent = event->eventNumber;
myRootSartreEvent.t = event->t;
myRootSartreEvent.Q2 = event->Q2;
myRootSartreEvent.x = event->x;
myRootSartreEvent.y = event->y;
myRootSartreEvent.s = event->s;
myRootSartreEvent.W = event->W;
myRootSartreEvent.xpom = event->xpom;
myRootSartreEvent.pol = event->polarization == transverse ? 0 : 1;
myRootSartreEvent.dmode = event->diffractiveMode == coherent ? 0 : 1;
- eIn = &event->particles[0].p;
- pIn = &event->particles[1].p;
- eOut = &event->particles[2].p;
- pOut = &event->particles[6].p;
- vm = &event->particles[4].p;
- gamma = &event->particles[3].p;
+ *eIn = event->particles[0].p;
+ *pIn = event->particles[1].p;
+ *eOut = event->particles[2].p;
+ *pOut = event->particles[6].p;
+ *vm = event->particles[4].p;
+ *gamma = event->particles[3].p;
+ //
+ // Decay the vector meson and fill the decay products in the tree.
+ // For the decay we use the version that is sensitive to the
+ // polarization of the virtual photon (on a statistical basis).
+ // Also update "event" with the decay product.
+ //
+
+ if (doPerformDecay) {
+ daughters = decayEngine.decayVectorMeson(*vm, *event, daughterID);
+ *vmDaughter1 = daughters.first;
+ *vmDaughter2 = daughters.second;
+
+ Particle vmD1, vmD2;
+ vmD1.index = 7;
+ vmD1.pdgId = daughterID;
+ vmD1.status = 1;
+ vmD1.p = daughters.first;
+ vmD1.parents.push_back(4);
+
+ vmD2.index = 8;
+ vmD2.pdgId = -daughterID;
+ vmD2.status = 1;
+ vmD2.p = daughters.second;
+ vmD2.parents.push_back(4);
+
+ event->particles[4].status = 2; // mark J/psi as decayed
+ event->particles[4].daughters.push_back(7);
+ event->particles[4].daughters.push_back(8);
+
+ //
+ // If there are no breakup fragments we simply
+ // add the decay daughters to the particle list.
+ // Otherwise we insert them before the fragments.
+ //
+ if (settings->enableNuclearBreakup() && event->diffractiveMode == incoherent) {
+ vector<Particle>::iterator it;
+ it = event->particles.begin();
+ event->particles.insert(it+7, vmD1);
+ it = event->particles.begin(); // it might have become invalid if reallocated
+ event->particles.insert(it+8, vmD2);
+ }
+ else {
+ event->particles.push_back(vmD1);
+ event->particles.push_back(vmD2);
+ }
+ }
+
//If the event is incoherent, and nuclear breakup is enabled, fill the remnants to the tree
- if(settings->enableNuclearBreakup() and event->diffractiveMode == incoherent){
+ if (settings->enableNuclearBreakup() && event->diffractiveMode == incoherent){
protons.Clear();
neutrons.Clear();
remnants.Clear();
- for(unsigned int iParticle=7; iParticle < event->particles.size(); iParticle++){
+ unsigned int startIndex = doPerformDecay ? 9 : 7;
+ for(unsigned int iParticle=startIndex; iParticle < event->particles.size(); iParticle++){
if(event->particles[iParticle].status == 1) { // Final-state particle
const Particle& particle = event->particles[iParticle];
switch (abs(particle.pdgId)) {
case 2212: // (Anti-)proton
new(protons[protons.GetEntries()]) TLorentzVector(particle.p);
break;
case 2112: // (Anti-)neutron
new(neutrons[neutrons.GetEntries()]) TLorentzVector(particle.p);
break;
default: // Any other remnant
new(remnants[remnants.GetEntries()]) TLorentzVector(particle.p);
break;
} // switch
} // if
} // for
} // if
//
- // Decay the vector meson and fill the decay products in the tree.
- // For the decay we use the version that is sensitive to the
- // polarization of the virtual photon (on a statistical basis).
+ // Print out (here only for the first few events)
+ //
+ // if (iEvent < 10)
+ event->list();
+
+ //
+ // Fill and write event to file
//
- if (doPerformDecay) {
- daughters = decayEngine.decayVectorMeson(*vm, *event, daughterID);
- vmDaughter1 = &daughters.first;
- vmDaughter2 = &daughters.second;
- }
-
tree.Fill();
+#if defined(EIC_SMEAR_OUTPUT)
+ eicSmearWriter.writeEvent(event);
+#endif
}
cout << "All events processed\n" << endl;
//
// That's it, finish up
//
+#if defined(EIC_SMEAR_OUTPUT)
+ eicSmearWriter.close();
+#endif
double totalCS=sartre.totalCrossSection();
TH1D* histoForCSandNumberOfEvents = new TH1D("histoForCSandNumberOfEvents", "Cross-section and Number of Events", 2, 0., 1.);
histoForCSandNumberOfEvents->SetBinContent(1, totalCS);
histoForCSandNumberOfEvents->SetBinContent(2, maxEvents);
double runTime = sartre.runTime();
hfile->Write();
cout << rootfile.c_str() << " written." << endl;
+#if defined(EIC_SMEAR_OUTPUT)
+ cout << eicSmearWriter.filename() << " written." << endl;
+#endif
cout << "Total cross-section: " << totalCS << " nb" << endl;
sartre.listStatus();
cout << "CPU Time/event: " << 1000*runTime/maxEvents << " msec/evt" << endl;
return 0;
}
// UPC only
void randomlyReverseBeams(Event* myEvent)
{
TRandom3 *random = EventGeneratorSettings::randomGenerator();
if(random->Uniform(1) > 0.5){
for(unsigned int i=0; i<myEvent->particles.size(); i++)
myEvent->particles.at(i).p.RotateX(M_PI);
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:42 PM (6 h, 1 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023850
Default Alt Text
(13 KB)

Event Timeline