SND@LHC Software
Loading...
Searching...
No Matches
sndlhcGSimpleNtpConverter.cxx
Go to the documentation of this file.
1#include <iostream>
2#include <fstream>
3#include <set>
4
5#include "TFile.h"
6#include "TTree.h"
7#include "TMath.h"
8#include "TRandom3.h"
9
10#include "ShipUnit.h"
11
12#include "Tools/Flux/GSimpleNtpFlux.h"
13
14/*
15 Converts neutrino rays generated by FLUKA into GSimpleNtpFlux format for GENIE
16*/
17
18int main(int argc, char **argv)
19{
20
21 if (argc != 4) {
22 std::cout << "Three arguments required: path to FLUKA file AND output file name AND number of pp collisions used "
23 "to generate FLUKA file."
24 << std::endl;
25 return -1;
26 }
27
28 std::string inFileName = std::string(argv[1]);
29 std::string outFileName = std::string(argv[2]);
30 double pp_collision_number = std::stod(argv[3]);
31
32 // Convert FLUKA to PDG particle IDs. (ONLY NEUTRINOS IMPLEMENTED!!!)
33 std::map<int, int> FLUKAtoPDG{{27, 14}, {28, -14}, {5, 12}, {6, -12}, {43, 16}, {44, -16}};
34
35 // To generate a random seed which is stored in the metadata.
36 TRandom3 *ran = new TRandom3();
37
38 bool verbose = false;
39
40 std::cout << "Start converting" << std::endl;
41
42 // Set up input file
43 std::ifstream in_file(inFileName);
44
45 // Input variables
46 /*
47 # Scoring particles entering Region No 2935
48 # Col 1: FLUKA run number
49 # Col 2: primary event number
50 # -- Particle information --
51 # Col 3: FLUKA particle type ID
52 # Col 4: Generation number
53 # Col 5: Kinetic energy (GeV)
54 # Col 6: Statistical weight
55 # -- Crossing at scoring plane --
56 # Col 7: x coord (cm)
57 # Col 8: y coord (cm)
58 # Col 9: x dir cosine
59 # Col 10: y dir cosine
60 # Col 11: z coord (cm)
61 # Col 12: Particle age since primary event (sec)
62 # Col 13: Last decay x cooord (cm)
63 # Col 14: Last decay y cooord (cm)
64 # Col 15: Last decay z cooord (cm)
65 # Col 16: Last decay ID
66 # Col 17: Kinetic energy of decay parent (GeV)
67 # Col 18: Last interaction x cooord (cm)
68 # Col 19: Last interaction y cooord (cm)
69 # Col 20: Last interaction z cooord (cm)
70 # Col 21: Last interaction ID
71 # Col 22: Kinetic energy of parent (GeV)
72 */
73
74 int FlukaRun, evtNumber, FlukaID, generationN, last_decay_ID, last_interaction_ID;
75 double Ekin, weight, x, y, z, x_cos, y_cos, age, Ekin_parent;
76 double last_decay_x, last_decay_y, last_decay_z, Ekin_decay;
77 double last_interaction_x, last_interaction_y, last_interaction_z;
78
79 // Scoring plane info
80 // Find a Minimum and Maximum of the x, y, z axis (FLUKA coordinates) cm
81 // Max z, min z
82 double min_z = 999;
83 double max_z = -999;
84 // Max x, min x
85 double min_x = 999;
86 double max_x = -999;
87 // Max y, min y
88 double min_y = 999;
89 double max_y = -999;
90
91 TFile *fOut = new TFile(outFileName.c_str(), "RECREATE");
92 TTree *tOut = new TTree("flux", "a simple flux n-tuple");
93 TTree *metaOut = new TTree("meta", "metadata for flux n-tuple");
94
95 genie::flux::GSimpleNtpEntry *gsimple_entry = new genie::flux::GSimpleNtpEntry;
96 tOut->Branch("entry", &gsimple_entry);
97 genie::flux::GSimpleNtpMeta *meta_entry = new genie::flux::GSimpleNtpMeta;
98 metaOut->Branch("meta", &meta_entry);
99
100 genie::flux::GSimpleNtpAux *aux_entry = new genie::flux::GSimpleNtpAux;
101 tOut->Branch("aux", &aux_entry);
102
103 UInt_t metakey = TString::Hash(outFileName.c_str(), strlen(outFileName.c_str()));
104 metakey &= 0x7FFFFFFF;
105
106 // Metadata accumulators:
107 std::set<int> pdglist;
108 double min_weight = 1e10;
109 double max_weight = -1e10;
110 double max_energy = 0;
111
112 unsigned long int counter = 0;
113
114 if (in_file.is_open()) {
115 string in_line;
116
117 while (!in_file.eof()) {
118 getline(in_file, in_line);
119
120 // Skip lines containing "#"
121 if (in_line.find("#") == in_line.npos) {
122 in_file >> FlukaRun >> evtNumber >> FlukaID >> generationN >> Ekin >> weight >> x >> y >> x_cos >> y_cos >>
123 z >> age >> last_decay_x >> last_decay_y >> last_decay_z >> last_decay_ID >> Ekin_decay >>
124 last_interaction_x >> last_interaction_y >> last_interaction_z >> last_interaction_ID >> Ekin_parent;
125
126 gsimple_entry->Reset();
127 aux_entry->Reset();
128
129 // FLUKA z = 0 in sndsw coordinates : -48000 cm
130 z -= 48000; // In cm for consistency with the FLUKA file.
131
132 min_z = std::min(min_z, z);
133 max_z = std::max(max_z, z);
134 min_x = std::min(min_x, x);
135 max_x = std::max(max_x, x);
136 min_y = std::min(min_y, y);
137 max_y = std::max(max_y, y);
138
139 if (verbose)
140 std::cout << "Got entry " << counter++ << std::endl;
141
142 if (verbose) {
143 std::cout << FlukaRun << "\n"
144 << evtNumber << "\n"
145 << FlukaID << "\n"
146 << generationN << "\n"
147 << Ekin << "\n"
148 << weight << "\n"
149 << x << "\n"
150 << y << "\n"
151 << x_cos << "\n"
152 << y_cos << "\n"
153 << z << "\n"
154 << age << "\n"
155 << last_decay_x << "\n"
156 << last_decay_y << "\n"
157 << last_decay_z << "\n"
158 << last_decay_ID << "\n"
159 << Ekin_decay << "\n"
160 << last_interaction_x << "\n"
161 << last_interaction_y << "\n"
162 << last_interaction_z << "\n"
163 << last_interaction_ID << "\n"
164 << Ekin_parent << "\n"
165 << "--------------------------------------------------" << std::endl;
166 }
167
168 gsimple_entry->metakey = metakey;
169 gsimple_entry->pdg = FLUKAtoPDG[FlukaID];
170 gsimple_entry->wgt = weight;
171 gsimple_entry->vtxx = x * ShipUnit::cm / ShipUnit::m; // in m
172 gsimple_entry->vtxy = y * ShipUnit::cm / ShipUnit::m;
173 gsimple_entry->vtxz = z * ShipUnit::cm / ShipUnit::m;
174 gsimple_entry->dist = 0.; // Distance from hadron decay point to neutrino "vertex", to use for oscillations,
175 // for example. Don't use.
176
177 // I'm assuming x_cos, y_cos are normalized.
178 double z_cos = sqrt(1 - pow(x_cos, 2) - pow(y_cos, 2));
179 // And massless neutrinos
180 gsimple_entry->px = Ekin * x_cos; // in GeV/c
181 gsimple_entry->py = Ekin * y_cos;
182 gsimple_entry->pz = Ekin * z_cos;
183 gsimple_entry->E = Ekin;
184
185 // Set auxiliary data
186 aux_entry->auxint.push_back(FlukaRun);
187 aux_entry->auxint.push_back(evtNumber);
188 aux_entry->auxint.push_back(FlukaID);
189 aux_entry->auxint.push_back(generationN);
190 aux_entry->auxint.push_back(last_decay_ID);
191 aux_entry->auxint.push_back(last_interaction_ID);
192 aux_entry->auxdbl.push_back(age);
193 aux_entry->auxdbl.push_back(weight);
194 aux_entry->auxdbl.push_back(Ekin);
195 aux_entry->auxdbl.push_back(x);
196 aux_entry->auxdbl.push_back(y);
197 aux_entry->auxdbl.push_back(z);
198 aux_entry->auxdbl.push_back(x_cos);
199 aux_entry->auxdbl.push_back(y_cos);
200 aux_entry->auxdbl.push_back(last_decay_x);
201 aux_entry->auxdbl.push_back(last_decay_y);
202 aux_entry->auxdbl.push_back(last_decay_z);
203 aux_entry->auxdbl.push_back(Ekin_decay);
204 aux_entry->auxdbl.push_back(last_interaction_x);
205 aux_entry->auxdbl.push_back(last_interaction_y);
206 aux_entry->auxdbl.push_back(last_interaction_z);
207 aux_entry->auxdbl.push_back(Ekin_parent);
208
209 // Accumulate metadata
210 pdglist.insert(gsimple_entry->pdg);
211 min_weight = TMath::Min(min_weight, gsimple_entry->wgt);
212 max_weight = TMath::Max(max_weight, gsimple_entry->wgt);
213 max_energy = TMath::Max(max_energy, gsimple_entry->E);
214
215 // All done!
216 tOut->Fill();
217 }
218 }
219 }
220
221 // plane corner and plane direction of the scoring plane (The scoring plane is tilted)
222 double plane_corner[] = {min_x, min_y, min_z};
223 double plane_dir1[] = {max_x - min_x, 0, max_z - min_z};
224 double plane_dir2[] = {0, max_y - min_y, max_z - min_z};
225
226 // Sort out metadata
227 // Copy pdg list
228 meta_entry->pdglist.clear();
229 for (std::set<int>::const_iterator pdg_iterator = pdglist.begin(); pdg_iterator != pdglist.end(); ++pdg_iterator)
230 meta_entry->pdglist.push_back(*pdg_iterator);
231 meta_entry->maxEnergy = max_energy;
232 meta_entry->maxWgt = max_weight;
233 meta_entry->minWgt = min_weight;
234
235 meta_entry->protons = pp_collision_number; // Number of pp collisions.
236 for (int i = 0; i < 3; i++)
237 meta_entry->windowBase[i] = plane_corner[i] * ShipUnit::cm / ShipUnit::m;
238 for (int i = 0; i < 3; i++)
239 meta_entry->windowDir1[i] = plane_dir1[i] * ShipUnit::cm / ShipUnit::m;
240 for (int i = 0; i < 3; i++)
241 meta_entry->windowDir2[i] = plane_dir2[i] * ShipUnit::cm / ShipUnit::m;
242
243 meta_entry->infiles.push_back(inFileName);
244 meta_entry->seed = ran->GetSeed();
245 meta_entry->metakey = metakey;
246
247 meta_entry->auxintname.push_back("FlukaRun");
248 meta_entry->auxintname.push_back("evtNumber");
249 meta_entry->auxintname.push_back("FlukaID");
250 meta_entry->auxintname.push_back("generationN");
251 meta_entry->auxintname.push_back("Last decay ID");
252 meta_entry->auxintname.push_back("Last interaction ID");
253 meta_entry->auxdblname.push_back("age");
254 meta_entry->auxdblname.push_back("FLUKA_weight");
255 meta_entry->auxdblname.push_back("Kenetic energy");
256 meta_entry->auxdblname.push_back("X coordinate");
257 meta_entry->auxdblname.push_back("Y coordinate");
258 meta_entry->auxdblname.push_back("X cosine");
259 meta_entry->auxdblname.push_back("Y cosine");
260 meta_entry->auxdblname.push_back("Last decay X coordinate");
261 meta_entry->auxdblname.push_back("Last decay Y coordinate");
262 meta_entry->auxdblname.push_back("Last decay Z coordinate");
263 meta_entry->auxdblname.push_back("Kenetic Energy of decay parent");
264 meta_entry->auxdblname.push_back("Last interaction in Z coordinate");
265 meta_entry->auxdblname.push_back("Last in interaction Y coordinate");
266 meta_entry->auxdblname.push_back("Last in interaction Z coordinate");
267 meta_entry->auxdblname.push_back("Kenetic energy of parent");
268
269 metaOut->Fill();
270
271 fOut->cd();
272 metaOut->Write();
273 tOut->Write();
274 fOut->Close();
275 std::cout << "Done converting" << std::endl;
276}
counter
Int_t counter
Definition Pythia8Generator.cxx:14
main
int main()
Definition main.cc:133