NuWroGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 //
4 // NuWro neutrino event generator (parser)
5 //
6 // echurch@fnal.gov
7 //
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <sys/stat.h>
11 #include <cstdlib>
12 #include <string>
13 #include <iostream>
14 #include <iomanip>
15 #include <memory>
16 #include <unistd.h>
17 #include <stdio.h>
18 #include <fstream>
19 #include <iostream>
20 
21 // Framework includes
22 #include "art/Framework/Core/EDProducer.h"
23 #include "art/Framework/Core/ModuleMacros.h"
24 #include "art/Framework/Principal/Event.h"
25 #include "fhiclcpp/ParameterSet.h"
26 #include "art/Framework/Principal/Handle.h"
27 #include "art/Framework/Services/Registry/ServiceHandle.h"
28 #include "art_root_io/TFileService.h"
29 #include "art_root_io/TFileDirectory.h"
30 #include "messagefacility/MessageLogger/MessageLogger.h"
31 
32 // ROOT includes
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "TDatabasePDG.h"
36 #include "TSystem.h"
37 #include "TROOT.h"
38 #include "TBranch.h"
39 #include "TTree.h"
40 #include "TFile.h"
41 #include "TChain.h"
42 #include "TStopwatch.h"
43 
44 
45 // LArSoft includes
46 #include "nusimdata/SimulationBase/MCTruth.h"
47 #include "nusimdata/SimulationBase/MCParticle.h"
48 #include "larcore/Geometry/Geometry.h"
49 #include "larcoreobj/SummaryData/RunData.h"
50 
51 // #include "NWtree.h"
52 
53 class TH1F;
54 class TH2F;
55 
56 // NuWro include - event1dict.h includes event1.h,
57 // the definition for the NuWro event class
58 //#include "event1dict.h"
59 
60 namespace simb { class MCTruth; }
61 
62 namespace evgen {
63 
64 
65  /// A module to check the results from the Monte Carlo generator
66  class NuWroGen : public art::EDProducer {
67  public:
68  explicit NuWroGen(fhicl::ParameterSet const &pset);
69  virtual ~NuWroGen();
70 
71  void produce(art::Event& evt);
72  void beginJob();
73  void beginRun(art::Run& run);
74  void reconfigure(fhicl::ParameterSet const& p);
75  void endJob();
76 
77  private:
78 
79  std::string ParticleStatus(int StatusCode);
80  std::string ReactionChannel(int ccnc,int mode);
81 
82  void FillHistograms(const simb::MCTruth &mc);
83 
84  std::string fFileName;
85  std::string fTreeName;
86  std::ifstream *fEventFile;
87 
88  TStopwatch fStopwatch; ///keep track of how long it takes to run the job
89 
90  std::string fNuWroModuleLabel;
91  int fEventNumberOffset;
92 
93  std::vector<double> fxsecFluxWtd;
94  double fxsecTotal;
95 
96  TH1F* fGenerated[6]; ///< Spectra as generated
97 
98  TH1F* fVertexX; ///< vertex location of generated events in x
99  TH1F* fVertexY; ///< vertex location of generated events in y
100  TH1F* fVertexZ; ///< vertex location of generated events in z
101 
102  TH2F* fVertexXY; ///< vertex location in xy
103  TH2F* fVertexXZ; ///< vertex location in xz
104  TH2F* fVertexYZ; ///< vertex location in yz
105 
106  TH1F* fDCosX; ///< direction cosine in x
107  TH1F* fDCosY; ///< direction cosine in y
108  TH1F* fDCosZ; ///< direction cosine in z
109 
110  TH1F* fMuMomentum; ///< momentum of outgoing muons
111  TH1F* fMuDCosX; ///< direction cosine of outgoing mu in x
112  TH1F* fMuDCosY; ///< direction cosine of outgoing mu in y
113  TH1F* fMuDCosZ; ///< direction cosine of outgoing mu in z
114 
115  TH1F* fEMomentum; ///< momentum of outgoing electrons
116  TH1F* fEDCosX; ///< direction cosine of outgoing e in x
117  TH1F* fEDCosY; ///< direction cosine of outgoing e in y
118  TH1F* fEDCosZ; ///< direction cosine of outgoing e in z
119 
120  TH1F* fCCMode; ///< CC interaction mode
121  TH1F* fNCMode; ///< CC interaction mode
122  TH1F* fDyn; ///< mode in detail
123  TH1F* fWeight; ///< NuWro Wt
124  TH1F* fWeightNW; ///< NuWro Wt
125  TH1F* fDynNew;
126  TH1F* fDynNewThresh;
127  TH2F* f2DynNew;
128  TH2F* f2DynNewThresh;
129 
130  // for c2: fDeltaE is no longer used
131  //TH1F* fDeltaE; ///< difference in neutrino energy from MCTruth::Enu() vs TParticle
132  TH1F* fECons; ///< histogram to determine if energy is conserved in the event
133 
134 
135 
136  unsigned int countFile;
137  event *NuWroTTree;
138  TChain *ch;
139 
140  // here come all the NuWro variables...
141  /*
142  TTree *fChain; //!pointer to the analyzed TTree or TChain
143  Int_t fCurrent; //!current Tree number in a TChain
144 
145  UInt_t fUniqueID;
146 
147  UInt_t fBits;
148 
149  Bool_t flag_coh;
150 
151  Bool_t flag_qel;
152 
153  Bool_t flag_dis;
154 
155  Bool_t flag_res;
156 
157  Bool_t flag_nc;
158 
159  Bool_t flag_cc;
160 
161  Bool_t flag_anty;
162 
163  Int_t par_random_seed;
164 
165  Int_t par_number_of_events;
166 
167  Int_t par_number_of_test_events;
168 
169  Int_t par_user_events;
170 
171  Int_t par_beam_type;
172 
173  std::string par_beam_energy_string;
174 
175  Int_t par_beam_particle;
176 
177  Double_t par_beam_direction_x;
178 
179  Double_t par_beam_direction_y;
180 
181  Double_t par_beam_direction_z;
182 
183  std::string par_beam_content_string;
184 
185  std::string par_beam_folder;
186 
187  Int_t par_beam_file_first;
188 
189  Int_t par_beam_file_limit;
190 
191  Bool_t par_beam_weighted;
192 
193  std::string par_beam_file;
194 
195  Double_t par_beam_offset_x;
196 
197  Double_t par_beam_offset_y;
198 
199  Double_t par_beam_offset_z;
200 
201  Int_t par_beam_placement;
202 
203  Int_t par_beam_test_only;
204 
205  Int_t par_nucleus_p;
206 
207  Int_t par_nucleus_n;
208 
209  std::string par_nucleus_density;
210 
211  Double_t par_nucleus_E_b;
212 
213  Double_t par_nucleus_kf;
214 
215  Int_t par_nucleus_target;
216 
217  Int_t par_nucleus_model;
218 
219  Int_t par_target_type;
220 
221  std::string par_target_content_string;
222 
223  std::string par_geo_file;
224 
225  std::string par_geo_name;
226 
227  std::string par_geo_volume;
228 
229  Double_t par_geo_o_x;
230 
231  Double_t par_geo_o_y;
232 
233  Double_t par_geo_o_z;
234 
235  Double_t par_geo_d_x;
236 
237  Double_t par_geo_d_y;
238 
239  Double_t par_geo_d_z;
240 
241  Bool_t par_dyn_qel_cc;
242 
243  Bool_t par_dyn_qel_nc;
244 
245  Bool_t par_dyn_res_cc;
246 
247  Bool_t par_dyn_res_nc;
248 
249  Bool_t par_dyn_dis_cc;
250 
251  Bool_t par_dyn_dis_nc;
252 
253  Bool_t par_dyn_coh_cc;
254 
255  Bool_t par_dyn_coh_nc;
256 
257  Int_t par_qel_vector_ff_set;
258 
259  Int_t par_qel_axial_ff_set;
260 
261  Double_t par_qel_cc_vector_mass;
262 
263  Double_t par_qel_cc_axial_mass;
264 
265  Double_t par_qel_nc_axial_mass;
266 
267  Double_t par_qel_s_axial_mass;
268 
269  Int_t par_qel_strange;
270 
271  Int_t par_qel_strangeEM;
272 
273  Double_t par_delta_s;
274 
275  Bool_t par_flux_correction;
276 
277  Int_t par_delta_FF_set;
278 
279  Double_t par_pion_axial_mass;
280 
281  Double_t par_pion_C5A;
282 
283  Double_t par_res_dis_cut;
284 
285  Int_t par_spp_precision;
286 
287  Bool_t par_coh_mass_correction;
288 
289  Bool_t par_coh_new;
290 
291  Bool_t par_kaskada_on;
292 
293  Bool_t par_kaskada_newangle;
294 
295  Bool_t par_kaskada_redo;
296 
297  Bool_t par_pauli_blocking;
298 
299  std::string par_formation_zone;
300 
301  Bool_t par_first_step;
302 
303  Double_t par_step;
304 
305  Int_t par_xsec;
306 
307  Int_t par_sf_method;
308 
309  Bool_t par_cc_smoothing;
310 
311  Bool_t par_mixed_order;
312 
313  Double_t par_rmin;
314 
315  Double_t par_rmax;
316 
317  std::string par_path_to_data;
318 
319  Int_t in_;
320 
321  Double_t in_t[kMaxin]; //[in_]
322 
323  Double_t in_x[kMaxin]; //[in_]
324 
325  Double_t in_y[kMaxin]; //[in_]
326 
327  Double_t in_z[kMaxin]; //[in_]
328 
329  Double_t in__mass[kMaxin]; //[in_]
330 
331  Double_t in_r_t[kMaxin]; //[in_]
332 
333  Double_t in_r_x[kMaxin]; //[in_]
334 
335  Double_t in_r_y[kMaxin]; //[in_]
336 
337  Double_t in_r_z[kMaxin]; //[in_]
338 
339  Int_t in_pdg[kMaxin]; //[in_]
340 
341  Char_t in_ks[kMaxin]; //[in_]
342 
343  Char_t in_orgin[kMaxin]; //[in_]
344 
345  Double_t in_travelled[kMaxin]; //[in_]
346 
347  Int_t in_id[kMaxin]; //[in_]
348 
349  Int_t in_mother[kMaxin]; //[in_]
350 
351  Int_t in_endproc[kMaxin]; //[in_]
352 
353  Double_t in_fz[kMaxin]; //[in_]
354 
355  Double_t in_pt[kMaxin]; //[in_]
356 
357  Int_t in_mother_pdg[kMaxin]; //[in_]
358 
359  Int_t in_mother_proc[kMaxin]; //[in_]
360 
361  Double_t in_mother_momentum_x[kMaxin]; //[in_]
362 
363  Double_t in_mother_momentum_y[kMaxin]; //[in_]
364 
365  Double_t in_mother_momentum_z[kMaxin]; //[in_]
366 
367  Double_t in_mother_ek[kMaxin]; //[in_]
368 
369  Int_t in_his_nqel[kMaxin]; //[in_]
370 
371  Int_t in_his_nspp[kMaxin]; //[in_]
372 
373  Int_t in_his_ndpp[kMaxin]; //[in_]
374 
375  Int_t in_his_pqel[kMaxin]; //[in_]
376 
377  Int_t in_his_pcex[kMaxin]; //[in_]
378 
379  Int_t in_his_pspp[kMaxin]; //[in_]
380 
381  Int_t in_his_pdpp[kMaxin]; //[in_]
382 
383  Int_t in_his_ptpp[kMaxin]; //[in_]
384 
385  Int_t temp_;
386 
387  Double_t temp_t[kMaxtemp]; //[temp_]
388 
389  Double_t temp_x[kMaxtemp]; //[temp_]
390 
391  Double_t temp_y[kMaxtemp]; //[temp_]
392 
393  Double_t temp_z[kMaxtemp]; //[temp_]
394 
395  Double_t temp__mass[kMaxtemp]; //[temp_]
396 
397  Double_t temp_r_t[kMaxtemp]; //[temp_]
398 
399  Double_t temp_r_x[kMaxtemp]; //[temp_]
400 
401  Double_t temp_r_y[kMaxtemp]; //[temp_]
402 
403  Double_t temp_r_z[kMaxtemp]; //[temp_]
404 
405  Int_t temp_pdg[kMaxtemp]; //[temp_]
406 
407  Char_t temp_ks[kMaxtemp]; //[temp_]
408 
409  Char_t temp_orgin[kMaxtemp]; //[temp_]
410 
411  Double_t temp_travelled[kMaxtemp]; //[temp_]
412 
413  Int_t temp_id[kMaxtemp]; //[temp_]
414 
415  Int_t temp_mother[kMaxtemp]; //[temp_]
416 
417  Int_t temp_endproc[kMaxtemp]; //[temp_]
418 
419  Double_t temp_fz[kMaxtemp]; //[temp_]
420 
421  Double_t temp_pt[kMaxtemp]; //[temp_]
422 
423  Int_t temp_mother_pdg[kMaxtemp]; //[temp_]
424 
425  Int_t temp_mother_proc[kMaxtemp]; //[temp_]
426 
427  Double_t temp_mother_momentum_x[kMaxtemp]; //[temp_]
428 
429  Double_t temp_mother_momentum_y[kMaxtemp]; //[temp_]
430 
431  Double_t temp_mother_momentum_z[kMaxtemp]; //[temp_]
432 
433  Double_t temp_mother_ek[kMaxtemp]; //[temp_]
434 
435  Int_t temp_his_nqel[kMaxtemp]; //[temp_]
436 
437  Int_t temp_his_nspp[kMaxtemp]; //[temp_]
438 
439  Int_t temp_his_ndpp[kMaxtemp]; //[temp_]
440 
441  Int_t temp_his_pqel[kMaxtemp]; //[temp_]
442 
443  Int_t temp_his_pcex[kMaxtemp]; //[temp_]
444 
445  Int_t temp_his_pspp[kMaxtemp]; //[temp_]
446 
447  Int_t temp_his_pdpp[kMaxtemp]; //[temp_]
448 
449  Int_t temp_his_ptpp[kMaxtemp]; //[temp_]
450 
451  Int_t out_;
452 
453  Double_t out_t[kMaxout]; //[out_]
454 
455  Double_t out_x[kMaxout]; //[out_]
456 
457  Double_t out_y[kMaxout]; //[out_]
458 
459  Double_t out_z[kMaxout]; //[out_]
460 
461  Double_t out__mass[kMaxout]; //[out_]
462 
463  Double_t out_r_t[kMaxout]; //[out_]
464 
465  Double_t out_r_x[kMaxout]; //[out_]
466 
467  Double_t out_r_y[kMaxout]; //[out_]
468 
469  Double_t out_r_z[kMaxout]; //[out_]
470 
471  Int_t out_pdg[kMaxout]; //[out_]
472 
473  Char_t out_ks[kMaxout]; //[out_]
474 
475  Char_t out_orgin[kMaxout]; //[out_]
476 
477  Double_t out_travelled[kMaxout]; //[out_]
478 
479  Int_t out_id[kMaxout]; //[out_]
480 
481  Int_t out_mother[kMaxout]; //[out_]
482 
483  Int_t out_endproc[kMaxout]; //[out_]
484 
485  Double_t out_fz[kMaxout]; //[out_]
486 
487  Double_t out_pt[kMaxout]; //[out_]
488 
489  Int_t out_mother_pdg[kMaxout]; //[out_]
490 
491  Int_t out_mother_proc[kMaxout]; //[out_]
492 
493  Double_t out_mother_momentum_x[kMaxout]; //[out_]
494 
495  Double_t out_mother_momentum_y[kMaxout]; //[out_]
496 
497  Double_t out_mother_momentum_z[kMaxout]; //[out_]
498 
499  Double_t out_mother_ek[kMaxout]; //[out_]
500 
501  Int_t out_his_nqel[kMaxout]; //[out_]
502 
503  Int_t out_his_nspp[kMaxout]; //[out_]
504 
505  Int_t out_his_ndpp[kMaxout]; //[out_]
506 
507  Int_t out_his_pqel[kMaxout]; //[out_]
508 
509  Int_t out_his_pcex[kMaxout]; //[out_]
510 
511  Int_t out_his_pspp[kMaxout]; //[out_]
512 
513  Int_t out_his_pdpp[kMaxout]; //[out_]
514 
515  Int_t out_his_ptpp[kMaxout]; //[out_]
516 
517  Int_t post_;
518 
519  Double_t post_t[kMaxpost]; //[post_]
520 
521  Double_t post_x[kMaxpost]; //[post_]
522 
523  Double_t post_y[kMaxpost]; //[post_]
524 
525  Double_t post_z[kMaxpost]; //[post_]
526 
527  Double_t post__mass[kMaxpost]; //[post_]
528 
529  Double_t post_r_t[kMaxpost]; //[post_]
530 
531  Double_t post_r_x[kMaxpost]; //[post_]
532 
533  Double_t post_r_y[kMaxpost]; //[post_]
534 
535  Double_t post_r_z[kMaxpost]; //[post_]
536 
537  Int_t post_pdg[kMaxpost]; //[post_]
538 
539  Char_t post_ks[kMaxpost]; //[post_]
540 
541  Char_t post_orgin[kMaxpost]; //[post_]
542 
543  Double_t post_travelled[kMaxpost]; //[post_]
544 
545  Int_t post_id[kMaxpost]; //[post_]
546 
547  Int_t post_mother[kMaxpost]; //[post_]
548 
549  Int_t post_endproc[kMaxpost]; //[post_]
550 
551  Double_t post_fz[kMaxpost]; //[post_]
552 
553  Double_t post_pt[kMaxpost]; //[post_]
554 
555  Int_t post_mother_pdg[kMaxpost]; //[post_]
556 
557  Int_t post_mother_proc[kMaxpost]; //[post_]
558 
559  Double_t post_mother_momentum_x[kMaxpost]; //[post_]
560 
561  Double_t post_mother_momentum_y[kMaxpost]; //[post_]
562 
563  Double_t post_mother_momentum_z[kMaxpost]; //[post_]
564 
565  Double_t post_mother_ek[kMaxpost]; //[post_]
566 
567  Int_t post_his_nqel[kMaxpost]; //[post_]
568 
569  Int_t post_his_nspp[kMaxpost]; //[post_]
570 
571  Int_t post_his_ndpp[kMaxpost]; //[post_]
572 
573  Int_t post_his_pqel[kMaxpost]; //[post_]
574 
575  Int_t post_his_pcex[kMaxpost]; //[post_]
576 
577  Int_t post_his_pspp[kMaxpost]; //[post_]
578 
579  Int_t post_his_pdpp[kMaxpost]; //[post_]
580 
581  Int_t post_his_ptpp[kMaxpost]; //[post_]
582 
583  Int_t all_;
584 
585  Double_t all_t[kMaxall]; //[all_]
586 
587  Double_t all_x[kMaxall]; //[all_]
588 
589  Double_t all_y[kMaxall]; //[all_]
590 
591  Double_t all_z[kMaxall]; //[all_]
592 
593  Double_t all__mass[kMaxall]; //[all_]
594 
595  Double_t all_r_t[kMaxall]; //[all_]
596 
597  Double_t all_r_x[kMaxall]; //[all_]
598 
599  Double_t all_r_y[kMaxall]; //[all_]
600 
601  Double_t all_r_z[kMaxall]; //[all_]
602 
603  Int_t all_pdg[kMaxall]; //[all_]
604 
605  Char_t all_ks[kMaxall]; //[all_]
606 
607  Char_t all_orgin[kMaxall]; //[all_]
608 
609  Double_t all_travelled[kMaxall]; //[all_]
610 
611  Int_t all_id[kMaxall]; //[all_]
612 
613  Int_t all_mother[kMaxall]; //[all_]
614 
615  Int_t all_endproc[kMaxall]; //[all_]
616 
617  Double_t all_fz[kMaxall]; //[all_]
618 
619  Double_t all_pt[kMaxall]; //[all_]
620 
621  Int_t all_mother_pdg[kMaxall]; //[all_]
622 
623  Int_t all_mother_proc[kMaxall]; //[all_]
624 
625  Double_t all_mother_momentum_x[kMaxall]; //[all_]
626 
627  Double_t all_mother_momentum_y[kMaxall]; //[all_]
628 
629  Double_t all_mother_momentum_z[kMaxall]; //[all_]
630 
631  Double_t all_mother_ek[kMaxall]; //[all_]
632 
633  Int_t all_his_nqel[kMaxall]; //[all_]
634 
635  Int_t all_his_nspp[kMaxall]; //[all_]
636 
637  Int_t all_his_ndpp[kMaxall]; //[all_]
638 
639  Int_t all_his_pqel[kMaxall]; //[all_]
640 
641  Int_t all_his_pcex[kMaxall]; //[all_]
642 
643  Int_t all_his_pspp[kMaxall]; //[all_]
644 
645  Int_t all_his_pdpp[kMaxall]; //[all_]
646 
647  Int_t all_his_ptpp[kMaxall]; //[all_]
648 
649  Double_t weight;
650 
651  Double_t norm;
652 
653  Double_t r_x;
654 
655  Double_t r_y;
656 
657  Double_t r_z;
658 
659  Double_t density;
660 
661  Int_t dyn;
662 
663  Int_t nod[12];
664 
665  Int_t place[11][20];
666 
667  Double_t pabsen;
668 
669  Bool_t nopp;
670 
671  Bool_t abs;
672 
673  Int_t nofi;
674 
675  Double_t pen[10];
676 
677  Double_t absr;
678 
679  Double_t odl;
680 
681  Double_t density_hist[50];
682 
683  Double_t radius_hist;
684 
685  Int_t protons_hist;
686 
687  Int_t neutrons_hist;
688 
689  Int_t pr;
690 
691  Int_t nr;
692 
693 
694  // List of branches
695 
696  TBranch *b_e_fUniqueID; //!
697 
698  TBranch *b_e_fBits; //!
699 
700  TBranch *b_e_flag_coh; //!
701 
702  TBranch *b_e_flag_qel; //!
703 
704  TBranch *b_e_flag_dis; //!
705 
706  TBranch *b_e_flag_res; //!
707 
708  TBranch *b_e_flag_nc; //!
709 
710  TBranch *b_e_flag_cc; //!
711 
712  TBranch *b_e_flag_anty; //!
713 
714  TBranch *b_e_par_random_seed; //!
715 
716  TBranch *b_e_par_number_of_events; //!
717 
718  TBranch *b_e_par_number_of_test_events; //!
719 
720  TBranch *b_e_par_user_events; //!
721 
722  TBranch *b_e_par_beam_type; //!
723 
724  TBranch *b_e_par_beam_energy_string; //!
725 
726  TBranch *b_e_par_beam_particle; //!
727 
728  TBranch *b_e_par_beam_direction_x; //!
729 
730  TBranch *b_e_par_beam_direction_y; //!
731 
732  TBranch *b_e_par_beam_direction_z; //!
733 
734  TBranch *b_e_par_beam_content_string; //!
735 
736  TBranch *b_e_par_beam_folder; //!
737 
738  TBranch *b_e_par_beam_file_first; //!
739 
740  TBranch *b_e_par_beam_file_limit; //!
741 
742  TBranch *b_e_par_beam_weighted; //!
743 
744  TBranch *b_e_par_beam_file; //!
745 
746  TBranch *b_e_par_beam_offset_x; //!
747 
748  TBranch *b_e_par_beam_offset_y; //!
749 
750  TBranch *b_e_par_beam_offset_z; //!
751 
752  TBranch *b_e_par_beam_placement; //!
753 
754  TBranch *b_e_par_beam_test_only; //!
755 
756  TBranch *b_e_par_nucleus_p; //!
757 
758  TBranch *b_e_par_nucleus_n; //!
759 
760  TBranch *b_e_par_nucleus_density; //!
761 
762  TBranch *b_e_par_nucleus_E_b; //!
763 
764  TBranch *b_e_par_nucleus_kf; //!
765 
766  TBranch *b_e_par_nucleus_target; //!
767 
768  TBranch *b_e_par_nucleus_model; //!
769 
770  TBranch *b_e_par_target_type; //!
771 
772  TBranch *b_e_par_target_content_string; //!
773 
774  TBranch *b_e_par_geo_file; //!
775 
776  TBranch *b_e_par_geo_name; //!
777 
778  TBranch *b_e_par_geo_volume; //!
779 
780  TBranch *b_e_par_geo_o_x; //!
781 
782  TBranch *b_e_par_geo_o_y; //!
783 
784  TBranch *b_e_par_geo_o_z; //!
785 
786  TBranch *b_e_par_geo_d_x; //!
787 
788  TBranch *b_e_par_geo_d_y; //!
789 
790  TBranch *b_e_par_geo_d_z; //!
791 
792  TBranch *b_e_par_dyn_qel_cc; //!
793 
794  TBranch *b_e_par_dyn_qel_nc; //!
795 
796  TBranch *b_e_par_dyn_res_cc; //!
797 
798  TBranch *b_e_par_dyn_res_nc; //!
799 
800  TBranch *b_e_par_dyn_dis_cc; //!
801 
802  TBranch *b_e_par_dyn_dis_nc; //!
803 
804  TBranch *b_e_par_dyn_coh_cc; //!
805 
806  TBranch *b_e_par_dyn_coh_nc; //!
807 
808  TBranch *b_e_par_qel_vector_ff_set; //!
809 
810  TBranch *b_e_par_qel_axial_ff_set; //!
811 
812  TBranch *b_e_par_qel_cc_vector_mass; //!
813 
814  TBranch *b_e_par_qel_cc_axial_mass; //!
815 
816  TBranch *b_e_par_qel_nc_axial_mass; //!
817 
818  TBranch *b_e_par_qel_s_axial_mass; //!
819 
820  TBranch *b_e_par_qel_strange; //!
821 
822  TBranch *b_e_par_qel_strangeEM; //!
823 
824  TBranch *b_e_par_delta_s; //!
825 
826  TBranch *b_e_par_flux_correction; //!
827 
828  TBranch *b_e_par_delta_FF_set; //!
829 
830  TBranch *b_e_par_pion_axial_mass; //!
831 
832  TBranch *b_e_par_pion_C5A; //!
833 
834  TBranch *b_e_par_res_dis_cut; //!
835 
836  TBranch *b_e_par_spp_precision; //!
837 
838  TBranch *b_e_par_coh_mass_correction; //!
839 
840  TBranch *b_e_par_coh_new; //!
841 
842  TBranch *b_e_par_kaskada_on; //!
843 
844  TBranch *b_e_par_kaskada_newangle; //!
845 
846  TBranch *b_e_par_kaskada_redo; //!
847 
848  TBranch *b_e_par_pauli_blocking; //!
849 
850  TBranch *b_e_par_formation_zone; //!
851 
852  TBranch *b_e_par_first_step; //!
853 
854  TBranch *b_e_par_step; //!
855 
856  TBranch *b_e_par_xsec; //!
857 
858  TBranch *b_e_par_sf_method; //!
859 
860  TBranch *b_e_par_cc_smoothing; //!
861 
862  TBranch *b_e_par_mixed_order; //!
863 
864  TBranch *b_e_par_rmin; //!
865 
866  TBranch *b_e_par_rmax; //!
867 
868  TBranch *b_e_par_path_to_data; //!
869 
870  TBranch *b_e_in_; //!
871 
872  TBranch *b_in_t; //!
873 
874  TBranch *b_in_x; //!
875 
876  TBranch *b_in_y; //!
877 
878  TBranch *b_in_z; //!
879 
880  TBranch *b_in__mass; //!
881 
882  TBranch *b_in_r_t; //!
883 
884  TBranch *b_in_r_x; //!
885 
886  TBranch *b_in_r_y; //!
887 
888  TBranch *b_in_r_z; //!
889 
890  TBranch *b_in_pdg; //!
891 
892  TBranch *b_in_ks; //!
893 
894  TBranch *b_in_orgin; //!
895 
896  TBranch *b_in_travelled; //!
897 
898  TBranch *b_in_id; //!
899 
900  TBranch *b_in_mother; //!
901 
902  TBranch *b_in_endproc; //!
903 
904  TBranch *b_in_fz; //!
905 
906  TBranch *b_in_pt; //!
907 
908  TBranch *b_in_mother_pdg; //!
909 
910  TBranch *b_in_mother_proc; //!
911 
912  TBranch *b_in_mother_momentum_x; //!
913 
914  TBranch *b_in_mother_momentum_y; //!
915 
916  TBranch *b_in_mother_momentum_z; //!
917 
918  TBranch *b_in_mother_ek; //!
919 
920  TBranch *b_in_his_nqel; //!
921 
922  TBranch *b_in_his_nspp; //!
923 
924  TBranch *b_in_his_ndpp; //!
925 
926  TBranch *b_in_his_pqel; //!
927 
928  TBranch *b_in_his_pcex; //!
929 
930  TBranch *b_in_his_pspp; //!
931 
932  TBranch *b_in_his_pdpp; //!
933 
934  TBranch *b_in_his_ptpp; //!
935 
936  TBranch *b_e_temp_; //!
937 
938  TBranch *b_temp_t; //!
939 
940  TBranch *b_temp_x; //!
941 
942  TBranch *b_temp_y; //!
943 
944  TBranch *b_temp_z; //!
945 
946  TBranch *b_temp__mass; //!
947 
948  TBranch *b_temp_r_t; //!
949 
950  TBranch *b_temp_r_x; //!
951 
952  TBranch *b_temp_r_y; //!
953 
954  TBranch *b_temp_r_z; //!
955 
956  TBranch *b_temp_pdg; //!
957 
958  TBranch *b_temp_ks; //!
959 
960  TBranch *b_temp_orgin; //!
961 
962  TBranch *b_temp_travelled; //!
963 
964  TBranch *b_temp_id; //!
965 
966  TBranch *b_temp_mother; //!
967 
968  TBranch *b_temp_endproc; //!
969 
970  TBranch *b_temp_fz; //!
971 
972  TBranch *b_temp_pt; //!
973 
974  TBranch *b_temp_mother_pdg; //!
975 
976  TBranch *b_temp_mother_proc; //!
977 
978  TBranch *b_temp_mother_momentum_x; //!
979 
980  TBranch *b_temp_mother_momentum_y; //!
981 
982  TBranch *b_temp_mother_momentum_z; //!
983 
984  TBranch *b_temp_mother_ek; //!
985 
986  TBranch *b_temp_his_nqel; //!
987 
988  TBranch *b_temp_his_nspp; //!
989 
990  TBranch *b_temp_his_ndpp; //!
991 
992  TBranch *b_temp_his_pqel; //!
993 
994  TBranch *b_temp_his_pcex; //!
995 
996  TBranch *b_temp_his_pspp; //!
997 
998  TBranch *b_temp_his_pdpp; //!
999 
1000  TBranch *b_temp_his_ptpp; //!
1001 
1002  TBranch *b_e_out_; //!
1003 
1004  TBranch *b_out_t; //!
1005 
1006  TBranch *b_out_x; //!
1007 
1008  TBranch *b_out_y; //!
1009 
1010  TBranch *b_out_z; //!
1011 
1012  TBranch *b_out__mass; //!
1013 
1014  TBranch *b_out_r_t; //!
1015 
1016  TBranch *b_out_r_x; //!
1017 
1018  TBranch *b_out_r_y; //!
1019 
1020  TBranch *b_out_r_z; //!
1021 
1022  TBranch *b_out_pdg; //!
1023 
1024  TBranch *b_out_ks; //!
1025 
1026  TBranch *b_out_orgin; //!
1027 
1028  TBranch *b_out_travelled; //!
1029 
1030  TBranch *b_out_id; //!
1031 
1032  TBranch *b_out_mother; //!
1033 
1034  TBranch *b_out_endproc; //!
1035 
1036  TBranch *b_out_fz; //!
1037 
1038  TBranch *b_out_pt; //!
1039 
1040  TBranch *b_out_mother_pdg; //!
1041 
1042  TBranch *b_out_mother_proc; //!
1043 
1044  TBranch *b_out_mother_momentum_x; //!
1045 
1046  TBranch *b_out_mother_momentum_y; //!
1047 
1048  TBranch *b_out_mother_momentum_z; //!
1049 
1050  TBranch *b_out_mother_ek; //!
1051 
1052  TBranch *b_out_his_nqel; //!
1053 
1054  TBranch *b_out_his_nspp; //!
1055 
1056  TBranch *b_out_his_ndpp; //!
1057 
1058  TBranch *b_out_his_pqel; //!
1059 
1060  TBranch *b_out_his_pcex; //!
1061 
1062  TBranch *b_out_his_pspp; //!
1063 
1064  TBranch *b_out_his_pdpp; //!
1065 
1066  TBranch *b_out_his_ptpp; //!
1067 
1068  TBranch *b_e_post_; //!
1069 
1070  TBranch *b_post_t; //!
1071 
1072  TBranch *b_post_x; //!
1073 
1074  TBranch *b_post_y; //!
1075 
1076  TBranch *b_post_z; //!
1077 
1078  TBranch *b_post__mass; //!
1079 
1080  TBranch *b_post_r_t; //!
1081 
1082  TBranch *b_post_r_x; //!
1083 
1084  TBranch *b_post_r_y; //!
1085 
1086  TBranch *b_post_r_z; //!
1087 
1088  TBranch *b_post_pdg; //!
1089 
1090  TBranch *b_post_ks; //!
1091 
1092  TBranch *b_post_orgin; //!
1093 
1094  TBranch *b_post_travelled; //!
1095 
1096  TBranch *b_post_id; //!
1097 
1098  TBranch *b_post_mother; //!
1099 
1100  TBranch *b_post_endproc; //!
1101 
1102  TBranch *b_post_fz; //!
1103 
1104  TBranch *b_post_pt; //!
1105 
1106  TBranch *b_post_mother_pdg; //!
1107 
1108  TBranch *b_post_mother_proc; //!
1109 
1110  TBranch *b_post_mother_momentum_x; //!
1111 
1112  TBranch *b_post_mother_momentum_y; //!
1113 
1114  TBranch *b_post_mother_momentum_z; //!
1115 
1116  TBranch *b_post_mother_ek; //!
1117 
1118  TBranch *b_post_his_nqel; //!
1119 
1120  TBranch *b_post_his_nspp; //!
1121 
1122  TBranch *b_post_his_ndpp; //!
1123 
1124  TBranch *b_post_his_pqel; //!
1125 
1126  TBranch *b_post_his_pcex; //!
1127 
1128  TBranch *b_post_his_pspp; //!
1129 
1130  TBranch *b_post_his_pdpp; //!
1131 
1132  TBranch *b_post_his_ptpp; //!
1133 
1134  TBranch *b_e_all_; //!
1135 
1136  TBranch *b_all_t; //!
1137 
1138  TBranch *b_all_x; //!
1139 
1140  TBranch *b_all_y; //!
1141 
1142  TBranch *b_all_z; //!
1143 
1144  TBranch *b_all__mass; //!
1145 
1146  TBranch *b_all_r_t; //!
1147 
1148  TBranch *b_all_r_x; //!
1149 
1150  TBranch *b_all_r_y; //!
1151 
1152  TBranch *b_all_r_z; //!
1153 
1154  TBranch *b_all_pdg; //!
1155 
1156  TBranch *b_all_ks; //!
1157 
1158  TBranch *b_all_orgin; //!
1159 
1160  TBranch *b_all_travelled; //!
1161 
1162  TBranch *b_all_id; //!
1163 
1164  TBranch *b_all_mother; //!
1165 
1166  TBranch *b_all_endproc; //!
1167 
1168  TBranch *b_all_fz; //!
1169 
1170  TBranch *b_all_pt; //!
1171 
1172  TBranch *b_all_mother_pdg; //!
1173 
1174  TBranch *b_all_mother_proc; //!
1175 
1176  TBranch *b_all_mother_momentum_x; //!
1177 
1178  TBranch *b_all_mother_momentum_y; //!
1179 
1180  TBranch *b_all_mother_momentum_z; //!
1181 
1182  TBranch *b_all_mother_ek; //!
1183 
1184  TBranch *b_all_his_nqel; //!
1185 
1186  TBranch *b_all_his_nspp; //!
1187 
1188  TBranch *b_all_his_ndpp; //!
1189 
1190  TBranch *b_all_his_pqel; //!
1191 
1192  TBranch *b_all_his_pcex; //!
1193 
1194  TBranch *b_all_his_pspp; //!
1195 
1196  TBranch *b_all_his_pdpp; //!
1197 
1198  TBranch *b_all_his_ptpp; //!
1199 
1200  TBranch *b_e_weight; //!
1201 
1202  TBranch *b_e_norm; //!
1203 
1204  TBranch *b_e_r_x; //!
1205 
1206  TBranch *b_e_r_y; //!
1207 
1208  TBranch *b_e_r_z; //!
1209 
1210  TBranch *b_e_density; //!
1211 
1212  TBranch *b_e_dyn; //!
1213 
1214  TBranch *b_e_nod; //!
1215 
1216  TBranch *b_e_place; //!
1217 
1218  TBranch *b_e_pabsen; //!
1219 
1220  TBranch *b_e_nopp; //!
1221 
1222  TBranch *b_e_abs; //!
1223 
1224  TBranch *b_e_nofi; //!
1225 
1226  TBranch *b_e_pen; //!
1227 
1228  TBranch *b_e_absr; //!
1229 
1230  TBranch *b_e_odl; //!
1231 
1232  TBranch *b_e_density_hist; //!
1233 
1234  TBranch *b_e_radius_hist; //!
1235 
1236  TBranch *b_e_protons_hist; //!
1237 
1238  TBranch *b_e_neutrons_hist; //!
1239 
1240  TBranch *b_e_pr; //!
1241 
1242  TBranch *b_e_nr; //!
1243  */
1244 
1245  };
1246 }
1247 
1248 
1249 namespace evgen{
1250 
1251  //____________________________________________________________________________
1252  NuWroGen::NuWroGen(fhicl::ParameterSet const& pset)
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  }
1265 
1266  //____________________________________________________________________________
1267 
1268  NuWroGen::~NuWroGen()
1269  {
1270  fStopwatch.Stop();
1271  }
1272 
1273  //____________________________________________________________________________
1274  void NuWroGen::reconfigure(fhicl::ParameterSet const& p)
1275  {
1276  fFileName = p.get< std::string >("NuWroFile","output.root");
1277  fTreeName = p.get< std::string >("TreeName","treeout");
1278 
1279  return;
1280  }
1281 //___________________________________________________________________________
1282 
1283  void NuWroGen::beginJob(){
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  }
1493 
1494  //____________________________________________________________________________
1495  void NuWroGen::beginRun(art::Run& run)
1496  {
1497  art::ServiceHandle<geo::Geometry const> geo;
1498  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
1499  }
1500 
1501  //____________________________________________________________________________
1502  void NuWroGen::endJob()
1503  {
1504  delete NuWroTTree;
1505  delete ch;
1506  fxsecFluxWtd.clear();
1507  }
1508  //____________________________________________________________________________
1509  void NuWroGen::produce(art::Event& evt)
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  }
1698 
1699 // //......................................................................
1700  std::string NuWroGen::ParticleStatus(int StatusCode)
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  }
1721 
1722 
1723 // //......................................................................
1724  std::string NuWroGen::ReactionChannel(int ccnc,int mode)
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  }
1750 
1751 // //......................................................................
1752  void NuWroGen::FillHistograms(const simb::MCTruth &mc)
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  }
1987 
1988 }
1989 
1990 
1991 namespace evgen{
1992 
1993  DEFINE_ART_MODULE(NuWroGen)
1994 
1995 }
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
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
simb::MCParticle::AddTrajectoryPoint
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:257
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
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
A module to check the results from the Monte Carlo generator.
Definition: NuWroGen_module.cc:66
evgen::NuWroGen::fWeightNW
TH1F * fWeightNW
NuWro Wt.
Definition: NuWroGen_module.cc:124
evgen::NuWroGen::beginRun
void beginRun(art::Run &run)
Definition: NuWroGen_module.cc:1495
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
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
MCParticle.h
Particle class.
evgen::NuWroGen::ParticleStatus
std::string ParticleStatus(int StatusCode)
Definition: NuWroGen_module.cc:1700
Geometry.h
art framework interface to geometry description
art::Run
Definition: Run.h:17
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::produce
void produce(art::Event &evt)
Definition: NuWroGen_module.cc:1509
evgen::NuWroGen::fNuWroModuleLabel
std::string fNuWroModuleLabel
keep track of how long it takes to run the job
Definition: NuWroGen_module.cc:90
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
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
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
breakpoints::beginJob
void beginJob()
Definition: Breakpoints.cc:14
RunHistoryService::reconfigure
virtual void reconfigure(fhicl::ParameterSet const &pset)
geo::GeometryCore::DetectorName
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:1761
evgen::NuWroGen::fTreeName
std::string fTreeName
Definition: NuWroGen_module.cc:85
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
wirecell.gen.depos.move
def move(depos, offset)
Definition: depos.py:107
create_imaplist.status
int status
Definition: create_imaplist.py:5
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:234
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
vector::size
size_t size
Definition: lodepng.cpp:55
evgen::NuWroGen::fWeight
TH1F * fWeight
NuWro Wt.
Definition: NuWroGen_module.cc:123
test.p
p
Definition: test.py:223
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
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
simb
Base utilities and modules for event generation and detector simulation.
Definition: SimulationDrawer.h:16
evgen::NuWroGen::fStopwatch
TStopwatch fStopwatch
Definition: NuWroGen_module.cc:88
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::fEventFile
std::ifstream * fEventFile
Definition: NuWroGen_module.cc:86
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
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
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
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::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
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
simb::MCNeutrino
Event generator information.
Definition: MCNeutrino.h:18
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
simb::MCNeutrino::Mode
int Mode() const
Definition: MCNeutrino.h:149
evgen
Event Generation using GENIE, cosmics or single particles.
Definition: CRTGen_module.cc:33
sim::ParticleStatusName
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
Definition: MCDumperUtils.cxx:126
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
simb::kBeamNeutrino
Beam neutrinos.
Definition: MCTruth.h:23
evgen::NuWroGen::fMuDCosZ
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
Definition: NuWroGen_module.cc:113
evgen::NuWroGen::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: NuWroGen_module.cc:1274
Generated by   doxygen 1.8.11