Public Member Functions | Private Member Functions | Private Attributes | List of all members
evgen::NuWroGen Class Reference

A module to check the results from the Monte Carlo generator. More...

Inheritance diagram for evgen::NuWroGen:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 NuWroGen (fhicl::ParameterSet const &pset)
 
virtual ~NuWroGen ()
 
void produce (art::Event &evt)
 
void beginJob ()
 
void beginRun (art::Run &run)
 
void reconfigure (fhicl::ParameterSet const &p)
 
void endJob ()
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Producer & operator= (Producer const &)=delete
 
Producer & operator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifier & operator= (Modifier const &)=delete
 
Modifier & operator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

std::string ParticleStatus (int StatusCode)
 
std::string ReactionChannel (int ccnc, int mode)
 
void FillHistograms (const simb::MCTruth &mc)
 

Private Attributes

std::string fFileName
 
std::string fTreeName
 
std::ifstream * fEventFile
 
TStopwatch fStopwatch
 
std::string fNuWroModuleLabel
 keep track of how long it takes to run the job More...
 
int fEventNumberOffset
 
std::vector< double > fxsecFluxWtd
 
double fxsecTotal
 
TH1F * fGenerated [6]
 Spectra as generated. More...
 
TH1F * fVertexX
 vertex location of generated events in x More...
 
TH1F * fVertexY
 vertex location of generated events in y More...
 
TH1F * fVertexZ
 vertex location of generated events in z More...
 
TH2F * fVertexXY
 vertex location in xy More...
 
TH2F * fVertexXZ
 vertex location in xz More...
 
TH2F * fVertexYZ
 vertex location in yz More...
 
TH1F * fDCosX
 direction cosine in x More...
 
TH1F * fDCosY
 direction cosine in y More...
 
TH1F * fDCosZ
 direction cosine in z More...
 
TH1F * fMuMomentum
 momentum of outgoing muons More...
 
TH1F * fMuDCosX
 direction cosine of outgoing mu in x More...
 
TH1F * fMuDCosY
 direction cosine of outgoing mu in y More...
 
TH1F * fMuDCosZ
 direction cosine of outgoing mu in z More...
 
TH1F * fEMomentum
 momentum of outgoing electrons More...
 
TH1F * fEDCosX
 direction cosine of outgoing e in x More...
 
TH1F * fEDCosY
 direction cosine of outgoing e in y More...
 
TH1F * fEDCosZ
 direction cosine of outgoing e in z More...
 
TH1F * fCCMode
 CC interaction mode. More...
 
TH1F * fNCMode
 CC interaction mode. More...
 
TH1F * fDyn
 mode in detail More...
 
TH1F * fWeight
 NuWro Wt. More...
 
TH1F * fWeightNW
 NuWro Wt. More...
 
TH1F * fDynNew
 
TH1F * fDynNewThresh
 
TH2F * f2DynNew
 
TH2F * f2DynNewThresh
 
TH1F * fECons
 histogram to determine if energy is conserved in the event More...
 
unsigned int countFile
 
event * NuWroTTree
 
TChain * ch
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

A module to check the results from the Monte Carlo generator.

Definition at line 66 of file NuWroGen_module.cc.

Constructor & Destructor Documentation

evgen::NuWroGen::NuWroGen ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 1252 of file NuWroGen_module.cc.

1253  : EDProducer{pset}
1254  {
1255 
1256  fEventNumberOffset = pset.get< int >("EventNumberOffset",0);
1257 
1258  this->reconfigure(pset);
1259  fStopwatch.Start();
1260 
1261  produces< std::vector<simb::MCTruth> >();
1262  produces< sumdata::RunData, art::InRun >();
1263 
1264  }
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
evgen::NuWroGen::fStopwatch
TStopwatch fStopwatch
Definition: NuWroGen_module.cc:88
evgen::NuWroGen::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: NuWroGen_module.cc:1274
evgen::NuWroGen::~NuWroGen ( )
virtual

Definition at line 1268 of file NuWroGen_module.cc.

1269  {
1270  fStopwatch.Stop();
1271  }
evgen::NuWroGen::fStopwatch
TStopwatch fStopwatch
Definition: NuWroGen_module.cc:88

Member Function Documentation

void evgen::NuWroGen::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 1283 of file NuWroGen_module.cc.

