MergeGenCollections_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 // \file MergeGenCollections.cxx
3 // \brief
4 // \authors janow@fnal.gov
5 ///////////////////////////////////////////////////////////////////////////
6 
7 #include <math.h>
8 #include <string>
9 
16 #include "fhiclcpp/ParameterSet.h"
18 
19 #include "dk2nu/tree/dk2nu.h"
20 #include "dk2nu/tree/NuChoice.h"
21 
22 #include "DAQChannelMap/DAQChannelMap.h"
26 
27 #include "Simulation/FLSHit.h"
28 #include "Simulation/FLSHitList.h"
29 #include "Simulation/Particle.h"
30 #include "SummaryData/POTSum.h"
31 #include "SummaryData/SpillData.h"
32 #include "Utilities/AssociationUtil.h"
33 
34 
35 /// Merge generation collections MCTruth/MCFlux/GTruth/Dk2Nu/NuChoice
37 {
39  public:
40 
41  explicit MergeGenCollections(fhicl::ParameterSet const &pset);
42  virtual ~MergeGenCollections();
43 
44  void produce(art::Event& evt);
45  void endSubRun(art::SubRun& sr);
46 
47  protected:
48 
49  // configuration settings
50  std::string fFirstCollection; /// MCTruth label for main collection
51  std::string fSecondCollection; /// Process label for the second collection
54  bool fMergeDk2Nu; /// Merge Dk2Nu/NuChoice objects?
55  int fSecondaryModeOffset; /// Flag a MCNeutrino as secondary
56  };
57 
58 
59  //----------------------------------------------------------------------
61  fFirstCollection (pset.get< std::string >("FirstGenCollection","")),
62  fSecondCollection (pset.get< std::string >("SecondCollection","")),
63  fFirstColIsMCNu (!pset.get< bool >("IsFirstCollectionCosmic")),
64  fSecondColIsMCNu (!pset.get< bool >("IsSecondCollectionCosmic")),
65  fMergeDk2Nu (pset.get< bool >("MergeDk2Nu")),
66  fSecondaryModeOffset (pset.get< int >("SecondaryModeOffset",0))
67  {
68  if(fFirstCollection=="" && fSecondCollection==""){
69  std::cerr << "MergeGenCollections: "
70  << "Neither collection is specified to merge from. Abort."
71  << std::endl;
72  std::abort();
73  }
74  if(!fFirstColIsMCNu && !fSecondColIsMCNu){
75  std::cerr << "MergeGenCollections: "
76  << "No MC Collections to merge. Abort."
77  << std::endl;
78  std::abort();
79  }
80 
81  produces< sumdata::POTSum, art::InSubRun >();
82  produces< sumdata::SpillData >();
83  produces< std::vector<simb::MCTruth> >();
84  produces< std::vector<simb::MCFlux> >();
85  produces< std::vector<simb::GTruth> >();
86  produces< std::vector<bsim::Dk2Nu> >();
87  produces< std::vector<bsim::NuChoice> >();
88  produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
89  produces< art::Assns<simb::MCTruth, simb::GTruth> >();
90  }
91 
92  //----------------------------------------------------------------------
93  MergeGenCollections::~MergeGenCollections()
94  {
95  }
96 
97  //----------------------------------------------------------------------
98  void MergeGenCollections::produce(art::Event& evt)
99  {
100  std::unique_ptr< sumdata::SpillData > spillcol (new sumdata::SpillData);
101  std::unique_ptr< std::vector<simb::MCTruth> > truthcol (new std::vector<simb::MCTruth >);
102  std::unique_ptr< std::vector<simb::MCFlux> > fluxcol (new std::vector<simb::MCFlux >);
103  std::unique_ptr< std::vector<simb::GTruth> > gtruthcol(new std::vector<simb::GTruth >);
104  std::unique_ptr< std::vector<bsim::Dk2Nu> > dk2nucol (new std::vector<bsim::Dk2Nu >);
105  std::unique_ptr< std::vector<bsim::NuChoice> > nuccol (new std::vector<bsim::NuChoice>);
106 
107  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> > assns (new art::Assns<simb::MCTruth, simb::MCFlux>);
108  std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > gassns(new art::Assns<simb::MCTruth, simb::GTruth>);
109 
110 
111 
112  //// First Collection
113  //
114  if(fFirstCollection!=""){ // Dont look for MC if not there
115 
116  // Pull SpillData from the first collection only
118  evt.getByLabel(fFirstCollection, spill);
119  if(spill.isValid()) *spillcol = *spill;
120 
121  // MCTruth from generator -- not necessarily neutrino MC
123  evt.getByLabel(fFirstCollection, mclist);
124 
129  if(fFirstColIsMCNu) {
130  evt.getByLabel(fFirstCollection, fllist);
131  evt.getByLabel(fFirstCollection, glist);
132  if(fMergeDk2Nu){
133  evt.getByLabel(fFirstCollection, dklist);
134  evt.getByLabel(fFirstCollection, nulist);
135  }
136  }
137 
138  // Loop over MCTruths and associate to Flux/GTruth
139  for(unsigned int i = 0; i < mclist->size(); ++i){
140  simb::MCTruth truth ((*mclist)[i]);
141  truthcol->push_back(truth);
142 
143  if(fFirstColIsMCNu){
144  if(i < fllist->size() && i < glist->size()){
145  simb::MCFlux flux ((*fllist)[i]);
146  simb::GTruth gtruth ((*glist)[i]);
147  fluxcol ->push_back(flux);
148  gtruthcol->push_back(gtruth);
149 
150  util::CreateAssn(*this, evt, *truthcol, *fluxcol, *assns, fluxcol ->size()-1, fluxcol ->size());
151  util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size());
152  }
153  if(i < dklist->size() && i < nulist->size()){
154  bsim::Dk2Nu dk2nu ((*dklist)[i]);
155  dk2nucol->push_back(dk2nu);
156  bsim::NuChoice nuc ((*nulist)[i]);
157  nuccol ->push_back(nuc);
158  }
159  } // end fFirstColIsMCNu
160  } // end loop over MCTruths
161  } // If FirstCollection specified
162 
163 
164 
165  //// Second Collection
166  //
167  if(fSecondCollection!=""){ // dont use label if not there
169  evt.getByLabel(fSecondCollection, mc2list);
170 
175  if(fSecondColIsMCNu) {
176  evt.getByLabel(fSecondCollection, g2list);
177  evt.getByLabel(fSecondCollection, fl2list);
178  evt.getByLabel(fSecondCollection, dk2list);
179  evt.getByLabel(fSecondCollection, nu2list);
180  }
181 
182  for(unsigned int i = 0; i < mc2list->size(); ++i){
183  simb::MCTruth tru2 ((*mc2list)[i]);
184 
185  if(fSecondaryModeOffset){
186  // Make new truth so MCNeutrino mode can be
187  // offset and used to flag overlaid MCTruth
188  simb::MCNeutrino nu2 = tru2.GetNeutrino();
189  simb::MCTruth truth;
190  truth.SetOrigin(tru2.Origin());
191  for( int p=0; p<tru2.NParticles(); p++ ){
192  simb::MCParticle mcp2 = tru2.GetParticle(p);
193  truth.Add(mcp2);
194  }
195  // set nu must be after MCParticle list is made
196  truth.SetNeutrino( nu2.CCNC(),
197  nu2.Mode() + fSecondaryModeOffset,
198  nu2.InteractionType(),
199  nu2.Target(), nu2.HitNuc(),
200  nu2.HitQuark(), nu2.W(),
201  nu2.X(), nu2.Y(),nu2.QSqr() );
202  truthcol->push_back(truth);
203 
204  } else {
205  // don't bother if no offset
206  truthcol->push_back(tru2);
207  }
208 
209 
210  if(fSecondColIsMCNu) {
211  if(i < fl2list->size() && i < g2list->size()){
212  simb::MCFlux flux ((*fl2list)[i]);
213  simb::GTruth gtruth ((*g2list) [i]);
214  gtruthcol->push_back(gtruth);
215  fluxcol->push_back(flux);
216 
217  util::CreateAssn(*this, evt, *truthcol, *fluxcol, *assns, fluxcol ->size()-1, fluxcol ->size());
218  util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size());
219  }
220  if(i < dk2list->size() && i < nu2list->size()){
221  bsim::Dk2Nu dk2nu ((*dk2list)[i]);
222  dk2nucol->push_back(dk2nu);
223  bsim::NuChoice nuc ((*nu2list)[i]);
224  nuccol ->push_back(nuc);
225  }
226  } // end fSecondColIsMCNu
227  } // end loop over MCTruths
228  } // If SecondCollection specified
229 
230  // put the collections in the event
231  evt.put(std::move(spillcol));
232  evt.put(std::move(truthcol));
233  evt.put(std::move(fluxcol));
234  evt.put(std::move(gtruthcol));
235  evt.put(std::move(dk2nucol));
236  evt.put(std::move(nuccol));
237  evt.put(std::move(assns));
238  evt.put(std::move(gassns));
239  }
240 
241  //--------------------------------------------------------------------
242  void MergeGenCollections::endSubRun(art::SubRun &sr)
243  {
244  std::unique_ptr< sumdata::POTSum > potcol(new sumdata::POTSum);
245 
246  if(fFirstColIsMCNu) {
248  sr.getByLabel(fFirstCollection, pots);
249 
250  if(pots.isValid()) *potcol = *pots;
251  }
252  else if(fSecondColIsMCNu) {
253  mf::LogWarning("MergeGenCollections")
254  << "Ignoring second-collection POTSum, using only first";
255  }
256  sr.put(std::move(potcol));
257  }
258 
260 } // end namespace merge
261 /////////////////////////////////////////////////////////////////////////
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
int CCNC() const
Definition: MCNeutrino.h:148
double QSqr() const
Definition: MCNeutrino.h:157
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
int HitQuark() const
Definition: MCNeutrino.h:153
Merge generation collections MCTruth/MCFlux/GTruth/Dk2Nu/NuChoice.
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
int fSecondaryModeOffset
Merge Dk2Nu/NuChoice objects?
int pots
const char * p
Definition: xmltok.h:285
simb::Origin_t Origin() const
Definition: MCTruth.h:73
int HitNuc() const
Definition: MCNeutrino.h:152
OStream cerr
Definition: OStream.cxx:7
std::string fSecondCollection
MCTruth label for main collection.
bool fFirstColIsMCNu
Process label for the second collection.
int NParticles() const
Definition: MCTruth.h:74
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
DEFINE_ART_MODULE(TestTMapFile)
bsim::Dk2Nu * dk2nu
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
Loaders::FluxType flux
object containing MC flux information
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
bool isValid() const
Definition: Handle.h:189
int InteractionType() const
Definition: MCNeutrino.h:150
double W() const
Definition: MCNeutrino.h:154
double Y() const
Definition: MCNeutrino.h:156
int evt
caf::StandardRecord * sr
double X() const
Definition: MCNeutrino.h:155
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
int Target() const
Definition: MCNeutrino.h:151
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
ProductID put(std::unique_ptr< PROD > &&)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Event generator information.
Definition: MCTruth.h:32
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149
size_t size() const
enum BeamMode string