MergeCollections_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 // \file MergeCollections.cxx
3 // \brief
4 // \authors janow@fnal.gov
5 ///////////////////////////////////////////////////////////////////////////
6 //#include "NovaSimMixer/MergeCollections.h"
7 
8 #include <math.h>
9 #include <string>
10 
14 #include "fhiclcpp/ParameterSet.h"
18 
19 
20 #include "DAQChannelMap/DAQChannelMap.h"
24 
25 
26 #include "Utilities/AssociationUtil.h"
27 #include "Simulation/FLSHit.h"
28 #include "Simulation/FLSHitList.h"
29 #include "Simulation/Particle.h"
30 
31 
32 ///group hits according to time and space
33 namespace MergeCollections {
34 
36  public:
37  explicit MergeCollections(fhicl::ParameterSet const &pset);
38  virtual ~MergeCollections();
39 
40  void produce(art::Event& evt);
41 
42  protected:
43  // internal methods
44 
45  // configuration settings
46  std::string fFirstGenCollection; /// Label on the main collection
47  std::string fFirstG4Collection; /// Label on the main collection
48  std::string fSecondCollection; /// Label on the main collection
50 
51 
52  };
53 
54 
55  //----------------------------------------------------------------------
57  EDProducer(pset),
58  fFirstGenCollection (pset.get< std::string >("FirstGenCollection")),
59  fFirstG4Collection (pset.get< std::string >("FirstG4Collection")),
60  fSecondCollection (pset.get< std::string >("SecondCollection")),
61  fIsFirstCollectionCosmic (pset.get< int >("IsFirstCollectionCosmic"))
62 
63  { produces< std::vector<simb::MCTruth> >();
64  produces< std::vector<simb::MCFlux> >();
65  produces< std::vector<sim::FLSHitList > >();
66  produces< std::vector<sim::Particle> >();
67  produces< std::vector<simb::GTruth> >();
68  produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
69  produces< art::Assns<simb::MCTruth, simb::GTruth> >();
70  produces< art::Assns<sim::Particle, simb::MCTruth> >();
71 
72 
73  }
74 
75  //----------------------------------------------------------------------
76  MergeCollections::~MergeCollections()
77  {
78  }
79 
80  //----------------------------------------------------------------------
81  void MergeCollections::produce(art::Event& evt)
82  {
83 
84 
85  std::unique_ptr< std::vector<simb::MCTruth> > truthcol (new std::vector<simb::MCTruth>);
86  std::unique_ptr< std::vector<simb::MCFlux> > fluxcol (new std::vector<simb::MCFlux >);
87  std::unique_ptr< std::vector<simb::GTruth> > gtruthcol(new std::vector<simb::GTruth>);
88  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> > assns (new art::Assns<simb::MCTruth, simb::MCFlux>);
89  std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > gassns(new art::Assns<simb::MCTruth, simb::GTruth>);
90  std::unique_ptr< art::Assns<sim::Particle, simb::MCTruth> > passns(new art::Assns<sim::Particle, simb::MCTruth>);
91 
92 
93  // Define the FLSHitList and Particle vectors.
94  std::unique_ptr<std::vector<sim::FLSHitList> > flshlcol(new std::vector<sim::FLSHitList> );
95  std::unique_ptr<std::vector<sim::Particle> > pcol (new std::vector<sim::Particle> );
96 
97 
98  /// Current event number
99 // fCurrentEvent = evt.id().event();
100 
101  // Pull the MC generator information out of the event
103  evt.getByLabel(fFirstGenCollection, mclist);
104 
107 
108  if(!fIsFirstCollectionCosmic) {
109  evt.getByLabel(fFirstGenCollection, fllist);
110  evt.getByLabel(fFirstGenCollection, glist);
111 
112  }
113 
115  evt.getByLabel(fFirstG4Collection, flslist);
116 
117 
119  evt.getByLabel(fFirstG4Collection, partlist);
120 
121 
122  int IdOffset=0;
123 
124  for(unsigned int i = 0; i < partlist->size(); ++i){
125 
126  sim::Particle part ((*partlist)[i]);
127  IdOffset=part.TrackId();
128  // std::cout << " Primary Collection ID " << IdOffset << std::endl;
129  pcol->push_back(part);
130 
131  }
132 
133 
134  /// Loop over neutrino interactions
135 
136  // std::cout << " primary mc list size " << mclist->size() << std::endl;
137 
138  for(unsigned int i = 0; i < mclist->size(); ++i){
139 
140  simb::MCTruth truth ((*mclist)[i]);
141  truthcol->push_back(truth);
142 
143  if(! fIsFirstCollectionCosmic){
144  simb::MCFlux flux ((*fllist)[i]);
145  fluxcol ->push_back(flux);
146 
147  simb::GTruth gtruth ((*glist)[i]);
148  gtruthcol->push_back(gtruth);
149 
150  util::CreateAssn(evt, *truthcol, *fluxcol, *assns, fluxcol ->size()-1, fluxcol ->size());
151  util::CreateAssn(evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size());
152 
153  }
154 
155  // IdOffset=truth.NParticles();
156 
157  // std::cout<<"First collection:"<< flslist->size() << std::endl;
158 
159  if(i<flslist->size()){
160  sim::FLSHitList mcfls ((*flslist)[i]);
161  flshlcol->push_back(mcfls);
162  }
163  }
164 
165  ///Second collection
166  // std::cout<<"Id Offset = "<<IdOffset<<std::endl;
167 
169  evt.getByLabel(fSecondCollection, mc2list);
170 
172  evt.getByLabel(fSecondCollection, g2list);
173 
175  evt.getByLabel(fSecondCollection, fl2list);
176 
178  evt.getByLabel(fSecondCollection, fls2list);
179 
180 
182  evt.getByLabel(fSecondCollection, part2list);
183 
184  int LastPartId;
185  int *Offset = new int[1000];
186  int *OffsetParticles = new int[1000];
187  int OffsetIndex;
188 
189  LastPartId=0;
190  for(int i=0;i<10;i++){
191  Offset[i]=0;
192  OffsetParticles[i]=0;
193  }
194  Offset[0]=IdOffset;
195  OffsetParticles[0]=partlist->size();
196  OffsetIndex=1;
197 
198  for(unsigned int i = 0; i < part2list->size(); ++i){
199 
200  OffsetParticles[OffsetIndex]++;
201 
202  sim::Particle part ((*part2list)[i]);
203  // std::cout<<"partID="<<part.TrackId()<<std::endl;
204  if(LastPartId > part.TrackId() ){
205  Offset[OffsetIndex]= Offset[OffsetIndex-1] + LastPartId;
206  OffsetIndex++ ;
207 
208  // std::cout<<"osfet for i="<<i<<" is "<<Offset[OffsetIndex-1]<<std::endl;
209  // std::cin.get();
210  }
211 
212  LastPartId=part.TrackId();
213 
214  int partmother =0;
215  if (part.Mother()!=0)partmother = part.Mother() + Offset[OffsetIndex-1];
216 
217  sim::Particle gpart(part.TrackId()+ Offset[OffsetIndex-1],
218  part.PdgCode(),
219  part.Process(),
220  partmother,
221  part.Mass(),
222  part.StatusCode());
223  // std::cout <<part.TrackId()<< "\t" << gpart.TrackId()<<"\t"<< Offset[OffsetIndex-1] << " mother " << part.Mother() << " mother offset " << partmother << std::endl;//std::cin.get();
224 
225  gpart.SetGvtx(part.Gvx(), part.Gvy(),part.Gvz(),part.Gvt());
226  gpart.SetPolarization(part.Polarization());
227  gpart.SetWeight(part.Weight());
228  gpart.SetRescatter(part.Rescatter());
229  for(int d=0; d<part.NumberDaughters(); d++){
230  gpart.AddDaughter(part.Daughter(d)+ Offset[OffsetIndex-1]);
231  }
232  for(unsigned int t=0; t<part.NumberTrajectoryPoints(); t++){
233  //std::cout<<part.Vx(t)<<"\t"<<part.EndX()<<std::endl;
234  gpart.AddTrajectoryPoint(part.Position(t),part.Momentum(t));
235  }
236  pcol->push_back(gpart);
237  // std::cout<< "mc particle " << part.TrackId() << " " << Offset[OffsetIndex-1] << std::endl;
238  }
239 
240  std::cout<<"mc1 size="<<mclist->size()<<std::endl;
241  std::cout<<"mc2 size="<<mc2list->size()<<std::endl;
242 
243  for(unsigned int i = 0; i < mc2list->size(); ++i){
244 
245  simb::MCTruth truth ((*mc2list)[i]);
246  simb::MCFlux flux ((*fl2list)[i]);
247  simb::GTruth gtruth ((*g2list) [i]);
248 
249  simb::MCTruth tmc(truth);
250 
251 
252  /// Offset ids for FLSHits
253 
254  // std::cout<<"size="<<fls2list->size()<< " mclist " << i << std::endl;
255  if(i<fls2list->size()){
256 
257  sim::FLSHitList mcfls ((*fls2list)[i]);
258  sim::FLSHitList tmcfls;
259  sim::FLSHit fHit;
260 
261  const std::vector<sim::FLSHit>& hits = mcfls.fHits;
262  for (unsigned int j=0; j<hits.size(); ++j)
263  {
264  fHit.Clear();
265  fHit=hits[j];
266  fHit.SetTrackId(hits[j].GetTrackID() + Offset[i]);
267  tmcfls.fHits.push_back(fHit);
268 
269  // std::cout << " fls ID " << hits[j].GetTrackID() << " " << Offset[i] << " " << hits[j].GetTrackID()+Offset[i] << std::endl;
270  }
271  flshlcol->push_back(tmcfls);
272 
273  }
274  /// End of the FLS offset
275 
276  truthcol->push_back(tmc);
277  gtruthcol->push_back(gtruth);
278  fluxcol ->push_back(flux);
279 
280  util::CreateAssn(evt, *truthcol, *fluxcol, *assns, fluxcol->size()-1, fluxcol->size());
281  util::CreateAssn(evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size());
282 
283 
284  // std::cout<<"end of i mcsize"<<std::endl;
285 
286 
287  }
288 
289  // Define the FLSHitList and Particle vectors.
290 
291 
292  std::unique_ptr< std::vector<simb::MCTruth> > temptruelist (new std::vector<simb::MCTruth>);
293  std::unique_ptr< std::vector<sim::Particle> > temppartlist (new std::vector<sim::Particle>);
294 
295  bool NewMCList = false;
296  int MCListIndex = -1;
297  for(unsigned int j = 0; j < pcol->size(); ++j){
298  sim::Particle part ((*pcol)[j]);
299  if(part.Mother()!=0 )NewMCList=false;
300  if(part.Mother()==0 && ! NewMCList){
301  MCListIndex++;
302  // std::cout << " pushing truthcol index " << MCListIndex << " onto temptrue " << std::endl;
303 
304  simb::MCTruth truth ((*truthcol)[MCListIndex]);
305  temptruelist->push_back(truth);
306  NewMCList=true;
307  }
308  // std::cout << j << " " << temptruelist->size()-1 << std::endl;
309 
310  temppartlist->push_back(part);
311  util::CreateAssn(evt, *temppartlist, *temptruelist, *passns, temptruelist->size()-1,temptruelist->size());
312 
313 
314  }
315 
316 
317  // put the collections in the event
318  evt.put(std::move(temptruelist));
319  evt.put(std::move(fluxcol));
320  evt.put(std::move(gtruthcol));
321  evt.put(std::move(assns));
322  evt.put(std::move(gassns));
323  evt.put(std::move(flshlcol));
324  evt.put(std::move(temppartlist));
325  evt.put(std::move(passns));
326  // int dum;
327  // std::cin >> dum;
328 
329  delete [] Offset ;
330  delete [] OffsetParticles;
331 
332  }
333 
334 
335 
337 } // end namespace merge
338 /////////////////////////////////////////////////////////////////////////
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
const TVector3 & Polarization() const
Definition: MCParticle.h:213
void SetTrackId(const int trackid)
Definition: FLSHit.h:101
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
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 PdgCode() const
Definition: MCParticle.h:211
std::vector< sim::FLSHit > fHits
Definition: FLSHitList.h:21
double Gvz() const
Definition: MCParticle.h:247
double Weight() const
Definition: MCParticle.h:253
std::string fFirstG4Collection
Label on the main collection.
int Mother() const
Definition: MCParticle.h:212
double Mass() const
Definition: MCParticle.h:238
int StatusCode() const
Definition: MCParticle.h:210
double Gvx() const
Definition: MCParticle.h:245
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
DEFINE_ART_MODULE(TestTMapFile)
double Gvy() const
Definition: MCParticle.h:246
int fIsFirstCollectionCosmic
Label on the main collection.
std::string Process() const
Definition: MCParticle.h:214
Loaders::FluxType flux
int NumberDaughters() const
Definition: MCParticle.h:216
object containing MC flux information
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
group hits according to time and space
TString part[npart]
Definition: Style.C:32
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void hits()
Definition: readHits.C:15
int evt
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
std::string fSecondCollection
Label on the main collection.
OStream cout
Definition: OStream.cxx:6
size_t Offset(bool allow_default)
Value passed to –offset, or 0 if not specified.
Definition: Utilities.cxx:395
void SetGvtx(double *v)
Definition: MCParticle.cxx:120
A vector of FLSHit from single neutrino interaction.
Definition: FLSHitList.h:13
double Gvt() const
Definition: MCParticle.h:248
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void Clear()
Clear the FLS hit.
Definition: FLSHit.cxx:39
int Rescatter() const
Definition: MCParticle.h:251
Event generator information.
Definition: MCTruth.h:32
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
enum BeamMode string