1283  {
1284 
1285  // Get access to the TFile service.
1286  art::ServiceHandle<art::TFileService const> tfs;
1287 
1288  fGenerated[0] = tfs->make<TH1F>("fGenerated_necc","", 100, 0.0, 20.0);
1289  fGenerated[1] = tfs->make<TH1F>("fGenerated_nebcc","", 100, 0.0, 20.0);
1290  fGenerated[2] = tfs->make<TH1F>("fGenerated_nmcc","", 100, 0.0, 20.0);
1291  fGenerated[3] = tfs->make<TH1F>("fGenerated_nmbcc","", 100, 0.0, 20.0);
1292  fGenerated[4] = tfs->make<TH1F>("fGenerated_nnc","", 100, 0.0, 20.0);
1293  fGenerated[5] = tfs->make<TH1F>("fGenerated_nbnc","", 100, 0.0, 20.0);
1294 
1295  fDCosX = tfs->make<TH1F>("fDCosX", ";dx/ds", 200, -1., 1.);
1296  fDCosY = tfs->make<TH1F>("fDCosY", ";dy/ds", 200, -1., 1.);
1297  fDCosZ = tfs->make<TH1F>("fDCosZ", ";dz/ds", 200, -1., 1.);
1298 
1299  fMuMomentum = tfs->make<TH1F>("fMuMomentum", ";p_{#mu} (GeV/c)", 500, 0., 50.);
1300  fMuDCosX = tfs->make<TH1F>("fMuDCosX", ";dx/ds;", 200, -1., 1.);
1301  fMuDCosY = tfs->make<TH1F>("fMuDCosY", ";dy/ds;", 200, -1., 1.);
1302  fMuDCosZ = tfs->make<TH1F>("fMuDCosZ", ";dz/ds;", 200, -1., 1.);
1303 
1304  fEMomentum = tfs->make<TH1F>("fEMomentum", ";p_{e} (GeV/c)", 500, 0., 50.);
1305  fEDCosX = tfs->make<TH1F>("fEDCosX", ";dx/ds;", 200, -1., 1.);
1306  fEDCosY = tfs->make<TH1F>("fEDCosY", ";dy/ds;", 200, -1., 1.);
1307  fEDCosZ = tfs->make<TH1F>("fEDCosZ", ";dz/ds;", 200, -1., 1.);
1308 
1309  fCCMode = tfs->make<TH1F>("fCCMode", ";CC Interaction Mode;", 5, 0., 5.);
1310  fCCMode->GetXaxis()->SetBinLabel(1, "QE");
1311  fCCMode->GetXaxis()->SetBinLabel(2, "Res");
1312  fCCMode->GetXaxis()->SetBinLabel(3, "DIS");
1313  fCCMode->GetXaxis()->SetBinLabel(4, "Coh");
1314  fCCMode->GetXaxis()->SetBinLabel(5, "kInverseMuDecay");
1315  fCCMode->GetXaxis()->CenterLabels();
1316 
1317  fNCMode = tfs->make<TH1F>("fNCMode", ";NC Interaction Mode;", 5, 0., 5.);
1318  fNCMode->GetXaxis()->SetBinLabel(1, "QE");
1319  fNCMode->GetXaxis()->SetBinLabel(2, "Res");
1320  fNCMode->GetXaxis()->SetBinLabel(3, "DIS");
1321  fNCMode->GetXaxis()->SetBinLabel(4, "Coh");
1322  fNCMode->GetXaxis()->SetBinLabel(5, "kNuElectronElastic");
1323  fNCMode->GetXaxis()->CenterLabels();
1324 
1325  fWeight = tfs->make<TH1F>("fWeight", ";Weight1;", 20, 1.E-42, 1.E-38);
1326  fWeightNW = tfs->make<TH1F>("fWeightNW", ";Weight2;", 20, 1.E-42, 1.E-38);
1327 
1328  fDyn = tfs->make<TH1F>("fDyn", ";Canonical Interaction Mode;", 13, 0., 13.);
1329  fDyn->GetXaxis()->SetBinLabel(1, "CCQE");
1330  fDyn->GetXaxis()->SetBinLabel(2, "NCelastic");
1331  fDyn->GetXaxis()->SetBinLabel(3, "CCResp2ppi+");
1332  fDyn->GetXaxis()->SetBinLabel(4, "CCResn2ppi0");
1333  fDyn->GetXaxis()->SetBinLabel(5, "CCResn2npi+");
1334  fDyn->GetXaxis()->SetBinLabel(6, "NCResp2ppi0");
1335  fDyn->GetXaxis()->SetBinLabel(7, "NCResp2npi+");
1336  fDyn->GetXaxis()->SetBinLabel(8, "NCResn2npi0");
1337  fDyn->GetXaxis()->SetBinLabel(9, "NCResn2ppi-");
1338  fDyn->GetXaxis()->SetBinLabel(10, "CC-DIS");
1339  fDyn->GetXaxis()->SetBinLabel(11, "NC-DIS");
1340  fDyn->GetXaxis()->SetBinLabel(12, "NC-COH");
1341  fDyn->GetXaxis()->SetBinLabel(13, "CC-COH");
1342  fDyn->GetXaxis()->CenterLabels();
1343 
1344  fDynNew = tfs->make<TH1F>("fDynNew", ";New Style Accounting Mode;", 10, 0., 10.);
1345  fDynNew->GetXaxis()->SetBinLabel(1, "1mu0p0pi");
1346  fDynNew->GetXaxis()->SetBinLabel(2, "1mu1p0pi");
1347  fDynNew->GetXaxis()->SetBinLabel(3, "1muge2p0pi");
1348  fDynNew->GetXaxis()->SetBinLabel(4, "1mu0p1pi");
1349  fDynNew->GetXaxis()->SetBinLabel(5, "1mu1p1pi");
1350  fDynNew->GetXaxis()->SetBinLabel(6, "1muge2p1pi");
1351  fDynNew->GetXaxis()->SetBinLabel(7, "1mu0p2pi");
1352  fDynNew->GetXaxis()->SetBinLabel(8, "1muge1p2pi");
1353  fDynNew->GetXaxis()->SetBinLabel(9, "NC");
1354  fDynNew->GetXaxis()->SetBinLabel(10, "Other");
1355  fDynNew->GetXaxis()->CenterLabels();
1356 
1357  fDynNewThresh = tfs->make<TH1F>("fDynNewThresh", ";New Style Accounting Mode (Tp>50MeV);", 10, 0., 10.);
1358  fDynNewThresh->GetXaxis()->SetBinLabel(1, "1mu0p0pi");
1359  fDynNewThresh->GetXaxis()->SetBinLabel(2, "1mu1p0pi");
1360  fDynNewThresh->GetXaxis()->SetBinLabel(3, "1muge2p0pi");
1361  fDynNewThresh->GetXaxis()->SetBinLabel(4, "1mu0p1pi");
1362  fDynNewThresh->GetXaxis()->SetBinLabel(5, "1mu1p1pi");
1363  fDynNewThresh->GetXaxis()->SetBinLabel(6, "1muge2p1pi");
1364  fDynNewThresh->GetXaxis()->SetBinLabel(7, "1mu0p2pi");
1365  fDynNewThresh->GetXaxis()->SetBinLabel(8, "1muge1p2pi");
1366  fDynNewThresh->GetXaxis()->SetBinLabel(9, "NC");
1367  fDynNewThresh->GetXaxis()->SetBinLabel(10, "Other");
1368  fDynNewThresh->GetXaxis()->CenterLabels();
1369 
1370  f2DynNew = tfs->make<TH2F>("f2DynNew", ";Old vs New Style Accounting Mode;", 13, 0.,13., 10,0.,10.);
1371 
1372  f2DynNew->GetXaxis()->SetBinLabel(1, "CCQE");
1373  f2DynNew->GetXaxis()->SetBinLabel(2, "NCelastic");
1374  f2DynNew->GetXaxis()->SetBinLabel(3, "CCResp2ppi+");
1375  f2DynNew->GetXaxis()->SetBinLabel(4, "CCResn2ppi0");
1376  f2DynNew->GetXaxis()->SetBinLabel(5, "CCResn2npi+");
1377  f2DynNew->GetXaxis()->SetBinLabel(6, "NCResp2ppi0");
1378  f2DynNew->GetXaxis()->SetBinLabel(7, "NCResp2npi+");
1379  f2DynNew->GetXaxis()->SetBinLabel(8, "NCResn2npi0");
1380  f2DynNew->GetXaxis()->SetBinLabel(9, "NCResn2ppi-");
1381  f2DynNew->GetXaxis()->SetBinLabel(10, "CC-DIS");
1382  f2DynNew->GetXaxis()->SetBinLabel(11, "NC-DIS");
1383  f2DynNew->GetXaxis()->SetBinLabel(12, "NC-COH");
1384  f2DynNew->GetXaxis()->SetBinLabel(13, "CC-COH");
1385  f2DynNew->GetXaxis()->CenterLabels();
1386  f2DynNew->GetYaxis()->SetBinLabel(1, "1mu0p0pi");
1387  f2DynNew->GetYaxis()->SetBinLabel(2, "1mu1p0pi");
1388  f2DynNew->GetYaxis()->SetBinLabel(3, "1muge2p0pi");
1389  f2DynNew->GetYaxis()->SetBinLabel(4, "1mu0p1pi");
1390  f2DynNew->GetYaxis()->SetBinLabel(5, "1mu1p1pi");
1391  f2DynNew->GetYaxis()->SetBinLabel(6, "1muge2p1pi");
1392  f2DynNew->GetYaxis()->SetBinLabel(7, "1mu0p2pi");
1393  f2DynNew->GetYaxis()->SetBinLabel(8, "1muge1p2pi");
1394  f2DynNew->GetYaxis()->SetBinLabel(9, "NC");
1395  f2DynNew->GetYaxis()->SetBinLabel(10, "Other");
1396  f2DynNew->GetYaxis()->CenterLabels();
1397 
1398  f2DynNewThresh = tfs->make<TH2F>("f2DynNewThresh", ";Old vs New Style Accounting Mode;", 13, 0.,13., 10, 0.,10.);
1399 
1400  f2DynNewThresh->GetXaxis()->SetBinLabel(1, "CCQE");
1401  f2DynNewThresh->GetXaxis()->SetBinLabel(2, "NCelastic");
1402  f2DynNewThresh->GetXaxis()->SetBinLabel(3, "CCResp2ppi+");
1403  f2DynNewThresh->GetXaxis()->SetBinLabel(4, "CCResn2ppi0");
1404  f2DynNewThresh->GetXaxis()->SetBinLabel(5, "CCResn2npi+");
1405  f2DynNewThresh->GetXaxis()->SetBinLabel(6, "NCResp2ppi0");
1406  f2DynNewThresh->GetXaxis()->SetBinLabel(7, "NCResp2npi+");
1407  f2DynNewThresh->GetXaxis()->SetBinLabel(8, "NCResn2npi0");
1408  f2DynNewThresh->GetXaxis()->SetBinLabel(9, "NCResn2ppi-");
1409  f2DynNewThresh->GetXaxis()->SetBinLabel(10, "CC-DIS");
1410  f2DynNewThresh->GetXaxis()->SetBinLabel(11, "NC-DIS");
1411  f2DynNewThresh->GetXaxis()->SetBinLabel(12, "NC-COH");
1412  f2DynNewThresh->GetXaxis()->SetBinLabel(13, "CC-COH");
1413  f2DynNewThresh->GetXaxis()->CenterLabels();
1414  f2DynNewThresh->GetYaxis()->SetBinLabel(1, "1mu0p0pi");
1415  f2DynNewThresh->GetYaxis()->SetBinLabel(2, "1mu1p0pi");
1416  f2DynNewThresh->GetYaxis()->SetBinLabel(3, "1muge2p0pi");
1417  f2DynNewThresh->GetYaxis()->SetBinLabel(4, "1mu0p1pi");
1418  f2DynNewThresh->GetYaxis()->SetBinLabel(5, "1mu1p1pi");
1419  f2DynNewThresh->GetYaxis()->SetBinLabel(6, "1muge2p1pi");
1420  f2DynNewThresh->GetYaxis()->SetBinLabel(7, "1mu0p2pi");
1421  f2DynNewThresh->GetYaxis()->SetBinLabel(8, "1muge1p2pi");
1422  f2DynNewThresh->GetYaxis()->SetBinLabel(9, "NC");
1423  f2DynNewThresh->GetYaxis()->SetBinLabel(10, "Other");
1424  f2DynNewThresh->GetYaxis()->CenterLabels();
1425 
1426  //fDeltaE = tfs->make<TH1F>("fDeltaE", ";#Delta E_{#nu} (GeV);", 200, -1., 1.);
1427  fECons = tfs->make<TH1F>("fECons", ";#Delta E(#nu,lepton);", 500, -5., 5.);
1428 
1429  art::ServiceHandle<geo::Geometry const> geo;
1430  double x = 2.1*geo->DetHalfWidth();
1431  double y = 2.1*geo->DetHalfHeight();
1432  double z = 2.*geo->DetLength();
1433  int xdiv = TMath::Nint(2*x/5.);
1434  int ydiv = TMath::Nint(2*y/5.);
1435  int zdiv = TMath::Nint(2*z/5.);
1436 
1437  fVertexX = tfs->make<TH1F>("fVertexX", ";x (cm)", xdiv, -x, x);
1438  fVertexY = tfs->make<TH1F>("fVertexY", ";y (cm)", ydiv, -y, y);
1439  fVertexZ = tfs->make<TH1F>("fVertexZ", ";z (cm)", zdiv, -0.2*z, z);
1440 
1441  fVertexXY = tfs->make<TH2F>("fVertexXY", ";x (cm);y (cm)", xdiv, -x, x, ydiv, -y, y);
1442  fVertexXZ = tfs->make<TH2F>("fVertexXZ", ";z (cm);x (cm)", zdiv, -0.2*z, z, xdiv, -x, x);
1443  fVertexYZ = tfs->make<TH2F>("fVertexYZ", ";z (cm);y (cm)", zdiv, -0.2*z, z, ydiv, -y, y);
1444 
1445 
1446  TClass::GetClass("line")->GetStreamerInfo(1);
1447  ch = new TChain(fTreeName.c_str());
1448  ch->Add(fFileName.c_str());
1449 
1450 
1451  std::cout << " Num entries in TTree is " << ch->GetEntries() << std::endl;
1452 
1453  NuWroTTree = new event();
1454  ch->SetBranchAddress ("e", &NuWroTTree);
1455 
1456  // initiate flux-wtd XSections.
1457  std::cout << "NuWroGen: Here's the output of the .txt file" << std::endl;
1458  std::ifstream xsecTxtFile((fFileName+".txt").c_str());
1459  unsigned int cntline(0);
1460  std::string line;
1461  if (xsecTxtFile.is_open())
1462  {
1463  while ( xsecTxtFile.good() )
1464  {
1465  getline (xsecTxtFile,line);
1466  cout << line << endl;
1467  if (cntline==0) { cntline++; continue;} // first line is header.
1468  stringstream ss(line); // Insert the string into a stream
1469  vector<std::string> tokens; // Create vector to hold our words
1470  string buf;
1471  while (ss >> buf)
1472  {
1473  tokens.push_back(buf);
1474  }
1475  // want last element in line. That's the xsec.
1476  if (tokens.size() && line.length())
1477  fxsecFluxWtd.push_back(atof(tokens.back().c_str()));
1478  else xsecTxtFile.close();
1479  tokens.clear();
1480  }
1481  }
1482  else std::cout << "Unable to open file";
1483 
1484  fxsecTotal=0.0;
1485  for (unsigned int ii=0;ii<fxsecFluxWtd.size();ii++)
1486  {
1487  fxsecTotal+=fxsecFluxWtd.at(ii);
1488  }
1489 
1490  countFile = fEventNumberOffset;
1491  std::cout << " Start this job on event " << countFile << std::endl;
1492  }
evgen::NuWroGen::fEDCosZ
TH1F * fEDCosZ
direction cosine of outgoing e in z
Definition: NuWroGen_module.cc:118
geo::GeometryCore::DetHalfWidth
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
Definition: GeometryCore.cxx:752
evgen::NuWroGen::fVertexYZ
TH2F * fVertexYZ
vertex location in yz
Definition: NuWroGen_module.cc:104
string
std::string string
Definition: nybbler.cc:12
evgen::NuWroGen::fCCMode
TH1F * fCCMode
CC interaction mode.
Definition: NuWroGen_module.cc:120
evgen::NuWroGen::fVertexXY
TH2F * fVertexXY
vertex location in xy
Definition: NuWroGen_module.cc:102
evgen::NuWroGen::fWeightNW
TH1F * fWeightNW
NuWro Wt.
Definition: NuWroGen_module.cc:124
evgen::NuWroGen::fDCosX
TH1F * fDCosX
direction cosine in x
Definition: NuWroGen_module.cc:106
evgen::NuWroGen::fGenerated
TH1F * fGenerated[6]
Spectra as generated.
Definition: NuWroGen_module.cc:96
evgen::NuWroGen::fVertexXZ
TH2F * fVertexXZ
vertex location in xz
Definition: NuWroGen_module.cc:103
evgen::NuWroGen::fNCMode
TH1F * fNCMode
CC interaction mode.
Definition: NuWroGen_module.cc:121
evgen::NuWroGen::fEDCosY
TH1F * fEDCosY
direction cosine of outgoing e in y
Definition: NuWroGen_module.cc:117
geo::GeometryCore::DetHalfHeight
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
Definition: GeometryCore.cxx:758
evgen::NuWroGen::fTreeName
std::string fTreeName
Definition: NuWroGen_module.cc:85
vector::size
size_t size
Definition: lodepng.cpp:55
evgen::NuWroGen::fWeight
TH1F * fWeight
NuWro Wt.
Definition: NuWroGen_module.cc:123
geo::GeometryCore::DetLength
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
Definition: GeometryCore.cxx:764
evgen::NuWroGen::fVertexY
TH1F * fVertexY
vertex location of generated events in y
Definition: NuWroGen_module.cc:99
evgen::NuWroGen::fVertexZ
TH1F * fVertexZ
vertex location of generated events in z
Definition: NuWroGen_module.cc:100
evgen::NuWroGen::fDCosY
TH1F * fDCosY
direction cosine in y
Definition: NuWroGen_module.cc:107
evgen::NuWroGen::fECons
TH1F * fECons
histogram to determine if energy is conserved in the event
Definition: NuWroGen_module.cc:132
evgen::NuWroGen::fEDCosX
TH1F * fEDCosX
direction cosine of outgoing e in x
Definition: NuWroGen_module.cc:116
evgen::NuWroGen::fDCosZ
TH1F * fDCosZ
direction cosine in z
Definition: NuWroGen_module.cc:108
pduneana::line
void line(double t, double *p, double &x, double &y, double &z)
Definition: PDSPAnalyzer_module.cc:4741
evgen::NuWroGen::fxsecFluxWtd
std::vector< double > fxsecFluxWtd
Definition: NuWroGen_module.cc:93
train.x
list x
Definition: train.py:276
evgen::NuWroGen::fMuDCosX
TH1F * fMuDCosX
direction cosine of outgoing mu in x
Definition: NuWroGen_module.cc:111
evgen::NuWroGen::fDyn
TH1F * fDyn
mode in detail
Definition: NuWroGen_module.cc:122
evgen::NuWroGen::fFileName
std::string fFileName
Definition: NuWroGen_module.cc:84
evgen::NuWroGen::fMuDCosY
TH1F * fMuDCosY
direction cosine of outgoing mu in y
Definition: NuWroGen_module.cc:112
evgen::NuWroGen::countFile
unsigned int countFile
Definition: NuWroGen_module.cc:136
evgen::NuWroGen::fEMomentum
TH1F * fEMomentum
momentum of outgoing electrons
Definition: NuWroGen_module.cc:115
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
evgen::NuWroGen::fMuMomentum
TH1F * fMuMomentum
momentum of outgoing muons
Definition: NuWroGen_module.cc:110
evgen::NuWroGen::fVertexX
TH1F * fVertexX
vertex location of generated events in x
Definition: NuWroGen_module.cc:98
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
evgen::NuWroGen::fMuDCosZ
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
Definition: NuWroGen_module.cc:113
void evgen::NuWroGen::beginRun ( art::Run &  run)
virtual

