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  fFirstGenCollection (pset.get< std::string >("FirstGenCollection")),
58  fFirstG4Collection (pset.get< std::string >("FirstG4Collection")),
59  fSecondCollection (pset.get< std::string >("SecondCollection")),
60  fIsFirstCollectionCosmic (pset.get< int >("IsFirstCollectionCosmic"))
61 
62  { produces< std::vector<simb::MCTruth> >();
63  produces< std::vector<simb::MCFlux> >();
64  produces< std::vector<sim::FLSHitList > >();
65  produces< std::vector<sim::Particle> >();
66  produces< std::vector<simb::GTruth> >();
67  produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
68  produces< art::Assns<simb::MCTruth, simb::GTruth> >();
69  produces< art::Assns<sim::Particle, simb::MCTruth> >();
70 
71 
72  }
73 
74  //----------------------------------------------------------------------
75  MergeCollections::~MergeCollections()
76  {
77  }
78 
79  //----------------------------------------------------------------------
80  void MergeCollections::produce(art::Event& evt)
81  {
82 
83 
84  std::unique_ptr< std::vector<simb::MCTruth> > truthcol (new std::vector<simb::MCTruth>);
85  std::unique_ptr< std::vector<simb::MCFlux> > fluxcol (new std::vector<simb::MCFlux >);
86  std::unique_ptr< std::vector<simb::GTruth> > gtruthcol(new std::vector<simb::GTruth>);
87  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> > assns (new art::Assns<simb::MCTruth, simb::MCFlux>);
88  std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > gassns(new art::Assns<simb::MCTruth, simb::GTruth>);
89  std::unique_ptr< art::Assns<sim::Particle, simb::MCTruth> > passns(new art::Assns<sim::Particle, simb::MCTruth>);
90 
91 
92  // Define the FLSHitList and Particle vectors.
93  std::unique_ptr<std::vector<sim::FLSHitList> > flshlcol(new std::vector<sim::FLSHitList> );
94  std::unique_ptr<std::vector<sim::Particle> > pcol (new std::vector<sim::Particle> );
95 
96 
97  /// Current event number
98 // fCurrentEvent = evt.id().event();
99 
100  // Pull the MC generator information out of the event
102  evt.getByLabel(fFirstGenCollection, mclist);
103 
106 
107  if(!fIsFirstCollectionCosmic) {
108  evt.getByLabel(fFirstGenCollection, fllist);
109  evt.getByLabel(fFirstGenCollection, glist);
110 
111  }
112 
114  evt.getByLabel(fFirstG4Collection, flslist);
115 
116 
118  evt.getByLabel(fFirstG4Collection, partlist);
119 
120 
121  int IdOffset=0;
122 
123  for(unsigned int i = 0; i < partlist->size(); ++i){
124 
125  sim::Particle part ((*partlist)[i]);
126  IdOffset=part.TrackId();
127  // std::cout << " Primary Collection ID " << IdOffset << std::endl;
128  pcol->push_back(part);
129 
130  }
131 
132 
133  /// Loop over neutrino interactions
134 
135  // std::cout << " primary mc list size " << mclist->size() << std::endl;
136 
137  for(unsigned int i = 0; i < mclist->size(); ++i){
138 
139  simb::MCTruth truth ((*mclist)[i]);
140  truthcol->push_back(truth);
141 
142  if(! fIsFirstCollectionCosmic){
143  simb::MCFlux flux ((*fllist)[i]);
144  fluxcol ->push_back(flux);
145 
146  simb::GTruth gtruth ((*glist)[i]);
147  gtruthcol->push_back(gtruth);
148 
149  util::CreateAssn(*this, evt, *truthcol, *fluxcol, *assns, fluxcol ->size()-1, fluxcol ->size());
150  util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size());
151 
152  }
153 
154  // IdOffset=truth.NParticles();
155 
156  // std::cout<<"First collection:"<< flslist->size() << std::endl;
157 
158  if(i<flslist->size()){
159  sim::FLSHitList mcfls ((*flslist)[i]);
160  flshlcol->push_back(mcfls);
161  }
162  }
163 
164  ///Second collection
165  // std::cout<<"Id Offset = "<<IdOffset<<std::endl;
166 
168  evt.getByLabel(fSecondCollection, mc2list);
169 
171  evt.getByLabel(fSecondCollection, g2list);
172 
174  evt.getByLabel(fSecondCollection, fl2list);
175 
177  evt.getByLabel(fSecondCollection, fls2list);
178 
179 
181  evt.getByLabel(fSecondCollection, part2list);
182 
183  int LastPartId;
184  int *Offset = new int[1000];
185  int *OffsetParticles = new int[1000];
186  int OffsetIndex;
187 
188  LastPartId=0;
189  for(int i=0;i<10;i++){
190  Offset[i]=0;
191  OffsetParticles[i]=0;
192  }
193  Offset[0]=IdOffset;
194  OffsetParticles[0]=partlist->size();
195  OffsetIndex=1;
196 
197  for(unsigned int i = 0; i < part2list->size(); ++i){
198 
199  OffsetParticles[OffsetIndex]++;
200 
201  sim::Particle part ((*part2list)[i]);
202  // std::cout<<"partID="<<part.TrackId()<<std::endl;
203  if(LastPartId > part.TrackId() ){
204  Offset[OffsetIndex]= Offset[OffsetIndex-1] + LastPartId;
205  OffsetIndex++ ;
206 
207  // std::cout<<"osfet for i="<<i<<" is "<<Offset[OffsetIndex-1]<<std::endl;
208  // std::cin.get();
209  }
210 
211  LastPartId=part.TrackId();
212 
213  int partmother =0;
214  if (part.Mother()!=0)partmother = part.Mother() + Offset[OffsetIndex-1];
215 
216  sim::Particle gpart(part.TrackId()+ Offset[OffsetIndex-1],
217  part.PdgCode(),
218  part.Process(),
219  partmother,
220  part.Mass(),
221  part.StatusCode());
222  // std::cout <<part.TrackId()<< "\t" << gpart.TrackId()<<"\t"<< Offset[OffsetIndex-1] << " mother " << part.Mother() << " mother offset " << partmother << std::endl;//std::cin.get();
223 
224  gpart.SetGvtx(part.Gvx(), part.Gvy(),part.Gvz(),part.Gvt());
225  gpart.SetPolarization(part.Polarization());
226  gpart.SetWeight(part.Weight());
227  gpart.SetRescatter(part.Rescatter());
228  for(int d=0; d<part.NumberDaughters(); d++){
229  gpart.AddDaughter(part.Daughter(d)+ Offset[OffsetIndex-1]);
230  }
231  for(unsigned int t=0; t<part.NumberTrajectoryPoints(); t++){
232  //std::cout<<part.Vx(t)<<"\t"<<part.EndX()<<std::endl;
233  gpart.AddTrajectoryPoint(part.Position(t),part.Momentum(t));
234  }
235  pcol->push_back(gpart);
236  // std::cout<< "mc particle " << part.TrackId() << " " << Offset[OffsetIndex-1] << std::endl;
237  }
238 
239  std::cout<<"mc1 size="<<mclist->size()<<std::endl;
240  std::cout<<"mc2 size="<<mc2list->size()<<std::endl;
241 
242  for(unsigned int i = 0; i < mc2list->size(); ++i){
243 
244  simb::MCTruth truth ((*mc2list)[i]);
245  simb::MCFlux flux ((*fl2list)[i]);
246  simb::GTruth gtruth ((*g2list) [i]);
247 
248  simb::MCTruth tmc(truth);
249 
250 
251  /// Offset ids for FLSHits
252 
253  // std::cout<<"size="<<fls2list->size()<< " mclist " << i << std::endl;
254  if(i<fls2list->size()){
255 
256  sim::FLSHitList mcfls ((*fls2list)[i]);
257  sim::FLSHitList tmcfls;
258  sim::FLSHit fHit;
259 
260  const std::vector<sim::FLSHit>& hits = mcfls.fHits;
261  for (unsigned int j=0; j<hits.size(); ++j)
262  {
263  fHit.Clear();
264  fHit=hits[j];
265  fHit.SetTrackId(hits[j].GetTrackID() + Offset[i]);
266  tmcfls.fHits.push_back(fHit);
267 
268  // std::cout << " fls ID " << hits[j].GetTrackID() << " " << Offset[i] << " " << hits[j].GetTrackID()+Offset[i] << std::endl;
269  }
270  flshlcol->push_back(tmcfls);
271 
272  }
273  /// End of the FLS offset
274 
275  truthcol->push_back(tmc);
276  gtruthcol->push_back(gtruth);
277  fluxcol ->push_back(flux);
278 
279  util::CreateAssn(*this, evt, *truthcol, *fluxcol, *assns, fluxcol->size()-1, fluxcol->size());
280  util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size());
281 
282 
283  // std::cout<<"end of i mcsize"<<std::endl;
284 
285 
286  }
287 
288  // Define the FLSHitList and Particle vectors.
289 
290 
291  std::unique_ptr< std::vector<simb::MCTruth> > temptruelist (new std::vector<simb::MCTruth>);
292  std::unique_ptr< std::vector<sim::Particle> > temppartlist (new std::vector<sim::Particle>);
293 
294  bool NewMCList = false;
295  int MCListIndex = -1;
296  for(unsigned int j = 0; j < pcol->size(); ++j){
297  sim::Particle part ((*pcol)[j]);
298  if(part.Mother()!=0 )NewMCList=false;
299  if(part.Mother()==0 && ! NewMCList){
300  MCListIndex++;
301  // std::cout << " pushing truthcol index " << MCListIndex << " onto temptrue " << std::endl;
302 
303  simb::MCTruth truth ((*truthcol)[MCListIndex]);
304  temptruelist->push_back(truth);
305  NewMCList=true;
306  }
307  // std::cout << j << " " << temptruelist->size()-1 << std::endl;
308 
309  temppartlist->push_back(part);
310  util::CreateAssn(*this, evt, *temppartlist, *temptruelist, *passns, temptruelist->size()-1,temptruelist->size());
311 
312 
313  }
314 
315 
316  // put the collections in the event
317  evt.put(std::move(temptruelist));
318  evt.put(std::move(fluxcol));
319  evt.put(std::move(gtruthcol));
320  evt.put(std::move(assns));
321  evt.put(std::move(gassns));
322  evt.put(std::move(flshlcol));
323  evt.put(std::move(temppartlist));
324  evt.put(std::move(passns));
325  // int dum;
326  // std::cin >> dum;
327 
328  delete [] Offset ;
329  delete [] OffsetParticles;
330 
331  }
332 
333 
334 
336 } // end namespace merge
337 /////////////////////////////////////////////////////////////////////////
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
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
TString part[npart]
Definition: Style.C:32
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
void SetGvtx(double *v)
Definition: MCParticle.cxx:120
A vector of FLSHit from single neutrino interaction.
Definition: FLSHitList.h:13
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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
size_t size() const
enum BeamMode string