68{
69
71
72
73
74
75
76 const int nfates = 9;
77 int countfate [nfates];
78 double sigma [nfates];
79 double sigma_err [nfates];
80 string fatestr [nfates] = " ";
82
92
93 for (int k=0; k<nfates; k++) {
94 countfate[k] = 0;
95 sigma [k] = 0.;
96 sigma_err[k] = 0.;
98 }
99
100
101 int nev = 0;
102 int probe_pdg = 0;
103 int target_pdg = 0;
104 int displayno = 100;
105 double kin_energy = 0.;
106
107
108
109
110
111 TTree * tree = 0;
112 TTree * ginuke = 0;
114
116
117 tree = dynamic_cast <TTree *> ( file.Get("gtree") );
118 ginuke = dynamic_cast <TTree *> ( file.Get("ginuke") );
120
121
122
123
124
125
126
127 if (tree) {
128
130 tree->SetBranchAddress("gmcrec", &mcrec);
131
132
133 nev = (int) tree->GetEntries();
135 << "Processing " << nev << " events";
136
137
138
139
140
141 for(int ievent = 0; ievent < nev; ievent++) {
142
143
144 tree->GetEntry(ievent);
145
146
148
149
150 if(ievent==0) {
151 kin_energy = event.Particle(0)->KinE();
152 probe_pdg = event.Particle(0)->Pdg();
153 target_pdg = event.Particle(1)->Pdg();
154 }
155
156
159 if(ievent<displayno) {
162 }
163
164
165 switch (fate){
166
167 case 0: countfate[0]++; break;
168 case 1: countfate[1]++; break;
169 case 2: countfate[2]++; break;
170 case 3: countfate[3]++; break;
171 case 4: countfate[4]++; break;
172 case 5: countfate[5]++; break;
173 case 6: countfate[6]++; break;
174 case 13: countfate[8]++; break;
175 default:
176 if (7<=fate && fate<=12) countfate[7]++;
177 else {
179 << "Undefined fate from FindhAFate() : " << fate;
180 }
181 break;
182 }
183
184
186
187 }
188 }
189 else if ( ginuke ) {
190
191
193 << "Found ginuke type file";
194
195 nev = (int) ginuke->GetEntries();
197 << "Processing " << nev << " events";
198
199 int kmax = 250;
200 int index = 0;
201 int numh = 0;
202 int numpip = 0;
203 int numpi0 = 0;
204 int numpim = 0;
205 int pdg_had[kmax];
206 double E_had [kmax];
207 double energy = 0.0;
208
209 ginuke->SetBranchAddress("ke", &kin_energy);
210 ginuke->SetBranchAddress("probe",&probe_pdg );
211 ginuke->SetBranchAddress("tgt", &target_pdg);
212 ginuke->SetBranchAddress("nh", &numh );
213 ginuke->SetBranchAddress("npip", &numpip );
214 ginuke->SetBranchAddress("npi0", &numpi0 );
215 ginuke->SetBranchAddress("npim", &numpim );
216 ginuke->SetBranchAddress("pdgh", &pdg_had );
217 ginuke->SetBranchAddress("Eh", &E_had );
218 ginuke->SetBranchAddress("e", &energy );
219
220
221 for(int ievent = 0; ievent < nev; ievent++) {
222
223
224 ginuke->GetEntry(ievent);
225
226
227 if (energy==E_had[0] && numh==1)
228 { index=1; }
229 else if (energy!=E_had[0] && numh==1)
230 { index=3; }
231 else if (
pdg::IsPion(probe_pdg) && numpip+numpi0+numpim==0)
232 { index=5; }
233 else if ( (
pdg::IsNucleon(probe_pdg) && numpip+numpi0+numpim==0 && numh>2 )
234 || (probe_pdg==
kPdgGamma && energy!=E_had[0] && numpip+numpi0+numpim==0))
235 { index=6; }
236 else if ( numpip+numpi0+numpim> (
pdg::IsPion(probe_pdg) ? 1 : 0) )
237 { index=7; }
238 else if ( numh==2 )
239 {
240 if ( (
pdg::IsPion(probe_pdg) && ( probe_pdg==pdg_had[0] || probe_pdg==pdg_had[1] ))
242 else index=2;
243 }
244 else
245 {
246 bool undef = true;
248 {
249 for (int iter = 0; iter < numh; iter++)
250 {
251 if (probe_pdg==211 && pdg_had[iter]==-211) { index=8; undef=false; }
252 else if (probe_pdg==-211 && pdg_had[iter]==211) { index=8; undef=false; }
253 }
254 }
255 if (undef) { index=0; }
256 }
257 countfate[index]++;
258 if (ievent<displayno) {
261 }
262
263 }
264 }
265 else {
267 << "Could not read input file!";
268 return 1;
269 }
270
271
272
273
274
276 const double dnev = (double) nev;
277 const int NR = 3;
278 const double R0 = 1.4;
279
282 double nuclear_radius = NR * R0 * TMath::Power(A, 1./3.);
283 double area = TMath::Pi() * TMath::Power(nuclear_radius,2);
284
286 string probe_name = pdglib->
Find(probe_pdg)->GetName();
287 string target_name = pdglib->
Find(target_pdg)->GetName();
288
290 << " Probe = " << probe_name
291 << ", KE = " << kin_energy << " GeV";
293 << " Target = " << target_name
294 << " (Z,A) = (" << Z << ", " << A
295 << "), nuclear radius = " << nuclear_radius
296 << " fm, area = " << area << " fm**2 " << '\n';
297
298 int cnttot = 0;
299 int nullint = countfate[1];
300 double sigtot = 0;
301 double sigtoterr = 0;
302 double sigtotScat = 0;
303 double sigtotAbs = 0;
304
305 for(int k=0; k<nfates; k++) {
306 if(k!=1) {
307 cnttot += countfate[k];
308 double ratio = countfate[k]/dnev;
309 sigma[k] = fm2tomb * area * ratio;
310 sigma_err[k] = fm2tomb * area * TMath::Sqrt(ratio*(1-ratio)/dnev);
311 if(sigma_err[k]==0) {
312 sigma_err[k] = fm2tomb * area * TMath::Sqrt(countfate[k])/dnev;
313 }
314 if(countfate[k]>0) {
316 << " --> " << setw(26) << fatestr[k]
317 << ": " << setw(7) << countfate[k] << " events -> "
318 << setw(7) << sigma[k] << " +- " << sigma_err[k] << " (mb)" << '\n';
319 }
320 if(k==5) {
321 sigtotAbs += sigma[k];
322 }
323 else
324 if (k!=1) {
325 sigtotScat += sigma[k];
326 }
327 }
328 }
329
330 sigtot = fm2tomb * area * cnttot/dnev;
331 sigtoterr = fm2tomb * area * TMath::Sqrt(cnttot)/dnev;
332
333 double sigtot_noelas = fm2tomb * area * (cnttot-countfate[3])/dnev;
334 double sigtoterr_noelas = fm2tomb * area * TMath::Sqrt(cnttot-countfate[3])/dnev;
335
336 double ratio_as = (sigtotScat==0) ? 0 : sigtotAbs/(double)sigtotScat;
337
339 << "\n\n --------------------------------------------------- "
340 << "\n ==> " << setw(28) << " Total: " << setw(7) << cnttot
341 << " events -> " << setw(7) << sigtot << " +- " << sigtoterr << " (mb)"
342 << "\n (-> " << setw(28) << " Hadrons escaped nucleus: "
343 << setw(7) << nullint << " events)"
344 << "\n ==> " << setw(28) << " Ratio (abs/scat) = "
345 << setw(7) << ratio_as
346 << "\n ==> " << setw(28) << " avg. num of int. = "
347 << setw(7) << cnttot/dnev
348 << "\n ==> " << setw(28) << " no interaction = "
349 << setw(7) << (dnev-cnttot)/dnev
350 << "\n ------------------------------------------------------- \n";
351
353 {
354 ifstream test_file;
358 test_file.close();
359 ofstream xsec_file;
361 if (!file_exists)
362 {
363 xsec_file << "#KE" << "\t" << "Undef" << "\t"
364 << "sig" << "\t" << "CEx" << "\t"
365 << "sig" << "\t" << "Elas" << "\t"
366 << "sig" << "\t" << "Inelas"<< "\t"
367 << "sig" << "\t" << "Abs" << "\t"
368 << "sig" << "\t" << "KO" << "\t"
369 << "sig" << "\t" << "PiPro" << "\t"
370 << "sig" << "\t" << "DCEx" << "\t"
371 << "sig" << "\t" << "Reac" << "\t"
372 << "sig" << "\t" << "Tot" << "\t" << "sig" << endl;
373 }
374 xsec_file << kin_energy;
375 for(int k=0; k<nfates; k++) {
376 if (k==1) continue;
377 xsec_file << "\t" << sigma[k] << "\t" << sigma_err[k];
378 }
379 xsec_file << "\t" << sigtot_noelas << "\t" << sigtoterr_noelas;
380 xsec_file << "\t" << sigtot << "\t" << sigtoterr << endl;
381 xsec_file.close();
382 }
383
384 return 0;
385}
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
GENIE's GHEP MC event record.
static string AsString(INukeFateHN_t fate)
MINOS-style ntuple record. Each such ntuple record holds a generated EventRecord object....
void Clear(Option_t *opt="")
Singleton class to load & serve a TDatabasePDG.
static PDGLibrary * Instance(void)
TParticlePDG * Find(int pdgc, bool must_exist=true)
void GetCommandLineArgs(int argc, char **argv)
INukeFateHA_t FindhAFate(const GHepRecord *evrec)
bool file_exists(const std::string &file_name)
Returns true if a given file exists and is accessible, or false otherwise.
int IonPdgCodeToZ(int pdgc)
int IonPdgCodeToA(int pdgc)
static constexpr double mb
static constexpr double fm2
enum genie::EINukeFateHA_t INukeFateHA_t