Reimplemented from art::EDProducer.

Definition at line 1495 of file NuWroGen_module.cc.

1496  {
1497  art::ServiceHandle<geo::Geometry const> geo;
1498  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
1499  }
geo::GeometryCore::DetectorName
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:1761
art::DataViewImpl::put
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
void evgen::NuWroGen::endJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 1502 of file NuWroGen_module.cc.

1503  {
1504  delete NuWroTTree;
1505  delete ch;
1506  fxsecFluxWtd.clear();
1507  }
evgen::NuWroGen::fxsecFluxWtd
std::vector< double > fxsecFluxWtd
Definition: NuWroGen_module.cc:93
void evgen::NuWroGen::FillHistograms ( const simb::MCTruth &  mc)
private

look for the outgoing lepton in the particle stack just interested in the first one

Definition at line 1752 of file NuWroGen_module.cc.

1753  {
1754  // Decide which histograms to put the spectrum in
1755  int id = -1;
1756  if (mc.GetNeutrino().CCNC()==simb::kCC) {
1757  fCCMode->Fill(mc.GetNeutrino().Mode());
1758  if (mc.GetNeutrino().Nu().PdgCode() == 12) id = 0;
1759  else if (mc.GetNeutrino().Nu().PdgCode() == -12) id = 1;
1760  else if (mc.GetNeutrino().Nu().PdgCode() == 14) id = 2;
1761  else if (mc.GetNeutrino().Nu().PdgCode() == -14) id = 3;
1762  else return;
1763  }
1764  else {
1765  fNCMode->Fill(mc.GetNeutrino().Mode());
1766  if (mc.GetNeutrino().Nu().PdgCode() > 0) id = 4;
1767  else id = 5;
1768  }
1769  if (id==-1) abort();
1770 
1771  // Fill the specta histograms
1772  fGenerated[id]->Fill(mc.GetNeutrino().Nu().E() );
1773 
1774  //< fill the vertex histograms from the neutrino - that is always
1775  //< particle 0 in the list
1776  simb::MCNeutrino mcnu = mc.GetNeutrino();
1777  const simb::MCParticle nu = mcnu.Nu();
1778 
1779  fVertexX->Fill(nu.Vx());
1780  fVertexY->Fill(nu.Vy());
1781  fVertexZ->Fill(nu.Vz());
1782 
1783  fVertexXY->Fill(nu.Vx(), nu.Vy());
1784  fVertexXZ->Fill(nu.Vz(), nu.Vx());
1785  fVertexYZ->Fill(nu.Vz(), nu.Vy());
1786 
1787  double mom = nu.P();
1788  if(std::abs(mom) > 0.){
1789  fDCosX->Fill(nu.Px()/mom);
1790  fDCosY->Fill(nu.Py()/mom);
1791  fDCosZ->Fill(nu.Pz()/mom);
1792  }
1793 
1794 
1795 // MF_LOG_DEBUG("GENIEInteractionInformation")
1796 // << std::endl
1797 // << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode())
1798 // << std::endl
1799 // << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl
1800 // << std::setiosflags(std::ios::left)
1801 // << std::setw(20) << "PARTICLE"
1802 // << std::setiosflags(std::ios::left)
1803 // << std::setw(32) << "STATUS"
1804 // << std::setw(18) << "E (GeV)"
1805 // << std::setw(18) << "m (GeV/c2)"
1806 // << std::setw(18) << "Ek (GeV)"
1807 // << std::endl << std::endl;
1808 
1809 // const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
1810 
1811 // // Loop over the particle stack for this event
1812 // for(int i = 0; i < mc.NParticles(); ++i){
1813 // simb::MCParticle part(mc.GetParticle(i));
1814 // std::string name = databasePDG->GetParticle(part.PdgCode())->GetName();
1815 // int code = part.StatusCode();
1816 // std::string status = ParticleStatus(code);
1817 // double mass = part.Mass();
1818 // double energy = part.E();
1819 // double Ek = (energy-mass); // Kinetic Energy (GeV)
1820 // if(status=="kIStFinalStB4Interactions")
1821 // MF_LOG_DEBUG("GENIEFinalState")
1822 // << std::setiosflags(std::ios::left) << std::setw(20) << name
1823 // << std::setiosflags(std::ios::left) << std::setw(32) <<status
1824 // << std::setw(18)<< energy
1825 // << std::setw(18)<< mass
1826 // << std::setw(18)<< Ek <<std::endl;
1827 // else
1828 // MF_LOG_DEBUG("GENIEFinalState")
1829 // << std::setiosflags(std::ios::left) << std::setw(20) << name
1830 // << std::setiosflags(std::ios::left) << std::setw(32) << status
1831 // << std::setw(18) << energy
1832 // << std::setw(18) << mass <<std::endl;
1833 
1834  std::cout << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode()) << std::endl;
1835  std::cout << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl;
1836  std::cout << std::setiosflags(std::ios::left)
1837  << std::setw(20) << "PARTICLE"
1838  << std::setiosflags(std::ios::left)
1839  << std::setw(32) << "STATUS"
1840  << std::setw(18) << "E (GeV)"
1841  << std::setw(18) << "m (GeV/c2)"
1842  << std::setw(18) << "Ek (GeV)"
1843  << std::endl << std::endl;
1844 
1845  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
1846 
1847  // Loop over the particle stack for this event
1848  for(int i = 0; i < mc.NParticles(); ++i){
1849  simb::MCParticle part(mc.GetParticle(i));
1850  std::string name;
1851  if (part.PdgCode() == 18040)
1852  name = "Ar40 18040";
1853  else if (part.PdgCode() != -99999 )
1854  {
1855  name = databasePDG->GetParticle(part.PdgCode())->GetName();
1856  }
1857 
1858  int code = part.StatusCode();
1859  std::string status = ParticleStatus(code);
1860  double mass = part.Mass();
1861  double energy = part.E();
1862  double Ek = (energy-mass); // Kinetic Energy (GeV)
1863 
1864  std::cout << std::setiosflags(std::ios::left) << std::setw(20) << name
1865  << std::setiosflags(std::ios::left) << std::setw(32) <<status
1866  << std::setw(18)<< energy
1867  << std::setw(18)<< mass
1868  << std::setw(18)<< Ek <<std::endl;
1869  }
1870 
1871  if(mc.GetNeutrino().CCNC() == simb::kCC){
1872 
1873  ///look for the outgoing lepton in the particle stack
1874  ///just interested in the first one
1875  for(int i = 0; i < mc.NParticles(); ++i){
1876  simb::MCParticle part(mc.GetParticle(i));
1877  if(std::abs(part.PdgCode()) == 11){
1878  fEMomentum->Fill(part.P());
1879  fEDCosX->Fill(part.Px()/part.P());
1880  fEDCosY->Fill(part.Py()/part.P());
1881  fEDCosZ->Fill(part.Pz()/part.P());
1882  fECons->Fill(nu.E() - part.E());
1883  break;
1884  }
1885  else if(std::abs(part.PdgCode()) == 13){
1886  fMuMomentum->Fill(part.P());
1887  fMuDCosX->Fill(part.Px()/part.P());
1888  fMuDCosY->Fill(part.Py()/part.P());
1889  fMuDCosZ->Fill(part.Pz()/part.P());
1890  fECons->Fill(nu.E() - part.E());
1891  break;
1892  }
1893  }// end loop over particles
1894  }//end if CC interaction
1895 
1896  // fill fDyn
1897  double bin(0.0);
1898  double binNew(0.0);
1899  double binNewThresh(0.0);
1900  if (NuWroTTree->flag.qel && NuWroTTree->flag.cc) bin = 1.;
1901  else if (NuWroTTree->flag.qel && NuWroTTree->flag.nc) bin = 2.;
1902  else if (NuWroTTree->flag.dis && NuWroTTree->flag.cc) bin = 10.;
1903  else if (NuWroTTree->flag.dis && NuWroTTree->flag.nc) bin = 11.;
1904  else if (NuWroTTree->flag.coh && NuWroTTree->flag.cc) bin = 12.;
1905  else if (NuWroTTree->flag.coh && NuWroTTree->flag.nc) bin = 13.;
1906  else if (NuWroTTree->flag.res)
1907  {
1908  unsigned int ii(0);
1909  unsigned int p(0), n(0), pip(0), pim(0), pi0(0);
1910  while(ii<NuWroTTree->post.size())
1911  {
1912  if (NuWroTTree->out[ii].pdg==211) pip++;
1913  if (NuWroTTree->out[ii].pdg==-211) pim++;
1914  if (NuWroTTree->out[ii].pdg==111) pi0++;
1915  if (NuWroTTree->out[ii].pdg==2112) n++;
1916  if (NuWroTTree->out[ii].pdg==2212) p++;
1917  ii++;
1918  }
1919  if (NuWroTTree->flag.cc && pip && p) bin = 3.;
1920  else if (NuWroTTree->flag.cc && pi0 && p) bin = 4.;
1921  else if (NuWroTTree->flag.cc && pip && n) bin = 5.;
1922  else if (NuWroTTree->flag.nc && pi0 && p) bin = 6.;
1923  else if (NuWroTTree->flag.nc && pip && n) bin = 7.;
1924  else if (NuWroTTree->flag.nc && pi0 && n) bin = 8.;
1925  else if (NuWroTTree->flag.nc && pim && p) bin = 9.;
1926  else
1927  std::cout << "NuWroGen: bin=0 events, cc?" << NuWroTTree->flag.cc << " nc?: " << NuWroTTree->flag.nc << " Num protons: " << p << " Num pi+s: " << pip << " Num pi-s: " << pim << " Num pi0s: " << pi0 << " Num neutrons: " << n << std::endl;
1928  }
1929 
1930  unsigned int ii(0);
1931  unsigned int p(0), n(0), pip(0), pim(0), pi0(0), pThresh(0.);
1932  while(ii<NuWroTTree->post.size())
1933  {
1934  if (NuWroTTree->out[ii].pdg==211) pip++;
1935  if (NuWroTTree->out[ii].pdg==-211) pim++;
1936  if (NuWroTTree->out[ii].pdg==111) pi0++;
1937  if (NuWroTTree->out[ii].pdg==2112) n++;
1938  if (NuWroTTree->out[ii].pdg==2212) p++;
1939  if (NuWroTTree->out[ii].pdg==2212 &&
1940  (NuWroTTree->out[ii].t/1000.-0.939)>0.050) pThresh++;
1941  ii++;
1942  }
1943 
1944  if (NuWroTTree->flag.cc && p==0 && (pip+pim)==0) binNew = 1;
1945  else if (NuWroTTree->flag.cc && p==1 && (pip+pim)==0) binNew = 2;
1946  else if (NuWroTTree->flag.cc && p>=2 && (pip+pim)==0) binNew = 3;
1947  else if (NuWroTTree->flag.cc && p==0 && (pip+pim)==1) binNew = 4;
1948  else if (NuWroTTree->flag.cc && p==1 && (pip+pim)==1) binNew = 5;
1949  else if (NuWroTTree->flag.cc && p>=2 && (pip+pim)==1) binNew = 6;
1950  else if (NuWroTTree->flag.cc && p==0 && (pip+pim)>=2) binNew = 7;
1951  else if (NuWroTTree->flag.cc && p>=1 && (pip+pim)>=2) binNew = 8;
1952  else if (NuWroTTree->flag.nc) binNew = 9;
1953  else binNew = 10;
1954 
1955  if (NuWroTTree->flag.cc && pThresh==0 && (pip+pim)==0) binNewThresh = 1;
1956  else if (NuWroTTree->flag.cc && pThresh==1 && (pip+pim)==0) binNewThresh = 2;
1957  else if (NuWroTTree->flag.cc && pThresh>=2 && (pip+pim)==0) binNewThresh = 3;
1958  else if (NuWroTTree->flag.cc && pThresh==0 && (pip+pim)==1) binNewThresh = 4;
1959  else if (NuWroTTree->flag.cc && pThresh==1 && (pip+pim)==1) binNewThresh = 5;
1960  else if (NuWroTTree->flag.cc && pThresh>=2 && (pip+pim)==1) binNewThresh = 6;
1961  else if (NuWroTTree->flag.cc && pThresh==0 && (pip+pim)>=2) binNewThresh = 7;
1962  else if (NuWroTTree->flag.cc && pThresh>=1 && (pip+pim)>=2) binNewThresh = 8;
1963  else if (NuWroTTree->flag.nc ) binNewThresh = 9;
1964  else binNewThresh = 10;
1965 
1966  // Needs to be normalized to 70 tonnes, 6e20pot
1967  double N_ArAtoms(70.*1000*1000/40*6.022e23);
1968  // arXiv:pdf/0806.1449v2.pdf adjusted to uBooNE distance
1969  double FluxNorm(5.19e-10 * (540./460.)*(540./460.));
1970  double nucleiPerAtom((double)(NuWroTTree->par.nucleus_p+NuWroTTree->par.nucleus_n));
1971  double NumEvtsRunThisJob(10000.);
1972  // The weight coming out of NuWro is proportional to the xsection for that process. Instead use the flux-wtd avg's for each dyn mechanism as given in the ouutput.root.txt file.
1973  // double wt = NuWroTTree->weight * N_ArAtoms * nucleiPerAtom * 6.e20 * FluxNorm / NumEvtsRunThisJob;
1974 
1975  // double wt = fxsecFluxWtd.at(NuWroTTree->dyn) * N_ArAtoms * nucleiPerAtom * 6.e20 * FluxNorm / NumEvtsRunThisJob;
1976  double wt = fxsecTotal * N_ArAtoms * nucleiPerAtom * 6.e20 * FluxNorm / NumEvtsRunThisJob;
1977  fDyn->Fill(bin-0.5,wt);
1978  fDynNew->Fill(binNew-0.5,wt);
1979  fDynNewThresh->Fill(binNewThresh-0.5,wt);
1980  f2DynNew->Fill(bin-0.5,binNew-0.5,wt);
1981  f2DynNewThresh->Fill(bin-0.5,binNewThresh-0.5,wt);
1982 
1983  fWeight->Fill(wt);
1984  fWeightNW->Fill(NuWroTTree->weight);
1985  return;
1986  }
name
static QCString name
Definition: declinfo.cpp:673
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
evgen::NuWroGen::fEDCosZ
TH1F * fEDCosZ
direction cosine of outgoing e in z
Definition: NuWroGen_module.cc:118
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
simb::MCNeutrino::CCNC
int CCNC() const
Definition: MCNeutrino.h:148
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
generate_CCQE_events.energy
string energy
Definition: generate_CCQE_events.py:15
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
evgen::NuWroGen::fVertexYZ
TH2F * fVertexYZ
vertex location in yz
Definition: NuWroGen_module.cc:104
string
std::string string
Definition: nybbler.cc:12
evgen::NuWroGen::fCCMode
TH1F * fCCMode
CC interaction mode.
Definition: NuWroGen_module.cc:120
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
evgen::NuWroGen::fVertexXY
TH2F * fVertexXY
vertex location in xy
Definition: NuWroGen_module.cc:102
evgen::NuWroGen::ReactionChannel
std::string ReactionChannel(int ccnc, int mode)
Definition: NuWroGen_module.cc:1724
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:230
evgen::NuWroGen::fWeightNW
TH1F * fWeightNW
NuWro Wt.
Definition: NuWroGen_module.cc:124
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
evgen::NuWroGen::ParticleStatus
std::string ParticleStatus(int StatusCode)
Definition: NuWroGen_module.cc:1700
evgen::NuWroGen::fDCosX
TH1F * fDCosX
direction cosine in x
Definition: NuWroGen_module.cc:106
evgen::NuWroGen::fGenerated
TH1F * fGenerated[6]
Spectra as generated.
Definition: NuWroGen_module.cc:96
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
evgen::NuWroGen::fVertexXZ
TH2F * fVertexXZ
vertex location in xz
Definition: NuWroGen_module.cc:103
evgen::NuWroGen::fNCMode
TH1F * fNCMode
CC interaction mode.
Definition: NuWroGen_module.cc:121
evgen::NuWroGen::fEDCosY
TH1F * fEDCosY
direction cosine of outgoing e in y
Definition: NuWroGen_module.cc:117
e
const double e
Definition: gUpMuFluxGen.cxx:165
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
create_imaplist.status
int status
Definition: create_imaplist.py:5
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:234
evgen::NuWroGen::fWeight
TH1F * fWeight
NuWro Wt.
Definition: NuWroGen_module.cc:123
test.p
p
Definition: test.py:223
code
CodeOutputInterface * code
Definition: tclscanner.cpp:1114
evgen::NuWroGen::fVertexY
TH1F * fVertexY
vertex location of generated events in y
Definition: NuWroGen_module.cc:99
evgen::NuWroGen::fVertexZ
TH1F * fVertexZ
vertex location of generated events in z
Definition: NuWroGen_module.cc:100
setw
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
simb::MCTruth::GetParticle
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
evgen::NuWroGen::fDCosY
TH1F * fDCosY
direction cosine in y
Definition: NuWroGen_module.cc:107
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
evgen::NuWroGen::fECons
TH1F * fECons
histogram to determine if energy is conserved in the event
Definition: NuWroGen_module.cc:132
evgen::NuWroGen::fEDCosX
TH1F * fEDCosX
direction cosine of outgoing e in x
Definition: NuWroGen_module.cc:116
evgen::NuWroGen::fDCosZ
TH1F * fDCosZ
direction cosine in z
Definition: NuWroGen_module.cc:108
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
bin
QTextStream & bin(QTextStream &s)
Definition: qtextstream.cpp:2006
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
evgen::NuWroGen::fMuDCosX
TH1F * fMuDCosX
direction cosine of outgoing mu in x
Definition: NuWroGen_module.cc:111
evgen::NuWroGen::fDyn
TH1F * fDyn
mode in detail
Definition: NuWroGen_module.cc:122
evgen::NuWroGen::fMuDCosY
TH1F * fMuDCosY
direction cosine of outgoing mu in y
Definition: NuWroGen_module.cc:112
evgen::NuWroGen::fEMomentum
TH1F * fEMomentum
momentum of outgoing electrons
Definition: NuWroGen_module.cc:115
simb::MCNeutrino
Event generator information.
Definition: MCNeutrino.h:18
simb::MCNeutrino::Mode
int Mode() const
Definition: MCNeutrino.h:149
evgen::NuWroGen::fMuMomentum
TH1F * fMuMomentum
momentum of outgoing muons
Definition: NuWroGen_module.cc:110
evgen::NuWroGen::fVertexX
TH1F * fVertexX
vertex location of generated events in x
Definition: NuWroGen_module.cc:98
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:222
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
evgen::NuWroGen::fMuDCosZ
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
Definition: NuWroGen_module.cc:113
std::string evgen::NuWroGen::ParticleStatus ( int  StatusCode)
private

