DAQQuickClustering_module.cc
Go to the documentation of this file.
2 
3 void ClusterHitsInTime::DoIt(std::vector<recoHit> cHitVector)
4 {
5  fHitVector.clear();
6  fHitVector = cHitVector;
7 
8  //ORDER IN THE HIT VECTOR BY TIME.
9  std::sort(fHitVector.begin(), fHitVector.end(),
10  [](const recoHit& lhs, const recoHit& rhs){return lhs<rhs;});
11 
12  for(unsigned int i = 0; i < fHitVector.size()-1; i++)
13  {
14  //std::cout << "fHitVector.at(" << i << ").getHitTime() " << fHitVector.at(i).getHitTime()
15  // << std::endl;
16  std::vector<recoHit> vec_TempHits;
17  if(std::abs(fHitVector.at(i).getHitTime()-fHitVector.at(i+1).getHitTime())<=fTimeWindow)
18  {
19  int timeCount = 1;
20 
21  vec_TempHits.push_back(fHitVector.at(i));
22  vec_TempHits.push_back(fHitVector.at(i+1));
23 
24  while((i+timeCount+1)<fHitVector.size() &&
25  std::abs(fHitVector.at(i+timeCount).getHitTime()-fHitVector.at(i+timeCount+1).getHitTime()) <= fTimeWindow)
26  {
27  vec_TempHits.push_back(fHitVector.at(i + timeCount + 1));
28  timeCount++;
29  }
30 
31  i = i + timeCount;
32  cluster temp(fHitVector.at(0).getEvent(), vec_TempHits);
33  fVecClusters.push_back(temp);
34  }
35  }
36 
37  fNClusters = fVecClusters.size();
38 }
39 
40 //......................................................
41 cluster::cluster(int cEvent, std::vector<recoHit> cHitVector)
42 {
43 
44  int type(0);
45  std::set<int> channels;
46  fHitVector.clear();
47  fStartChan = 50000;
48  fEndChan = 0;
49  fFirstHitTime = 50000;
50  fLastHitTime = 0;
51  fEvent = cEvent;
52  fHitVector = cHitVector;
53  fNHits = fHitVector.size();
54  for(int i = 0; i < fNHits; i++)
55  {
56  fHitSADC += fHitVector.at(i).getHitSADC();
57  int chan = fHitVector.at(i).getHitChan();
58  float time = fHitVector.at(i).getHitTime();
59 
60  channels.insert(chan);
61 
62  if(chan < fStartChan)
63  fStartChan = chan;
64 
65  if(chan > fEndChan)
66  fEndChan = chan;
67 
68  if(time < fFirstHitTime)
69  fFirstHitTime = time;
70 
71  if(time > fLastHitTime)
72  fLastHitTime = time;
73 
74  if(fHitVector.at(i).getGenType()==1)
75  type++;
76  }
77 
78  fNChan = channels.size();
79 
80  //CALL THE CLUSTER MARLEY IF THERE ARE MORE THAN TWO MARLEY HITS IN IT.
81  if(type>=2)
82  fType = 1;
83  else
84  fType = 0;
85 
86  fChanWidth = fEndChan - fStartChan;
87  fTimeWidth = fLastHitTime - fFirstHitTime;
88 
89  //ORDER IN THE HIT VECTOR BY TIME.
90  std::sort(fHitVector.begin(), fHitVector.end(),
91  [](const recoHit& lhs, const recoHit& rhs){return lhs.getHitTime() < rhs.getHitTime();});
92 
93 }
94 
95 
96 //......................................................
98 
99  this->reconfigure(p);
100 
101 }
102 
103 
104 //......................................................
106 {
107 
108  fRawDigitLabel = p.get<std::string>("RawDigitLabel");
109  fHitLabel = p.get<std::string>("HitLabel" );
110 
111  fGEANTLabel = p.get<std::string>("GEANT4Label" );
112  fMARLLabel = p.get<std::string>("MARLEYLabel" );
113  fAPALabel = p.get<std::string>("APALabel" );
114  fCPALabel = p.get<std::string>("CPALabel" );
115  fAr39Label = p.get<std::string>("Argon39Label" );
116  fNeutLabel = p.get<std::string>("NeutronLabel" );
117  fKrypLabel = p.get<std::string>("KryptonLabel" );
118  fPlonLabel = p.get<std::string>("PoloniumLabel");
119  fRdonLabel = p.get<std::string>("RadonLabel" );
120  fAr42Label = p.get<std::string>("Argon42Label" );
121 
122 
123  cut_AdjChanTolerance = p.get<std::vector<int>> ("AdjChanTolerance");
124  cut_HitsInWindow = p.get<std::vector<int>> ("HitsInWindow");
125  cut_MinChannels = p.get<std::vector<int>> ("MinChannels");
126  cut_MinChanWidth = p.get<std::vector<int>> ("MinChanWidth");
127  cut_TimeWindowSize = p.get<std::vector<float>>("TimeWindowSize");
128  cut_TotalADC = p.get<std::vector<float>>("TotalADC");
129 
130  NConfigs = cut_AdjChanTolerance.size();
131  NCuts = 6;
132 
133  detectorScaling = p.get<double>("detectorScaling");
134 
135 } // Reconfigure
136 
137 
138 //......................................................
140 {
141  MarlParts.clear(); APAParts .clear(); CPAParts .clear(); Ar39Parts.clear();
142  NeutParts.clear(); KrypParts.clear(); PlonParts.clear(); RdonParts.clear();
143  Ar42Parts.clear();
144  trkIDToPType.clear();
145  Run = SubRun = Event = -1;
146 
149  TotGen_Ar42 = 0;
150 
151  out_HitView.clear();
152  out_GenType.clear();
153  out_HitChan.clear();
154  out_HitTime.clear();
155  out_HitSADC.clear();
156  out_HitRMS .clear();
157 
158  NTotHits = NColHits = NIndHits = 0;
160  for (int hh=0; hh<nMaxHits; ++hh)
161  {
162  HitView[hh] = HitSize[hh] = HitChan[hh] = GenType[hh] = 0;
163  HitTime[hh] = HitRMS [hh] = HitSADC[hh] = 0;
164  HitInt [hh] = HitPeak[hh] = HitTPC [hh] = 0;
165  NCorrespondingIDEs[hh] = 0;
166  Hit_X[hh] = 0;
167  Hit_Y[hh] = 0;
168  Hit_Z[hh] = 0;
169  Hit_Energy[hh] = 0;
170  Hit_NumElectrons[hh] = 0;
171  }
172 
173  MarlSample.clear();
174  MarlTime .clear();
175  MarlWeight.clear();
176  ENu = 0;
177  ENu_Lep = 0;
178  VertX = 0;
179  VertY = 0;
180  VertZ = 0;
181  VertexT = 0;
182  Px = 0;
183  Py = 0;
184  Pz = 0;
185  VertexChan = 0;
186  map_EventToMC.clear();
187 
188 
189 }
190 
191 
192 //......................................................
193 void DAQQuickClustering::clusterChannels(std::vector<recoHit> &vec_Hits,
194  std::vector<ClusterHitsInTime> &vec_ChannelCluster,
195  unsigned int const &config)
196 {
197  //HERE IT IS ASSUMED THAT THE HITS APPEAR SEQUENTIALLY BY CHANNEL.
198 
199  for(unsigned int i = 0; i < vec_Hits.size()-1; i++)
200  {
201  std::vector<recoHit> vec_TempHits;
202  double width = std::abs(vec_Hits.at(i ).getHitChan()-
203  vec_Hits.at(i+1).getHitChan());
204  //std::cout << cut_AdjChanTolerance.at(config) << std::endl;
205  if(width <= cut_AdjChanTolerance.at(config))
206  {
207  int channelCount = 1;
208 
209  vec_TempHits.push_back(vec_Hits.at(i ));
210  vec_TempHits.push_back(vec_Hits.at(i+1));
211 
212  while((i+channelCount+1)<vec_Hits.size() &&
213  std::abs(vec_Hits.at(i+channelCount).getHitChan()-vec_Hits.at(i+channelCount+1).getHitChan())<=cut_AdjChanTolerance.at(config))
214  {
215  vec_TempHits.push_back(vec_Hits.at(i + channelCount + 1));
216  channelCount++;
217  }
218 
219  i = i + channelCount;
221  temp.DoIt(vec_TempHits);
222 
223  vec_ChannelCluster.push_back(temp);
224  //std::cout << " vec_ChannelCluster.size() " << vec_ChannelCluster.size() << std::endl;
225 
226  }
227  }
228 
229  return;
230 }
231 
232 
233 //......................................................
234 void DAQQuickClustering::clusterCut(std::vector<cluster> &vec_Clusters, unsigned int const &config)
235 {
236  //REMEMBER WE NEED BOTH A MAXIMUM AND MINIMUM CHANNEL WIDTH DUE TO COSMICS.
237 
238  for(std::vector<cluster>::iterator it_Clusters=vec_Clusters.begin();
239  it_Clusters!=vec_Clusters.end();)
240  {
241  bool fail_MinChannels(false), fail_ChanWidth(false), fail_TotalADC(false);
242  if(it_Clusters->getNChan()<cut_MinChannels.at(config))
243  fail_MinChannels = true;
244 
245  if(it_Clusters->getChanWidth()<cut_MinChanWidth.at(config))
246  fail_ChanWidth = true;
247 
248  if(it_Clusters->getHitSADC()<cut_TotalADC.at(config))
249  fail_TotalADC = true;
250 
251  if(fail_MinChannels == true || fail_ChanWidth == true || fail_TotalADC == true)
252  it_Clusters = vec_Clusters.erase(it_Clusters);
253  else
254  it_Clusters++;
255  }
256  return;
257 }
258 
259 
260 //......................................................
261 void DAQQuickClustering::trigger(std::vector<cluster> &vec_Clusters, unsigned int const &config)
262 {
263 
264  for(unsigned int i = 0; i < vec_Clusters.size(); i++)
265  if(vec_Clusters.at(i).getNHits()>=cut_HitsInWindow.at(config))
266  vec_Clusters.at(i).setTriggerFlag(1);
267 
268  return;
269 }
270 
271 
272 //......................................................
274 {
276 
277  for(unsigned int i = 0; i < NConfigs; i++)
278  {
279  TGraph* g_config = tfs->make<TGraph>(NCuts);
280  g_config->SetName(Form("g_ConfigDefinitions%i",i));
281  g_config->SetPoint(0, 1, cut_AdjChanTolerance[i]);
282  g_config->SetPoint(1, 2, cut_HitsInWindow [i]);
283  g_config->SetPoint(2, 3, cut_MinChannels [i]);
284  g_config->SetPoint(3, 4, cut_MinChanWidth [i]);
285  g_config->SetPoint(4, 5, cut_TimeWindowSize [i]);
286  g_config->SetPoint(5, 6, cut_TotalADC [i]);
287  }
288 
289 }
290 
291 
292 //......................................................
294 {
295 
297 
298  makeConfigGraph();
299 
300  t_Output_unusedhits = tfs->make<TTree>("DAQQuickClustering_unusedhits",
301  "DAQQuickClustering_unusedhits");
302 
303  t_Output_unusedhits->Branch("Cluster", &out_Cluster, "Cluster/I");
304  t_Output_unusedhits->Branch("Event", &out_Event, "Event/I");
305  t_Output_unusedhits->Branch("Config", &out_Config, "Config/I");
306  t_Output_unusedhits->Branch("StartChan", &out_StartChan, "StartChan/I");
307  t_Output_unusedhits->Branch("EndChan", &out_EndChan, "EndChan/I");
308  t_Output_unusedhits->Branch("ChanWidth", &out_ChanWidth, "ChanWidth/I");
309  t_Output_unusedhits->Branch("NChan", &out_NChan, "NChan/I");
310  t_Output_unusedhits->Branch("Type", &out_Type, "Type/I");
311  t_Output_unusedhits->Branch("NHits", &out_NHits, "NHits/I");
312  t_Output_unusedhits->Branch("SumADC", &out_SumADC, "SumADC/F");
313  t_Output_unusedhits->Branch("FirstTimeHit", &out_FirstTimeHit,"FirstTimeHit/F");
314  t_Output_unusedhits->Branch("LastTimeHit", &out_LastTimeHit, "LastTimeHit/F");
315  t_Output_unusedhits->Branch("TimeWidth", &out_TimeWidth, "TimeWidth/F");
316  t_Output_unusedhits->Branch("ENu", &out_ENu, "ENu/D");
317  t_Output_unusedhits->Branch("ENu_Lep", &out_ENu_Lep, "ENu_Lep/D");
318  t_Output_unusedhits->Branch("MarlTime", &out_MarlTime, "MarlTime/D");
319  t_Output_unusedhits->Branch("HitView", &out_HitView);
320  t_Output_unusedhits->Branch("GenType", &out_GenType);
321  t_Output_unusedhits->Branch("HitChan", &out_HitChan);
322  t_Output_unusedhits->Branch("HitTime", &out_HitTime);
323  t_Output_unusedhits->Branch("HitSADC", &out_HitSADC);
324  t_Output_unusedhits->Branch("HitRMS", &out_HitRMS);
325 
326  t_Output_clusteredhits = tfs->make<TTree>("DAQQuickClustering_clusteredhits",
327  "DAQQuickClustering_clusteredhits");
328 
329  t_Output_clusteredhits->Branch("Cluster", &out_Cluster, "Cluster/I");
330  t_Output_clusteredhits->Branch("Event", &out_Event, "Event/I");
331  t_Output_clusteredhits->Branch("Config", &out_Config, "Config/I");
332  t_Output_clusteredhits->Branch("StartChan", &out_StartChan, "StartChan/I");
333  t_Output_clusteredhits->Branch("EndChan", &out_EndChan, "EndChan/I");
334  t_Output_clusteredhits->Branch("ChanWidth", &out_ChanWidth, "ChanWidth/I");
335  t_Output_clusteredhits->Branch("NChan", &out_NChan, "NChan/I");
336  t_Output_clusteredhits->Branch("Type", &out_Type, "Type/I");
337  t_Output_clusteredhits->Branch("NHits", &out_NHits, "NHits/I");
338  t_Output_clusteredhits->Branch("SumADC", &out_SumADC, "SumADC/F");
339  t_Output_clusteredhits->Branch("FirstTimeHit", &out_FirstTimeHit,"FirstTimeHit/F");
340  t_Output_clusteredhits->Branch("LastTimeHit", &out_LastTimeHit, "LastTimeHit/F");
341  t_Output_clusteredhits->Branch("TimeWidth", &out_TimeWidth, "TimeWidth/F");
342  t_Output_clusteredhits->Branch("ENu", &out_ENu, "ENu/D");
343  t_Output_clusteredhits->Branch("ENu_Lep", &out_ENu_Lep, "ENu_Lep/D");
344  t_Output_clusteredhits->Branch("MarlTime", &out_MarlTime, "MarlTime/D");
345  t_Output_clusteredhits->Branch("HitView", &out_HitView);
346  t_Output_clusteredhits->Branch("GenType", &out_GenType);
347  t_Output_clusteredhits->Branch("HitChan", &out_HitChan);
348  t_Output_clusteredhits->Branch("HitTime", &out_HitTime);
349  t_Output_clusteredhits->Branch("HitSADC", &out_HitSADC);
350  t_Output_clusteredhits->Branch("HitRMS", &out_HitRMS);
351  h_ENu_MC = tfs->make<TH1D>("h_ENu_MC","h_ENu_MC",35,0,50);
352  h_MarlTime_MC = tfs->make<TH1D>("h_MarlTime_MC","h_MarlTime_MC",100,-0.1,10.5);
353  h_TimeElapsed = tfs->make<TH1D>("h_TimeElapsed", "h_TimeElapsed", 50,0,0.5);
354  h_ENu_MC->Sumw2();
355  std::cout << "Finished beginJob" << std::endl;
356 }
357 
358 
360 {
361  PType ThisPType = kUnknown;
362  auto const& it=trkIDToPType.find(TrID);
363  if(it!=trkIDToPType.end()){
364  ThisPType=it->second;
365  }
366 
367  return ThisPType;
368 }
369 
370 
371 bool DAQQuickClustering::InMyMap( int TrID, std::map< int, simb::MCParticle> ParMap )
372 {
373 
375  ParIt = ParMap.find(TrID);
376 
377  if(ParIt != ParMap.end()) return true;
378  else return false;
379 
380 }
381 
382 
383 
385 {
386  std::cout << "Job ended." << std::endl;
387  std::cerr << "firstCatch " << firstCatch << std::endl;
388  std::cerr << "secondCatch " << secondCatch << std::endl;
389  std::cerr << "thirdCatch " << thirdCatch << std::endl;
390 }
391 
392 
393 void DAQQuickClustering::FillMyMaps(std::map<int, simb::MCParticle> &MyMap,
394  art::FindManyP<simb::MCParticle> Assn,
395  art::ValidHandle< std::vector<simb::MCTruth> > Hand)
396 {
397  for ( size_t L1=0; L1 < Hand->size(); ++L1 ) {
398  for ( size_t L2=0; L2 < Assn.at(L1).size(); ++L2 ) {
399  const simb::MCParticle ThisPar = (*Assn.at(L1).at(L2));
400  MyMap[ThisPar.TrackId()] = ThisPar;
401  }
402  }
403  return;
404 }
405 
406 //......................................................
408 {
409  ResetVariables();
410 
411  Run = evt.run();
412  SubRun = evt.subRun();
413  Event = evt.event();
414 
415  //GET INFORMATION ABOUT THE DETECTOR'S GEOMETRY.
416  auto const* geo = lar::providerFrom<geo::Geometry>();
417 
418 // GET THE RECO HITS.
419  auto reco_hits = evt.getValidHandle<std::vector<recob::Hit> >(fHitLabel);
420 
421  //LIFT OUT THE MARLEY PARTICLES.
422  auto MarlTrue = evt.getValidHandle<std::vector<simb::MCTruth> >(fMARLLabel);
423  art::FindManyP<simb::MCParticle> MarlAssn(MarlTrue,evt,fGEANTLabel);
424  FillMyMaps(MarlParts, MarlAssn, MarlTrue);
425  TotGen_Marl = MarlParts.size();
426 
427  //SUPERNOVA TRUTH.
428  art::FindManyP<sim::SupernovaTruth> SNTruth(MarlTrue, evt, fMARLLabel);
429 
430  double Px_(0), Py_(0), Pz_(0), Pnorm(1);
431  for(unsigned int i = 0; i < MarlTrue->size(); i++)
432  {
433  Nu_Type = MarlTrue->at(i).GetNeutrino().Nu().PdgCode();
434  ENu = MarlTrue->at(i).GetNeutrino().Nu().E();
435  Mode = MarlTrue->at(i).GetNeutrino().Mode();
436  CCNC = MarlTrue->at(i).GetNeutrino().CCNC();
437  Target = MarlTrue->at(i).GetNeutrino().Target();
438  HitNucleon = MarlTrue->at(i).GetNeutrino().HitNuc();
439  Nu_Lep_Type = MarlTrue->at(i).GetNeutrino().Lepton().PdgCode();
440  VertX = MarlTrue->at(i).GetNeutrino().Lepton().Vx();
441  VertY = MarlTrue->at(i).GetNeutrino().Lepton().Vy();
442  VertZ = MarlTrue->at(i).GetNeutrino().Lepton().Vz();
443  Px_ = MarlTrue->at(i).GetNeutrino().Lepton().Px();
444  Py_ = MarlTrue->at(i).GetNeutrino().Lepton().Py();
445  Pz_ = MarlTrue->at(i).GetNeutrino().Lepton().Pz();
446  ENu_Lep = MarlTrue->at(i).GetNeutrino().Lepton().E();
447  for (unsigned int j = 0; j < SNTruth.at(i).size(); j++)
448  {
449  const sim::SupernovaTruth ThisTr = (*SNTruth.at(i).at(j));
450  MarlTime .push_back(ThisTr.SupernovaTime);
451  MarlWeight.push_back(ThisTr.Weight);
452  MarlSample.push_back(ThisTr.SamplingMode);
453  }
454  }
455  Pnorm = std::sqrt(Px_*Px_+Py_*Py_+Pz_*Pz_);
456  Px = Px_/Pnorm;
457  Py = Py_/Pnorm;
458  Pz = Pz_/Pnorm;
459  bool caught = false;
460  double Vertex[3] = {VertX, VertY, VertZ};
462  geo::PlaneID Plane(geo->FindTPCAtPosition(Vertex),geo::kZ);
463  try
464  {
465  WireID = geo->NearestWireID(Vertex, Plane);
466  }
467  catch(...)
468  {
469  caught = true;
470  }
471 
472  if(caught)
473  VertexChan = -1;
474  else
475  VertexChan = geo->PlaneWireToChannel(WireID);
476 
477  //CM/MICROSECOND.
478  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
479  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
480  double drift_velocity = detProp.DriftVelocity(detProp.Efield(),detProp.Temperature());
481  //CM/TICK
482  drift_velocity = drift_velocity*0.5;
483  VertexT = VertX/drift_velocity;
484 
485  auto APATrue = evt.getValidHandle<std::vector<simb::MCTruth> >(fAPALabel);
486  art::FindManyP<simb::MCParticle> APAAssn(APATrue,evt,fGEANTLabel);
487  FillMyMaps( APAParts, APAAssn, APATrue );
488  TotGen_APA = APAParts.size();
489 
490  auto CPATrue = evt.getValidHandle<std::vector<simb::MCTruth> >(fCPALabel);
491  art::FindManyP<simb::MCParticle> CPAAssn(CPATrue,evt,fGEANTLabel);
492  FillMyMaps( CPAParts, CPAAssn, CPATrue );
493  TotGen_CPA = CPAParts.size();
494 
495  auto Ar39True = evt.getValidHandle<std::vector<simb::MCTruth> >(fAr39Label);
496  art::FindManyP<simb::MCParticle> Ar39Assn(Ar39True,evt,fGEANTLabel);
497  FillMyMaps( Ar39Parts, Ar39Assn, Ar39True );
498  TotGen_Ar39 = Ar39Parts.size();
499 
500  auto NeutTrue = evt.getValidHandle<std::vector<simb::MCTruth> >(fNeutLabel);
501  art::FindManyP<simb::MCParticle> NeutAssn(NeutTrue,evt,fGEANTLabel);
502  FillMyMaps( NeutParts, NeutAssn, NeutTrue );
503  TotGen_Neut = NeutParts.size();
504 
505  auto KrypTrue = evt.getValidHandle<std::vector<simb::MCTruth> >(fKrypLabel);
506  art::FindManyP<simb::MCParticle> KrypAssn(KrypTrue,evt,fGEANTLabel);
507  FillMyMaps( KrypParts, KrypAssn, KrypTrue );
508  TotGen_Kryp = KrypParts.size();
509 
510  auto PlonTrue = evt.getValidHandle<std::vector<simb::MCTruth> >(fPlonLabel);
511  art::FindManyP<simb::MCParticle> PlonAssn(PlonTrue,evt,fGEANTLabel);
512  FillMyMaps( PlonParts, PlonAssn, PlonTrue );
513  TotGen_Plon = PlonParts.size();
514 
515  auto RdonTrue = evt.getValidHandle<std::vector<simb::MCTruth> >(fRdonLabel);
516  art::FindManyP<simb::MCParticle> RdonAssn(RdonTrue,evt,fGEANTLabel);
517  FillMyMaps( RdonParts, RdonAssn, RdonTrue );
518  TotGen_Rdon = RdonParts.size();
519 
520  auto Ar42True = evt.getValidHandle<std::vector<simb::MCTruth> >(fAr42Label);
521  art::FindManyP<simb::MCParticle> Ar42Assn(Ar42True,evt,fGEANTLabel);
522  FillMyMaps( Ar42Parts, Ar42Assn, Ar42True );
523  TotGen_Ar42 = Ar42Parts.size();
524 
525  std::map<PType, std::map< int, simb::MCParticle >&> PTypeToMap{
526  { kMarl, MarlParts},
527  {kAPA, APAParts},
528  {kCPA, CPAParts},
529  {kAr39, Ar39Parts},
530  {kNeut, NeutParts},
531  {kKryp, KrypParts},
532  {kPlon, PlonParts},
533  {kRdon, RdonParts},
534  {kAr42, Ar42Parts}
535  };
536 
537  for(auto const& it : PTypeToMap){
538  const PType p=it.first;
539  auto const& m=it.second;
540  for(auto const& it2 : m){
541  trkIDToPType.insert(std::make_pair(it2.first, p));
542  }
543  }
544 
545  //std::cout << "THE EVENTS NUMBER IS: " << Event << std::endl;
546 
547  std::vector< recob::Hit > ColHits_Marl;
548  std::vector< recob::Hit > ColHits_CPA;
549  std::vector< recob::Hit > ColHits_APA;
550  std::vector< recob::Hit > ColHits_Ar39;
551  std::vector< recob::Hit > ColHits_Neut;
552  std::vector< recob::Hit > ColHits_Kryp;
553  std::vector< recob::Hit > ColHits_Plon;
554  std::vector< recob::Hit > ColHits_Rdon;
555  std::vector< recob::Hit > ColHits_Oth;
556  std::vector< recob::Hit > ColHits_Ar42;
557 
558  NTotHits = reco_hits->size();
559  int colHitCount(0);
560  int LoopHits = std::min( NTotHits, nMaxHits );
561  //if(LoopHits != NTotHits)
562  // std::cerr << "---- There are " << NTotHits
563  // << " hits in the event, but array is of size "
564  // << nMaxHits << ", so looping over first " << LoopHits
565  // << " hits." << std::endl;
566 
567  for(int hit = 0; hit < LoopHits; ++hit)
568  {
569  recob::Hit const& ThisHit = reco_hits->at(hit);
570  if(ThisHit.View() == geo::kU || ThisHit.View() == geo::kV) ++NIndHits;
571  else ++NColHits;
572  }
573 
574  for(int hit = 0; hit < LoopHits; ++hit)
575  {
576  recob::Hit const& ThisHit = reco_hits->at(hit);
577 
578  if (ThisHit.View() == 2) {
579  std::vector< sim::TrackIDE > ThisHitIDE;
580  //GETTING HOLD OF THE SIM::IDEs.
581  std::vector<const sim::IDE*> ThisSimIDE;
582  try
583  {
584  ThisHitIDE = bt_serv->HitToTrackIDEs(clockData, ThisHit );
585  }
586  catch(...)
587  {
588  //std::cout << "FIRST CATCH" << std::endl;
589  firstCatch++;
590  try
591  {
592  ThisSimIDE = bt_serv->HitToSimIDEs_Ps(clockData, ThisHit);
593  }
594  catch(...)
595  {
596  //std::cout << "SECOND CATCH" << std::endl;
597  secondCatch++;
598  continue;
599  }
600  continue;
601  }
602  try
603  {
604  ThisSimIDE = bt_serv->HitToSimIDEs_Ps(clockData, ThisHit);
605  }
606  catch(...)
607  {
608  //std::cout << "THIRD CATCH" << std::endl;
609  thirdCatch++;
610  continue;
611  }
612 
613  HitView[colHitCount] = ThisHit.View();
614  HitSize[colHitCount] = ThisHit.EndTick() - ThisHit.StartTick();
615  HitTPC [colHitCount] = ThisHit.WireID().TPC;
616  HitChan[colHitCount] = ThisHit.Channel();
617  HitTime[colHitCount] = ThisHit.PeakTime();
618  HitRMS [colHitCount] = ThisHit.RMS();
619  HitSADC[colHitCount] = ThisHit.SummedADC();
620  HitInt [colHitCount] = ThisHit.Integral();
621  HitPeak[colHitCount] = ThisHit.PeakAmplitude();
622  NCorrespondingIDEs[colHitCount] = ThisHitIDE.size();
623 
624  if(ThisHitIDE.size()==0)
626 
627  //WHICH PARTICLE CONTRIBUTED MOST TO THE HIT.
628  int MainTrID = -1;
629  double TopEFrac = -DBL_MAX;
630  for (size_t ideL=0; ideL < ThisHitIDE.size(); ++ideL)
631  {
632  if ( ThisHitIDE[ideL].energyFrac > TopEFrac )
633  {
634  TopEFrac = ThisHitIDE[ideL].energyFrac;
635  MainTrID = ThisHitIDE[ideL].trackID;
636  }
637  }
638 
639  PType ThisPType = WhichParType( MainTrID );
640  GenType[colHitCount] = ThisPType;
641 
642  if(MainTrID == -1)
643  {
644  Hit_X[colHitCount] = -1;
645  Hit_Y[colHitCount] = -1;
646  Hit_Z[colHitCount] = -1;
647  Hit_Energy[colHitCount] = -1;
648  Hit_NumElectrons[colHitCount] = -1;
649  }
650  else
651  {
652  for(unsigned int i = 0; i < ThisSimIDE.size(); i++)
653  {
654  if(ThisSimIDE.at(i)->trackID==MainTrID)
655  {
656  Hit_X[colHitCount] = ThisSimIDE.at(i)->x;
657  Hit_Y[colHitCount] = ThisSimIDE.at(i)->y;
658  Hit_Z[colHitCount] = ThisSimIDE.at(i)->z;
659  Hit_Energy[colHitCount] = ThisSimIDE.at(i)->energy;
660  Hit_NumElectrons[colHitCount] = ThisSimIDE.at(i)->numElectrons;
661  break;
662  }
663  }
664  }
665  switch (ThisPType){
666  case kUnknown: { ColHits_Oth .push_back( ThisHit ); break; }
667  case kMarl: { ColHits_Marl.push_back( ThisHit ); break; }
668  case kAPA: { ColHits_APA .push_back( ThisHit ); break; }
669  case kCPA: { ColHits_CPA .push_back( ThisHit ); break; }
670  case kAr39: { ColHits_Ar39.push_back( ThisHit ); break; }
671  case kNeut: { ColHits_Neut.push_back( ThisHit ); break; }
672  case kKryp: { ColHits_Kryp.push_back( ThisHit ); break; }
673  case kPlon: { ColHits_Plon.push_back( ThisHit ); break; }
674  case kRdon: { ColHits_Rdon.push_back( ThisHit ); break; }
675  case kAr42: { ColHits_Ar42.push_back( ThisHit ); break; }
676  default: break;
677  }
678 
679  colHitCount++;
680  }
681  }
682 
683  std::vector<int> vec_ClusterCount(NConfigs);
684 
685  h_ENu_MC ->Fill(1000*ENu);
686  h_MarlTime_MC->Fill(MarlTime.back());
687 
688  map_EventToMC[Event] = {ENu, ENu_Lep, MarlTime.back()};
689 
690  //MAKE RECOHIT OBJECTS EVENTWISE FROM THE TREE.
691  std::vector<recoHit> vec_Hits;
692  for(int j = 0; j < NColHits; j++)
693  {
694  recoHit hit(Event, HitView[j], GenType[j],
695  HitChan[j], HitTime[j], HitSADC[j], HitRMS[j]);
696  vec_Hits.push_back(hit);
697  }
698 
699  for(unsigned int j = 0; j < 1; j++)
700  {
701  std::vector<ClusterHitsInTime> vec_ChannelCluster;
702  std::vector<cluster> vec_Clusters;
703  TStopwatch *timeElapsed = new TStopwatch();
704  clusterChannels(vec_Hits, vec_ChannelCluster, j);
705 
706  for(unsigned int k = 0; k < vec_ChannelCluster.size(); k++)
707  {
708  std::vector<cluster> vec_Temp = vec_ChannelCluster.at(k).getClusterVector();
709  //std::cout << " vec_temp.size() " << vec_Temp.size() << std::endl;
710  for(unsigned int l = 0; l < vec_Temp.size(); l++)
711  vec_Clusters.push_back(vec_Temp.at(l));
712  }
713  clusterCut(vec_Clusters, j);
714  trigger (vec_Clusters, j);
715  h_TimeElapsed->Fill(timeElapsed->RealTime()*1000);
716  delete timeElapsed;
717  timeElapsed = NULL;
718  //std::cout << " vec_Cluster.size() " << vec_Clusters.size() << std::endl;
719 
720  //FILL THE OUTPUT TREE.
721  for(unsigned int k = 0; k < vec_Clusters.size(); k++)
722  {
723  if(vec_Clusters.at(k).getTriggerFlag()==1)
724  {
725  out_Config = j;
726  out_Cluster = vec_ClusterCount.at(j);
727  out_Event = vec_Clusters.at(k).getEvent();
728  out_StartChan = vec_Clusters.at(k).getStartChan();
729  out_EndChan = vec_Clusters.at(k).getEndChan();
730  out_ChanWidth = vec_Clusters.at(k).getChanWidth();
731  out_NChan = vec_Clusters.at(k).getNChan();
732  out_Type = vec_Clusters.at(k).getType();
733  out_NHits = vec_Clusters.at(k).getNHits();
734  out_SumADC = vec_Clusters.at(k).getHitSADC();
735  out_FirstTimeHit = vec_Clusters.at(k).getFirstTimeHit();
736  out_LastTimeHit = vec_Clusters.at(k).getLastTimeHit();
737  out_TimeWidth = vec_Clusters.at(k).getTimeWidth();
738  out_ENu = map_EventToMC[vec_Clusters.at(k).getEvent()].at(0);
739  out_ENu_Lep = map_EventToMC[vec_Clusters.at(k).getEvent()].at(1);
740  out_MarlTime = map_EventToMC[vec_Clusters.at(k).getEvent()].at(2);
741  for(unsigned int l = 0; l < vec_Clusters.at(k).getHits().size(); l++)
742  {
743  out_HitView.push_back(vec_Clusters.at(k).getHits().at(l).getHitView());
744  out_GenType.push_back(vec_Clusters.at(k).getHits().at(l).getGenType());
745  out_HitChan.push_back(vec_Clusters.at(k).getHits().at(l).getHitChan());
746  out_HitTime.push_back(vec_Clusters.at(k).getHits().at(l).getHitTime());
747  out_HitSADC.push_back(vec_Clusters.at(k).getHits().at(l).getHitSADC());
748  out_HitRMS .push_back(vec_Clusters.at(k).getHits().at(l).getHitRMS ());
749  }
750  t_Output_clusteredhits->Fill();
751  out_HitView.clear(); out_GenType.clear(); out_HitChan.clear();
752  out_HitTime.clear(); out_HitSADC.clear(); out_HitRMS.clear();
753 
754  vec_ClusterCount.at(j)++;
755  }
756  }
757  vec_ChannelCluster.clear();
758  vec_Clusters.clear();
759  }
760  vec_Hits.clear();
761 }
762 
intermediate_table::iterator iterator
std::vector< int > MarlSample
std::vector< int > cut_HitsInWindow
std::map< int, simb::MCParticle > NeutParts
std::vector< int > cut_MinChannels
EventNumber_t event() const
Definition: DataViewImpl.cc:85
std::map< int, PType > trkIDToPType
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
void FillMyMaps(std::map< int, simb::MCParticle > &MyMap, art::FindManyP< simb::MCParticle > Assn, art::ValidHandle< std::vector< simb::MCTruth > > Hand)
std::string string
Definition: nybbler.cc:12
Planes which measure V.
Definition: geo_types.h:130
geo::WireID WireID() const
Definition: Hit.h:233
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:220
std::vector< int > cut_MinChanWidth
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
std::vector< double > MarlTime
void analyze(art::Event const &evt) override
Planes which measure Z direction.
Definition: geo_types.h:132
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:232
std::vector< double > out_HitTime
Cluster finding and building.
std::map< int, std::vector< double > > map_EventToMC
static QStrList * l
Definition: config.cpp:1044
SupernovaSamplingMode_t SamplingMode
Method used to sample the supernova neutrino&#39;s energy and arrival time.
std::vector< double > MarlWeight
int TrackId() const
Definition: MCParticle.h:210
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
std::map< int, simb::MCParticle > CPAParts
std::map< int, simb::MCParticle > Ar42Parts
double Weight
Statistical weight for this neutrino vertex.
std::vector< double > out_HitRMS
T abs(T value)
Planes which measure U.
Definition: geo_types.h:129
std::vector< int > cut_AdjChanTolerance
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
static Config * config
Definition: config.cpp:1054
std::vector< float > cut_TimeWindowSize
std::map< int, simb::MCParticle > MarlParts
T get(std::string const &key) const
Definition: ParameterSet.h:271
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
p
Definition: test.py:223
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
std::map< int, simb::MCParticle > APAParts
RunNumber_t run() const
Definition: DataViewImpl.cc:71
void clusterChannels(std::vector< recoHit > &vec_Hits, std::vector< ClusterHitsInTime > &vec_ChannelCluster, unsigned int const &config)
raw::TDCtick_t StartTick() const
Initial tdc tick for hit.
Definition: Hit.h:216
void clusterCut(std::vector< cluster > &vec_Clusters, unsigned int const &config)
float getHitTime() const
Detector simulation of raw signals on wires.
raw::TDCtick_t EndTick() const
Final tdc tick for hit.
Definition: Hit.h:217
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
std::vector< float > cut_TotalADC
double SupernovaTime
Arrival time of the supernova neutrino (seconds)
std::vector< int > out_HitView
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
void DoIt(std::vector< recoHit > cHitVector)
static QCString type
Definition: declinfo.cpp:672
void trigger(std::vector< cluster > &vec_Clusters, unsigned int const &config)
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
Definition: Hit.h:223
void reconfigure(fhicl::ParameterSet const &p)
std::vector< cluster > fVecClusters
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
std::map< int, simb::MCParticle > PlonParts
const int nMaxHits
std::vector< const sim::IDE * > HitToSimIDEs_Ps(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
art::ServiceHandle< cheat::BackTrackerService > bt_serv
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
std::map< int, simb::MCParticle > Ar39Parts
TCEvent evt
Definition: DataStructs.cxx:7
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
std::map< int, simb::MCParticle > RdonParts
std::vector< double > out_HitSADC
recob::tracking::Plane Plane
Definition: TrackState.h:17
LArSoft geometry interface.
Definition: ChannelGeo.h:16
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:230
bool InMyMap(int TrID, std::map< int, simb::MCParticle > ParMap)
std::vector< int > out_GenType
QTextStream & endl(QTextStream &s)
std::map< int, simb::MCParticle > KrypParts
std::vector< int > out_HitChan
std::vector< recoHit > fHitVector
DAQQuickClustering(fhicl::ParameterSet const &p)