Skip to content

Commit 2db478b

Browse files
kjplowskjplows
committed
Flatten the GENIE event record in flat CAF
1 parent ccdf1cb commit 2db478b

1 file changed

Lines changed: 90 additions & 2 deletions

File tree

‎sbncode/CAFMaker/CAFMaker_module.cc‎

Lines changed: 90 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -119,6 +119,9 @@
119119
// GENIE
120120
#include "Framework/EventGen/EventRecord.h"
121121
#include "Framework/Ntuple/NtpMCEventRecord.h"
122+
#include "Framework/Conventions/Units.h"
123+
#include "Framework/GHEP/GHepParticle.h"
124+
122125
#include "nugen/EventGeneratorBase/GENIE/GENIE2ART.h"
123126

124127
#include "canvas/Persistency/Provenance/ProcessConfiguration.h"
@@ -246,6 +249,27 @@ class CAFMaker : public art::EDProducer {
246249
bool fSaveGENIEEventRecord;
247250
unsigned int fGenieEventCounter;
248251

252+
//TBits * fGenieEvtRec_brEvtFlags = 0; ////< Generator-specific event flags
253+
//TObjString * fGenieEvtRec_brEvtCode = 0; ////< Generator-specific string with 'event code'
254+
int fGenieEvtRec_brEvtNum = 0; ////< Event number
255+
double fGenieEvtRec_brEvtXSec = 0.0; ////< Cross section for selected event (1e-38 cm2)
256+
double fGenieEvtRec_brEvtDXSec = 0.0; ////< Cross section for selected event kinematics (1e-38 cm2 / {K^n})
257+
unsigned int fGenieEvtRec_brEvtKPS = 0; ////< Kinematic phase space variables. See $GENIE/src/Framework/Conventions/KinePhaseSpace.h -> KinePhaseSpace_t
258+
double fGenieEvtRec_brEvtWght = 0.0; ////< Weight for that event
259+
double fGenieEvtRec_brEvtProb = 0.0; ////< Probability for that event (given cross section, path lengths, etc)
260+
double fGenieEvtRec_brEvtVtx[4] = {0.0}; ////< Event vertex position in detector coord syst (SI)
261+
int fGenieEvtRec_brStdHepN = 0; ////< Number of particles in particle array
262+
int fGenieEvtRec_brStdHepPdg [250] = {0}; ////< Pdg codes (& generator specific codes for pseudoparticles)
263+
int fGenieEvtRec_brStdHepStatus[250] = {0}; ////< Generator-specific status code
264+
int fGenieEvtRec_brStdHepRescat[250] = {0}; ////< Hadron transport model-specific rescattering code
265+
double fGenieEvtRec_brStdHepX4 [250][4] = {{0.0}}; ////< 4-x (x, y, z, t) of particle in hit nucleus frame (fm)
266+
double fGenieEvtRec_brStdHepP4 [250][4] = {{0.0}}; ////< 4-p (px,py,pz,E) of particle in LAB frame (GeV)
267+
double fGenieEvtRec_brStdHepPolz [250][3] = {{0.0}}; ////< Polarization vector
268+
int fGenieEvtRec_brStdHepFd [250] = {0}; ////< First daughter
269+
int fGenieEvtRec_brStdHepLd [250] = {0}; ////< Last daughter
270+
int fGenieEvtRec_brStdHepFm [250] = {0}; ////< First mother
271+
int fGenieEvtRec_brStdHepLm [250] = {0}; ////< Last mother
272+
249273
flat::Flat<caf::StandardRecord>* fFlatRecord = 0;
250274
flat::Flat<caf::StandardRecord>* fFlatRecordb = 0;
251275
flat::Flat<caf::StandardRecord>* fFlatRecordp = 0;
@@ -1050,7 +1074,32 @@ void CAFMaker::InitializeOutfiles()
10501074

10511075
if (fSaveGENIEEventRecord){
10521076
fFlatGenieTree = new TTree( "GenieEvtRecTree", "GenieEvtRecTree" );
1053-
fFlatGenieTree->Branch("GenieEvtRec", &fFlatGenieEvtRec);
1077+
//fFlatGenieTree->Branch("GenieEvtRec", &fFlatGenieEvtRec);
1078+
1079+
// Flatten the event record.
1080+
// Follows the convention of the GENIE numi_rootracker (see $GENIE/src/Apps/gNtpConv.cxx).
1081+
//fFlatGenieTree->Branch("GenieEvtRec.EvtFlags", "TBits", fGenieEvtRec_brEvtFlags, 32000, 1);
1082+
//fFlatGenieTree->Branch("GenieEvtRec.EvtCode", "TObjString", fGenieEvtRec_brEvtCode, 32000, 1);
1083+
1084+
fFlatGenieTree->Branch("GenieEvtRec.EvtNum", &fGenieEvtRec_brEvtNum, "GenieEvtRec.EvtNum/I" );
1085+
fFlatGenieTree->Branch("GenieEvtRec.EvtXSec", &fGenieEvtRec_brEvtXSec, "GenieEvtRec.EvtXSec/D" );
1086+
fFlatGenieTree->Branch("GenieEvtRec.EvtDXSec", &fGenieEvtRec_brEvtDXSec, "GenieEvtRec.EvtDXSec/D" );
1087+
fFlatGenieTree->Branch("GenieEvtRec.EvtKPS", &fGenieEvtRec_brEvtKPS, "GenieEvtRec.EvtKPS/i" );
1088+
fFlatGenieTree->Branch("GenieEvtRec.EvtWght", &fGenieEvtRec_brEvtWght, "GenieEvtRec.EvtWght/D" );
1089+
fFlatGenieTree->Branch("GenieEvtRec.EvtProb", &fGenieEvtRec_brEvtProb, "GenieEvtRec.EvtProb/D" );
1090+
fFlatGenieTree->Branch("GenieEvtRec.EvtVtx", fGenieEvtRec_brEvtVtx, "GenieEvtRec.EvtVtx[4]/D" );
1091+
fFlatGenieTree->Branch("GenieEvtRec.StdHepN", &fGenieEvtRec_brStdHepN, "GenieEvtRec.StdHepN/I" );
1092+
fFlatGenieTree->Branch("GenieEvtRec.StdHepPdg", fGenieEvtRec_brStdHepPdg, "GenieEvtRec.StdHepPdg[GenieEvtRec.StdHepN]/I" );
1093+
fFlatGenieTree->Branch("GenieEvtRec.StdHepStatus", fGenieEvtRec_brStdHepStatus, "GenieEvtRec.StdHepStatus[GenieEvtRec.StdHepN]/I" );
1094+
fFlatGenieTree->Branch("GenieEvtRec.StdHepRescat", fGenieEvtRec_brStdHepRescat, "GenieEvtRec.StdHepRescat[GenieEvtRec.StdHepN]/I" );
1095+
fFlatGenieTree->Branch("GenieEvtRec.StdHepX4", fGenieEvtRec_brStdHepX4, "GenieEvtRec.StdHepX4[GenieEvtRec.StdHepN][4]/D" );
1096+
fFlatGenieTree->Branch("GenieEvtRec.StdHepP4", fGenieEvtRec_brStdHepP4, "GenieEvtRec.StdHepP4[GenieEvtRec.StdHepN][4]/D" );
1097+
fFlatGenieTree->Branch("GenieEvtRec.StdHepPolz", fGenieEvtRec_brStdHepPolz, "GenieEvtRec.StdHepPolz[GenieEvtRec.StdHepN][3]/D" );
1098+
fFlatGenieTree->Branch("GenieEvtRec.StdHepFd", fGenieEvtRec_brStdHepFd, "GenieEvtRec.StdHepFd[GenieEvtRec.StdHepN]/I" );
1099+
fFlatGenieTree->Branch("GenieEvtRec.StdHepLd", fGenieEvtRec_brStdHepLd, "GenieEvtRec.StdHepLd[GenieEvtRec.StdHepN]/I" );
1100+
fFlatGenieTree->Branch("GenieEvtRec.StdHepFm", fGenieEvtRec_brStdHepFm, "GenieEvtRec.StdHepFm[GenieEvtRec.StdHepN]/I" );
1101+
fFlatGenieTree->Branch("GenieEvtRec.StdHepLm", fGenieEvtRec_brStdHepLm, "GenieEvtRec.StdHepLm[GenieEvtRec.StdHepN]/I" );
1102+
10541103
fFlatGenieTree->Branch("GENIEEntry", &fGenieEventCounter, "GENIE/i");
10551104
fFlatGenieTree->Branch("SourceFileHash", &fSourceFileHash, "SourceFileHash/i");
10561105
}
@@ -1394,7 +1443,46 @@ void CAFMaker::produce(art::Event& evt) noexcept {
13941443
fGenieTree->Fill();
13951444
}
13961445
if(fFlatGenieTree){
1397-
fFlatGenieEvtRec->Fill(fGenieEventCounter, genie_rec);
1446+
//fFlatGenieEvtRec->Fill(fGenieEventCounter, genie_rec);
1447+
fGenieEvtRec_brEvtNum = fGenieEventCounter;
1448+
fGenieEvtRec_brEvtXSec = genie_rec->XSec() * (1e+38/genie::units::cm2);
1449+
fGenieEvtRec_brEvtDXSec = genie_rec->DiffXSec() * (1e+38/genie::units::cm2);
1450+
fGenieEvtRec_brEvtKPS = genie_rec->DiffXSecVars();
1451+
fGenieEvtRec_brEvtWght = genie_rec->Weight();
1452+
fGenieEvtRec_brEvtProb = genie_rec->Probability();
1453+
fGenieEvtRec_brEvtVtx[0] = genie_rec->Vertex()->X();
1454+
fGenieEvtRec_brEvtVtx[1] = genie_rec->Vertex()->Y();
1455+
fGenieEvtRec_brEvtVtx[2] = genie_rec->Vertex()->Z();
1456+
fGenieEvtRec_brEvtVtx[3] = genie_rec->Vertex()->T();
1457+
1458+
int iparticle=0;
1459+
genie::GHepParticle * p = 0;
1460+
while( genie_rec->Particle(iparticle) != 0 ) {
1461+
p = genie_rec->Particle(iparticle);
1462+
fGenieEvtRec_brStdHepPdg[iparticle] = p->Pdg();
1463+
fGenieEvtRec_brStdHepStatus[iparticle] = (int) p->Status();
1464+
fGenieEvtRec_brStdHepRescat[iparticle] = p->RescatterCode();
1465+
fGenieEvtRec_brStdHepX4 [iparticle][0] = p->X4()->X();
1466+
fGenieEvtRec_brStdHepX4 [iparticle][1] = p->X4()->Y();
1467+
fGenieEvtRec_brStdHepX4 [iparticle][2] = p->X4()->Z();
1468+
fGenieEvtRec_brStdHepX4 [iparticle][3] = p->X4()->T();
1469+
fGenieEvtRec_brStdHepP4 [iparticle][0] = p->P4()->Px();
1470+
fGenieEvtRec_brStdHepP4 [iparticle][1] = p->P4()->Py();
1471+
fGenieEvtRec_brStdHepP4 [iparticle][2] = p->P4()->Pz();
1472+
fGenieEvtRec_brStdHepP4 [iparticle][3] = p->P4()->E();
1473+
if(p->PolzIsSet()) {
1474+
fGenieEvtRec_brStdHepPolz [iparticle][0] = TMath::Sin(p->PolzPolarAngle()) * TMath::Cos(p->PolzAzimuthAngle());
1475+
fGenieEvtRec_brStdHepPolz [iparticle][1] = TMath::Sin(p->PolzPolarAngle()) * TMath::Sin(p->PolzAzimuthAngle());
1476+
fGenieEvtRec_brStdHepPolz [iparticle][2] = TMath::Cos(p->PolzPolarAngle());
1477+
}
1478+
fGenieEvtRec_brStdHepFd [iparticle] = p->FirstDaughter();
1479+
fGenieEvtRec_brStdHepLd [iparticle] = p->LastDaughter();
1480+
fGenieEvtRec_brStdHepFm [iparticle] = p->FirstMother();
1481+
fGenieEvtRec_brStdHepLm [iparticle] = p->LastMother();
1482+
iparticle++;
1483+
}
1484+
fGenieEvtRec_brStdHepN = iparticle;
1485+
13981486
fFlatGenieTree->Fill();
13991487
}
14001488
}

0 commit comments

Comments
 (0)