Definition at line 1700 of file NuWroGen_module.cc.

1701  {
1702  int code = StatusCode;
1703  std::string ParticleStatusName;
1704 
1705  switch(code)
1706  {
1707  case 0:
1708  ParticleStatusName = "kIStInitialState";
1709  break;
1710  case 1:
1711  ParticleStatusName = "kIStFinalState";
1712  break;
1713  case 11:
1714  ParticleStatusName = "kIStNucleonTarget";
1715  break;
1716  default:
1717  ParticleStatusName = "Status Unknown";
1718  }
1719  return ParticleStatusName;
1720  }
string
std::string string
Definition: nybbler.cc:12
code
CodeOutputInterface * code
Definition: tclscanner.cpp:1114
sim::ParticleStatusName
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
Definition: MCDumperUtils.cxx:126
void evgen::NuWroGen::produce ( art::Event &  evt)
virtual

Implements art::EDProducer.

Definition at line 1509 of file NuWroGen_module.cc.

1510  {
1511 
1512  std::cout << std::endl;
1513  std::cout<<"------------------------------------------------------------------------------"<<std::endl;
1514 // std::cout << "run : " << evt.Header().Run() << std::endl;
1515 // std::cout << "subrun : " << evt.Header().Subrun() << std::endl;
1516 // std::cout << "event : " << evt.Header().Event() << std::endl;
1517 // std::cout << "event : " << evt.id().event() << std::endl;
1518 
1519  std::string name, k, dollar;
1520  int partnumber = 0;
1521 
1522  int trackid = -1; // set track id to -i as these are all primary particles and have id <= 0
1523  std::string primary("primary");
1524  int FirstMother = -1;
1525 
1526  int Status = -9999;
1527 
1528 
1529  int ccnc = -9999;
1530  int mode = -9999;
1531  int targetnucleusPdg = -9999;
1532  int hitquarkPdg = -9999;
1533 
1534  TLorentzVector Neutrino;
1535  TLorentzVector Lepton;
1536  TLorentzVector Target;
1537  TLorentzVector q;
1538  TLorentzVector Hadron4mom;
1539  double Q2 = -9999;
1540 
1541  int Tpdg = 0; // for target
1542  double Tmass = 0;
1543  int Tstatus = 11;
1544  double Tcosx, Tcosy, Tcosz, Tenergy;
1545  TLorentzVector Tpos;
1546 
1547 
1548  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
1549  simb::MCTruth truth;
1550 
1551  /*
1552  if (countFile>=ch->GetEntries())
1553  {
1554  mf::LogInfo("NuWro (parser): You're on row ") << countFile << " of " << ch->GetEntries()<< ". Moving on ..." <<std::endl;
1555  return;
1556  }
1557  */
1558  ch->GetEntry(countFile++);
1559 
1560  //get the nuwro channel number
1561 
1562  //set the interaction type; CC or NC
1563  if (NuWroTTree->flag.cc) ccnc = simb::kCC;
1564  else if (NuWroTTree->flag.nc) ccnc = simb::kNC;
1565 
1566  //set the interaction mode; QE, Res, DIS, Coh, kNuElectronElastic, kInverseMuDecay
1567  if ( NuWroTTree->flag.qel )
1568  mode = simb::kQE;
1569  else if ( NuWroTTree->flag.res )
1570  mode = simb::kRes;
1571  else if ( NuWroTTree->flag.dis )
1572  mode = simb::kDIS;
1573  else if ( NuWroTTree->flag.coh )
1574  mode = simb::kCoh;
1575  if(partnumber == -1)
1576  Status = 0;
1577  else
1578  Status = 1;
1579 
1580 
1581  art::ServiceHandle<geo::Geometry const> geo;
1582  double X0 = NuWroTTree->par.geo_o[0] + geo->DetHalfWidth();
1583  double Y0 = NuWroTTree->par.geo_o[1];
1584  double Z0 = NuWroTTree->par.geo_o[2] + 0.25*geo->DetLength();
1585  TLorentzVector pos(X0, Y0, Z0, 0);
1586  Tpos = pos; // for target
1587 
1588 
1589  simb::MCParticle mcpartNu(trackid,
1590  NuWroTTree->in[0].pdg,
1591  primary,
1592  FirstMother,
1593  NuWroTTree->in[0].m()/1000.,
1594  Status
1595  );
1596  TLorentzVector mom(NuWroTTree->in[0].x/1000.,
1597  NuWroTTree->in[0].y/1000.,
1598  NuWroTTree->in[0].z/1000.,
1599  NuWroTTree->in[0].t/1000.);
1600  mcpartNu.AddTrajectoryPoint(pos,mom);
1601  truth.Add(mcpartNu);
1602 
1603 
1604  unsigned int ii(0);
1605  while(ii<NuWroTTree->post.size())
1606  {
1607  // loop over particles in an event
1608 
1609  simb::MCParticle mcpart(trackid,
1610  NuWroTTree->post[ii].pdg,
1611  primary,
1612  FirstMother,
1613  NuWroTTree->post[ii].m()/1000.,
1614  Status
1615  );
1616  TLorentzVector mom(NuWroTTree->post[ii].x/1000.,
1617  NuWroTTree->post[ii].y/1000.,
1618  NuWroTTree->post[ii].z/1000.,
1619  NuWroTTree->post[ii].t/1000.);
1620  mcpart.AddTrajectoryPoint(pos,mom);
1621  truth.Add(mcpart);
1622 
1623 
1624  ii++;
1625  }// loop over particles in an event
1626  // Incoming Neutrino is 0th element of in. Outgoing lepton is 0th of out.
1627  Neutrino.SetPxPyPzE(NuWroTTree->in[0].x/1000., NuWroTTree->in[0].y/1000., NuWroTTree->in[0].z/1000., NuWroTTree->in[0].t/1000. );
1628  Lepton.SetPxPyPzE(NuWroTTree->out[0].x/1000., NuWroTTree->out[0].y/1000., NuWroTTree->out[0].z/1000., NuWroTTree->out[0].t/1000. );
1629 
1630  /////////////////////////////////
1631 
1632  Tmass = NuWroTTree->par.nucleus_p + NuWroTTree->par.nucleus_n; // GeV
1633 
1634 
1635  Tenergy = NuWroTTree->in.back().t;
1636  Tcosx = NuWroTTree->in.back().x;
1637  Tcosy = NuWroTTree->in.back().y;
1638  Tcosz = NuWroTTree->in.back().z;
1639  Tmass = std::sqrt(std::abs(Tenergy*Tenergy - Tcosx*Tcosx
1640  - Tcosy*Tcosy - Tcosz*Tcosz))/1000.;
1641 
1642  Tenergy = Tmass-0.1; // force this negative, cuz kinetic energy>eps
1643  // seems to make G4 hang.
1644  Tcosx =0.; Tcosy = 0.; Tcosz = 0.;
1645 
1646 
1647  simb::MCParticle mcpartT(trackid,
1648  Tpdg,
1649  primary,
1650  FirstMother,
1651  Tmass,
1652  Tstatus
1653  );
1654 
1655 
1656  // Target = Hadron4mom - (Neutrino - Lepton); // commenting this out as target momentum no more is calculated by 4-momentum conservation
1657 
1658  TLorentzVector Tmom;
1659  Tmom.SetPxPyPzE(Tcosx/1000., Tcosy/1000., Tcosz/1000., Tenergy); // this makes literally |P| = 0 or 1 GeV/c for target, this affects only the Kinematic variables; X and Y
1660  // target |p| = 0 GeV/c if interaction is
1661  // DIS, Coh, nu-e Elastic Scattering, nu-e inverse mu decay (this comes from Nuwro),
1662  // else target |P| = 1 GeV/c
1663  Target = Tmom;
1664 
1665 
1666  mcpartT.AddTrajectoryPoint(Tpos,Tmom);
1667  // for now, do(n't) target onto stack. EC, 20-Jul-2012.
1668  truth.Add(mcpartT);
1669 
1670  q = Neutrino - Lepton;
1671  Q2 = -(q * q);
1672  // double W2 = Hadron4mom * Hadron4mom;
1673  // double InvariantMass = std::sqrt(W2);
1674 
1675  double x = Q2/((2*Target*q));
1676  double y = (Target*q)/(Neutrino*Target);
1677 
1678  truth.SetOrigin(simb::kBeamNeutrino);
1679  int channel(-999);
1680  targetnucleusPdg = NuWroTTree->in[NuWroTTree->in.size()-1].pdg;
1681  truth.SetNeutrino(ccnc, mode, channel,
1682  targetnucleusPdg,
1683  Tpdg,
1684  hitquarkPdg,
1685  //InvariantMass, x, y, Q2
1686  0, x, y, Q2
1687  );
1688 
1689  std::cout << truth.GetNeutrino() << std::endl;
1690 
1691  FillHistograms(truth);
1692 
1693  truthcol->push_back(truth);
1694  evt.put(std::move(truthcol));
1695 
1696  return;
1697  }
name
static QCString name
Definition: declinfo.cpp:673
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
geo::GeometryCore::DetHalfWidth
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
Definition: GeometryCore.cxx:752
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
simb::MCTruth::SetOrigin
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
string
std::string string
Definition: nybbler.cc:12
simb::MCTruth::SetNeutrino
void SetNeutrino(int CCNC, int mode, int interactionType, int target, int nucleon, int quark, double w, double x, double y, double qsqr)
Definition: MCTruth.cxx:30
channel
uint8_t channel
Definition: CRTFragment.hh:201
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
wirecell.gen.depos.move
def move(depos, offset)
Definition: depos.py:107
geo::GeometryCore::DetLength
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
Definition: GeometryCore.cxx:764
art::DataViewImpl::put
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
simb::MCTruth::Add
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
evgen::NuWroGen::FillHistograms
void FillHistograms(const simb::MCTruth &mc)
Definition: NuWroGen_module.cc:1752
train.x
list x
Definition: train.py:276
evgen::NuWroGen::countFile
unsigned int countFile
Definition: NuWroGen_module.cc:136
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
simb::kBeamNeutrino
Beam neutrinos.
Definition: MCTruth.h:23
std::string evgen::NuWroGen::ReactionChannel ( int  ccnc,
int  mode 
)
private

Definition at line 1724 of file NuWroGen_module.cc.

1725  {
1726  std::string ReactionChannelName=" ";
1727 
1728  if(ccnc==0)
1729  ReactionChannelName = "kCC";
1730  else if(ccnc==1)
1731  ReactionChannelName = "kNC";
1732  else std::cout<<"Current mode unknown!! "<<std::endl;
1733 
1734  if(mode==0)
1735  ReactionChannelName += "_kQE";
1736  else if(mode==1)
1737  ReactionChannelName += "_kRes";
1738  else if(mode==2)
1739  ReactionChannelName += "_kDIS";
1740  else if(mode==3)
1741  ReactionChannelName += "_kCoh";
1742  else if(mode==4)
1743  ReactionChannelName += "_kNuElectronElastic";
1744  else if(mode==5)
1745  ReactionChannelName += "_kInverseMuDecay";
1746  else std::cout<<"interaction mode unknown!! "<<std::endl;
1747 
1748  return ReactionChannelName;
1749  }
string
std::string string
Definition: nybbler.cc:12
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
void evgen::NuWroGen::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 1274 of file NuWroGen_module.cc.

1275  {
1276  fFileName = p.get< std::string >("NuWroFile","output.root");
1277  fTreeName = p.get< std::string >("TreeName","treeout");
1278 
1279  return;
1280  }
string
std::string string
Definition: nybbler.cc:12
evgen::NuWroGen::fTreeName
std::string fTreeName
Definition: NuWroGen_module.cc:85
test.p
p
Definition: test.py:223
evgen::NuWroGen::fFileName
std::string fFileName
Definition: NuWroGen_module.cc:84

Member Data Documentation

TChain* evgen::NuWroGen::ch
private

Definition at line 138 of file NuWroGen_module.cc.

unsigned int evgen::NuWroGen::countFile
private

Definition at line 136 of file NuWroGen_module.cc.

TH2F* evgen::NuWroGen::f2DynNew
private

Definition at line 127 of file NuWroGen_module.cc.

TH2F* evgen::NuWroGen::f2DynNewThresh
private

Definition at line 128 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fCCMode
private

CC interaction mode.

Definition at line 120 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fDCosX
private

direction cosine in x

Definition at line 106 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fDCosY
private

direction cosine in y

Definition at line 107 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fDCosZ
private

direction cosine in z

Definition at line 108 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fDyn
private

mode in detail

Definition at line 122 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fDynNew
private

Definition at line 125 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fDynNewThresh
private

Definition at line 126 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fECons
private

histogram to determine if energy is conserved in the event

Definition at line 132 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fEDCosX
private

direction cosine of outgoing e in x

Definition at line 116 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fEDCosY
private

direction cosine of outgoing e in y

Definition at line 117 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fEDCosZ
private

direction cosine of outgoing e in z

Definition at line 118 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fEMomentum
private

momentum of outgoing electrons

Definition at line 115 of file NuWroGen_module.cc.

std::ifstream* evgen::NuWroGen::fEventFile
private

Definition at line 86 of file NuWroGen_module.cc.

int evgen::NuWroGen::fEventNumberOffset
private

Definition at line 91 of file NuWroGen_module.cc.

std::string evgen::NuWroGen::fFileName
private

Definition at line 84 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fGenerated[6]
private

Spectra as generated.

Definition at line 96 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fMuDCosX
private

direction cosine of outgoing mu in x

Definition at line 111 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fMuDCosY
private

direction cosine of outgoing mu in y

Definition at line 112 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fMuDCosZ
private

direction cosine of outgoing mu in z

Definition at line 113 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fMuMomentum
private

momentum of outgoing muons

Definition at line 110 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fNCMode
private

CC interaction mode.

Definition at line 121 of file NuWroGen_module.cc.

std::string evgen::NuWroGen::fNuWroModuleLabel
private

keep track of how long it takes to run the job

Definition at line 90 of file NuWroGen_module.cc.

TStopwatch evgen::NuWroGen::fStopwatch
private

Definition at line 88 of file NuWroGen_module.cc.

std::string evgen::NuWroGen::fTreeName
private

Definition at line 85 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fVertexX
private

vertex location of generated events in x

Definition at line 98 of file NuWroGen_module.cc.

TH2F* evgen::NuWroGen::fVertexXY
private

vertex location in xy

Definition at line 102 of file NuWroGen_module.cc.

TH2F* evgen::NuWroGen::fVertexXZ
private

vertex location in xz

Definition at line 103 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fVertexY
private

vertex location of generated events in y

Definition at line 99 of file NuWroGen_module.cc.

TH2F* evgen::NuWroGen::fVertexYZ
private

vertex location in yz

Definition at line 104 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fVertexZ
private

vertex location of generated events in z

Definition at line 100 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fWeight
private

NuWro Wt.

Definition at line 123 of file NuWroGen_module.cc.

TH1F* evgen::NuWroGen::fWeightNW
private

NuWro Wt.

Definition at line 124 of file NuWroGen_module.cc.

std::vector<double> evgen::NuWroGen::fxsecFluxWtd
private

Definition at line 93 of file NuWroGen_module.cc.

double evgen::NuWroGen::fxsecTotal
private

Definition at line 94 of file NuWroGen_module.cc.

event* evgen::NuWroGen::NuWroTTree
private

Definition at line 137 of file NuWroGen_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.11