LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
DBclusterAna_module.cc
Go to the documentation of this file.
1 //
3 // DBSCANfinderAna class
4 //
5 // \author kinga.partyka@yale.edu
6 //
8 
9 //#ifndef DBCLUSTERANA_H
10 //#define DBCLUSTERANA_H
11 
12 #include <sstream>
13 #include <fstream>
14 #include <algorithm>
15 #include <functional>
16 #include <TH1.h>
17 #include <TH2.h>
18 #include <TH1F.h>
19 #include <TH2F.h>
20 #include "TDatabasePDG.h"
21 #include "TSystem.h"
22 #include <vector>
23 #include <string>
24 
25 //Framework includes
26 #include "art/Framework/Core/ModuleMacros.h"
27 #include "canvas/Persistency/Common/FindManyP.h"
28 #include "art/Framework/Principal/Event.h"
29 #include "fhiclcpp/ParameterSet.h"
30 #include "art/Framework/Principal/Handle.h"
31 #include "canvas/Persistency/Common/Ptr.h"
32 #include "canvas/Persistency/Common/PtrVector.h"
33 #include "art/Framework/Services/Registry/ServiceHandle.h"
34 #include "art/Framework/Services/Optional/TFileService.h"
35 #include "art/Framework/Services/Optional/TFileDirectory.h"
36 #include "messagefacility/MessageLogger/MessageLogger.h"
37 
38 //LArSoft includes
39 #include "larcore/Geometry/Geometry.h"
40 #include "larcorealg/Geometry/PlaneGeo.h"
41 #include "nusimdata/SimulationBase/MCTruth.h"
42 #include "nutools/ParticleNavigation/ParticleList.h"
43 #include "nutools/ParticleNavigation/EmEveIdCalculator.h"
44 #include "lardataobj/RecoBase/Hit.h"
45 #include "lardataobj/RecoBase/Cluster.h"
46 #include "lardataobj/RecoBase/Wire.h"
47 #include "lardataobj/RawData/RawDigit.h"
48 #include "larsim/MCCheater/BackTrackerService.h"
49 #include "larsim/MCCheater/ParticleInventoryService.h"
50 #include "lardata/Utilities/AssociationUtil.h"
51 
52 
53 #include "art/Framework/Core/EDAnalyzer.h"
54 
55 
56 
57 class TH1F;
58 class TH2F;
60 namespace cluster {
61 
62 
63  class DBclusterAna : public art::EDAnalyzer {
64 
65  public:
66 
67  explicit DBclusterAna(fhicl::ParameterSet const& pset);
68  virtual ~DBclusterAna();
69 
71  void analyze(const art::Event& evt);
72  void beginJob();
73 
74  private:
75  TH1F* fNoParticles_pdg;
76  TH1F* fNoParticles_trackid;
77  TH1F* fNoParticles_trackid_mother;
78  TH1F* fNoParticles_trackid_per_event;
79  TH1F* fNoParticles_pdg_per_event;
80  TH1F* fCl_for_Muon;
81  /* TH1F* fCl_for_Electron; */
82  /* TH1F* fCl_for_Positron; */
83  /* TH1F* fCl_for_Pion_111; */
84  /* TH1F* fCl_for_Pion_211; */
85  /* TH1F* fCl_for_Pion_m211;*/
86  /* TH1F* fCl_for_Proton; */
87  TH1F* fNoClustersInEvent;
88  TH1F* fPercentNoise;
89  TH1F* fno_of_clusters_per_track;
90  TH1F* fPercent_lost_muon_hits;
91  TH1F* fPercent_lost_electron_hits;
92  TH1F* fPercent_lost_positron_hits;
93  TH1F* fPercent_lost_111_hits;
94  TH1F* fPercent_lost_211_hits;
95  TH1F* fPercent_lost_m211_hits;
96  TH1F* fPercent_lost_2212_hits;
97  TH1F* fPercent_lost_2112_hits;
98 
99  TH1F* fPercent_lost_muon_energy;
100  TH1F* fPercent_lost_electron_energy;
101  TH1F* fPercent_lost_positron_energy;
102  TH1F* fPercent_lost_111_energy;
103  TH1F* fPercent_lost_211_energy;
104  TH1F* fPercent_lost_m211_energy;
105  TH1F* fPercent_lost_2212_energy;
106  TH1F* fPercent_lost_2112_energy;
107  TH1F* fEnergy;
108  TH2F* fbrian_in;
109  TH2F* fbrian_coll;
110 
111  std::string fDigitModuleLabel;
112  std::string fHitsModuleLabel;
113  std::string fLArG4ModuleLabel;
114  std::string fClusterFinderModuleLabel;
115  std::string fCalDataModuleLabel;
116  std::string fGenieGenModuleLabel;
117 
118 
119 
120  }; // class DBclusterAna
121 
122 }
123 
124 //#endif
125 
126 namespace cluster{
127 
128  //--------------------------------------------------------------------
129  DBclusterAna::DBclusterAna(fhicl::ParameterSet const& pset)
130  : EDAnalyzer(pset)
131  , fDigitModuleLabel (pset.get< std::string >("DigitModuleLabel") )
132  , fHitsModuleLabel (pset.get< std::string >("HitsModuleLabel") )
133  , fLArG4ModuleLabel (pset.get< std::string >("LArGeantModuleLabel") )
134  , fClusterFinderModuleLabel (pset.get< std::string >("ClusterFinderModuleLabel"))
135  , fCalDataModuleLabel (pset.get< std::string >("CalDataModuleLabel") )
136  , fGenieGenModuleLabel (pset.get< std::string >("GenieGenModuleLabel") )
137  {
138 
139 
140  }
141 
142  //------------------------------------------------------------------
143  DBclusterAna::~DBclusterAna()
144  {
145 
146  }
147 
148  void DBclusterAna::beginJob()
149  {
150 
151  // get access to the TFile service
152  art::ServiceHandle<art::TFileService> tfs;
153 
154  fNoParticles_pdg_per_event = tfs->make<TH1F>("fNoParticles_pdg_per_event","Average # of Particles per cluster for each event", 500,0 ,5);
155  fNoParticles_pdg=tfs->make<TH1F>("fNoParticles_pdg","Number of Particles in a Cluster for each cluster", 500,0 ,5);
156  fNoParticles_trackid=tfs->make<TH1F>("fNoParticles_trackid","Number of different TrackIDs in a Cluster", 300,0 ,30);
157 
158  fNoParticles_trackid_mother=tfs->make<TH1F>("fNoParticles_trackid_mother","Number of different TrackIDs in a Cluster(using mother)for each cluster", 300,0 ,30);
159 
160  fNoParticles_trackid_per_event=tfs->make<TH1F>("fNoParticles_trackid_per_event","Avg Number of different TrackIDs per Cluster per event", 300,0 ,30);
161  fCl_for_Muon=tfs->make<TH1F>("fCl_for_Muon","Number of Clusters for Muon per plane (pdg)", 1500,0 ,15);
162  // fCl_for_Electron=tfs->make<TH1F>("fCl_for_Electron","Number of Clusters for Electron (pdg)", 1500,0 ,15);
163  // fCl_for_Positron=tfs->make<TH1F>("fCl_for_Positron","Number of Clusters for Positron", 1500,0 ,15);
164  // fCl_for_Pion_111=tfs->make<TH1F>("fCl_for_Pion_111","Number of Clusters for Pion (111)", 1500,0 ,15);
165  // fCl_for_Pion_211=tfs->make<TH1F>("fCl_for_Pion_211","Number of Clusters for Pion (211)", 1500,0 ,15);
166  // fCl_for_Pion_m211=tfs->make<TH1F>("fCl_for_Pion_m211","Number of Clusters for Pion (-211)", 1500,0 ,15);
167  // fCl_for_Proton=tfs->make<TH1F>("fCl_for_Proton","Number of Clusters for Proton", 1500,0 ,15);
168 
169  fNoClustersInEvent=tfs->make<TH1F>("fNoClustersInEvent","Number of Clusters in an Event", 50,0 ,50);
170 
171  fPercentNoise=tfs->make<TH1F>("fPercentNoise","% of hits that were marked as Noise by DBSCAN",250,0 ,25);
172 
173  fno_of_clusters_per_track=tfs->make<TH1F>("fno_of_clusters_per_track","Number of Clusters per TrackID per plane", 1500,0 ,15);
174 
175  fPercent_lost_muon_hits=tfs->make<TH1F>("fPercent_lost_muon_hits","Number of muon hits excluded by dbscan in % (per Event)", 10000,0 ,100);
176  fPercent_lost_electron_hits=tfs->make<TH1F>("fPercent_lost_electron_hits","Number of electron hits excluded by dbscan in % (per Event)", 10000,0 ,100);
177  fPercent_lost_positron_hits=tfs->make<TH1F>("fPercent_lost_positron_hits","Number of positron hits excluded by dbscan in % (per Event)", 10000,0 ,100);
178  fPercent_lost_111_hits=tfs->make<TH1F>("fPercent_lost_111_hits","Number of pion(111) hits excluded by dbscan in % (per Event)", 10000,0 ,100);
179  fPercent_lost_211_hits=tfs->make<TH1F>("fPercent_lost_211_hits","Number of pion(211) hits excluded by dbscan in % (per Event)", 10000,0 ,100);
180  fPercent_lost_m211_hits=tfs->make<TH1F>("fPercent_lost_m211_hits","Number of pion(-211) hits excluded by dbscan in % (per Event)", 10000,0 ,100);
181  fPercent_lost_2212_hits=tfs->make<TH1F>("fPercent_lost_2212_hits","Number of proton hits excluded by dbscan in % (per Event)", 10000,0 ,100);
182  fPercent_lost_2112_hits=tfs->make<TH1F>("fPercent_lost_2112_hits","Number of neutron hits excluded by dbscan in % (per Event)", 10000,0 ,100);
183 
184  fPercent_lost_muon_energy=tfs->make<TH1F>("fPercent_lost_muon_energy"," muon energy excluded by dbscan in % (per Event)", 10000,0 ,100);
185  fPercent_lost_electron_energy=tfs->make<TH1F>("fPercent_lost_electron_energy","electron energy excluded by dbscan in % (per Event)", 10000,0 ,100);
186  fPercent_lost_positron_energy=tfs->make<TH1F>("fPercent_lost_positron_energy"," positron energy excluded by dbscan in % (per Event)", 10000,0 ,100);
187  fPercent_lost_111_energy=tfs->make<TH1F>("fPercent_lost_111_energy","pion(111) energy excluded by dbscan in % (per Event)", 10000,0 ,100);
188  fPercent_lost_211_energy=tfs->make<TH1F>("fPercent_lost_211_energy","pion(211) energy excluded by dbscan in % (per Event)", 10000,0 ,100);
189  fPercent_lost_m211_energy=tfs->make<TH1F>("fPercent_lost_m211_energy"," pion(-211) energy excluded by dbscan in % (per Event)", 10000,0 ,100);
190  fPercent_lost_2212_energy=tfs->make<TH1F>("fPercent_lost_2212_energy","proton energy excluded by dbscan in % (per Event)", 10000,0 ,100);
191  fPercent_lost_2112_energy=tfs->make<TH1F>("fPercent_lost_2112_energy","neutron energy excluded by dbscan in % (per Event)", 10000,0 ,100);
192 
193  fEnergy=tfs->make<TH1F>("fEnergy","energy for each voxel", 100000,0 ,0.0005);
194 
195  fbrian_in = tfs->make<TH2F>("fbrian_in", ";# Electrons deposited; # Electrons detected by hitfinder", 1000, 0, 10000000, 1000, 0, 10000000);
196  fbrian_coll = tfs->make<TH2F>("fbrian_coll", ";# Electrons deposited; # Electrons detected by hitfinder", 1000, 0, 10000000, 1000, 0, 10000000);
197 
198  }
199 
200  void DBclusterAna::analyze(const art::Event& evt)
201  {
202 
203  std::cout << "run : " << evt.run() << std::endl;
204  //std::cout << "subrun : " << evt.subRun() << std::endl; // Doesn't compile w. or w.o. id().
205  std::cout << "event : " << evt.id().event() << std::endl;
206  //----------------------------------------------------------------
207 
208  /* This is basically a module for studying MC efficiency/purity. Kick out now if not MC. EC, 8-Oct-2010 */
209  if (evt.isRealData())
210  {
211  std::cout<<"**** DBclusterAna: Bailing. Don't call this module if you're not MC. "<<std::endl;
212  exit (1);
213  }
214 
215  art::ServiceHandle<geo::Geometry> geom;
216  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
217  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
218 
219  art::Handle< std::vector<raw::RawDigit> > rdListHandle;
220  evt.getByLabel(fDigitModuleLabel,rdListHandle);
221  art::Handle< std::vector<recob::Hit> > hitListHandle;
222  evt.getByLabel(fHitsModuleLabel,hitListHandle);
223  art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
224  evt.getByLabel(fGenieGenModuleLabel,mctruthListHandle);
225  art::Handle< std::vector<recob::Cluster> > clusterListHandle;
226  evt.getByLabel(fClusterFinderModuleLabel,clusterListHandle);
227  art::Handle< std::vector<recob::Wire> > wireListHandle;
228  evt.getByLabel(fCalDataModuleLabel,wireListHandle);
229 
230  art::FindManyP<recob::Hit> fmh(clusterListHandle, evt, fClusterFinderModuleLabel);
231 
232  //----------------------------------------------------------------
233 
234  //----------------------------------------------------------------
235 
236  art::PtrVector<raw::RawDigit> rawdigits;
237 
238  for (size_t ii = 0; ii < rdListHandle->size(); ++ii){
239  art::Ptr<raw::RawDigit> rawdigit(rdListHandle,ii);
240  rawdigits.push_back(rawdigit);
241  }
242 
243  //get the simb::MCParticle collection from the art::Event and then use the
244  //Simulation/SimListUtils object to create a sim::ParticleList from the art::Event.
245  pi_serv->SetEveIdCalculator(new sim::EmEveIdCalculator);
246 
247  sim::ParticleList const& _particleList = pi_serv->ParticleList();
248 
249  std::vector<int> mc_trackids;
250 
251  //check_particleList(_particleList[0]);
252  //std::cout<<"checking trackID: ";
253  for ( auto i = _particleList.begin(); i != _particleList.end(); ++i ){
254  int trackID = (*i).first;
255  mc_trackids.push_back(trackID);
256  }
257 
258  // std::cout<<std::endl;
259  // std::cout<<" Track ID list from MC: "<<std::endl;
260  // for(unsigned int i=0;i<mc_trackids.size();++i){
261 
262  // std::cout<<mc_trackids[i]<<" ";
263 
264  // }
265  // std::cout<<"I have in total "<<mc_trackids.size()<<" different tracks"<<std::endl;
266 
267 
268  //----------------------------------------------------------------
269  art::PtrVector<simb::MCTruth> mclist;
270  for (size_t ii = 0; ii < mctruthListHandle->size(); ++ii){
271  art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle,ii);
272  mclist.push_back(mctparticle);
273  }
274 
275  // art::PtrVector<recob::Hit> hits;
276  // for (unsigned int ii = 0; ii < hitListHandle->size(); ++ii)
277  // {
278  // art::Ptr<recob::Hit> hitHolder(hitListHandle,ii);
279  // hits.push_back(hitHolder);
280  // }
281  //--------------------------------------------------
282  std::vector< art::Ptr<recob::Hit> > hits_vec;
283 
284  //--------------------------------------------------
285  std::vector< art::Ptr<recob::Hit> > hits;
286  art::fill_ptr_vector(hits, hitListHandle);
287  //---------------------------------------------------
288 
289  art::PtrVector<recob::Cluster> clusters;
290  for (unsigned int ii = 0; ii < clusterListHandle->size(); ++ii){
291  art::Ptr<recob::Cluster> clusterHolder(clusterListHandle,ii);
292  clusters.push_back(clusterHolder);
293  }
294 
295  std::cout<<"in Efficiency, clusters.size()= "<<clusters.size()<<std::endl;
296 
297  //---------------------------------------------------------------
298  art::PtrVector<recob::Wire> wirelist;
299 
300  for (size_t ii = 0; ii < wireListHandle->size(); ++ii){
301  art::Ptr<recob::Wire> wireHolder(wireListHandle,ii);
302 
303  wirelist.push_back(wireHolder);
304 
305  }
306 
307 
308 
309  //...........................................................................
310  // How many different particles do we have in a cluster???
311  // -- I will answer that by 3 methods
312  // 1) count different TrackIDs in each cluster
313  // 2) count different TrackIDs in each cluster with the use of "Mother" to get rid of the ones that just randomly got their TrackIDs changed by Geant4
314  // 3) count different pdg codes in each cluster <--probably not very good b/c if you have 2 "different" electrons in a cluster it's going to count them as one.
315  //..........................................................................
316  // How many clusters it takes to contain a particle???
317  // take each TrackID and count in how many clusters it appears
318  //.........................................................................
319 
320  double no_of_particles_in_cluster=0;
321  double sum_vec_trackid=0;
322  double no_of_clusters=0;
323  double total_no_hits_in_clusters=0;
324  //unsigned int plane=0;
325  art::Ptr<raw::RawDigit > _rawdigit;
326  art::Ptr<raw::RawDigit > _rawdigit2;
327  std::vector<int> vec_pdg;
328  std::vector<int> vec_trackid,vec_trackid_mother, vec_trackid_mother_en;
329  std::vector<int> all_trackids;
330  std::vector<int> ids;
331  std::vector<int>::iterator it,it2,it3,it4,it5,it6,it7,it8;
332  int no_cl_for_muon=0;
333  int no_cl_for_electron=0;
334  int no_cl_for_positron=0;
335  int no_cl_for_pion_111=0;
336  int no_cl_for_pion_211=0;
337  int no_cl_for_pion_m211=0;
338  int no_cl_for_proton=0;
339  double noCluster=0;
340  // int muon=0,electron=0,positron=0,pion=0;
341  double _hit_13=0,_hit_11=0,_hit_m_11=0,_hit_111=0,_hit_211=0,_hit_m211=0,_hit_2212=0,_hit_2112=0;
342  double _en_13=0,_en_11=0,_en_m11=0,_en_111=0,_en_211=0,_en_m211=0,_en_2212=0,_en_2112=0;
343  std::vector<double> diff_vec;
344 
345  double hit_energy=0;
346  double total_Q_cluster_hits=0;
347 
348 
349  /*
350  for(unsigned int i = 0; i < hits.size(); ++i) {
351  std::cout<<"channel: "<<hits[i]->Wire()->RawDigit()->Channel()<<" time= "<<(hits[i]->StartTime()+hits[i]->EndTime())/2.<<" X time= "<<hits[i]-> CrossingTime()<<std::endl;
352  }
353  */
354 
355  if(clusters.size()!=0 && hits.size()!=0){
356  for(unsigned int plane=0;plane<geom->Nplanes();++plane){
357  geo::View_t view = geom->Plane(plane).View();
358  //art::PtrVector<recob::Cluster>::const_iterator clusterIter = clusters.begin();
359  for(size_t j = 0; j < clusters.size(); ++j){
360 
361  // std::cout<<"I AM ON PLANE #"<<plane<<std::endl;
362  if( clusters[j]->View() == view){
363 
364  //std::cout<<"working on cluster # "<<j<<std::endl;
365  std::vector< art::Ptr<recob::Hit> > _hits = fmh.at(j);
366  art::Ptr<recob::Hit> _hits_ptr; //
367 
368  //delete
369  //std::cout<<"_hits.size()= "<<_hits.size()<<std::endl;
370  for(size_t p = 0; p<_hits.size(); ++p){
371  _hits_ptr=_hits[p];
372  hits_vec.push_back(_hits_ptr);
373  //std::cout<<"hit # "<<p<<" charge= "<<_hits[p]->Integral()<<std::endl;
374  total_Q_cluster_hits += _hits[p]->Integral();
375  }
376 
377 
378  std::vector< art::Ptr<recob::Hit> >::iterator itr = hits_vec.begin();
379 
380 
381  //std::cout<<"hits_vec.size()= "<<hits_vec.size()<<std::endl;
382  while(itr != hits_vec.end()) {
383  //std::cout<<"working on hit # "<<itr-hits_vec.begin()<<" charge= "<<_hits[itr-hits_vec.begin()]->Integral()<<std::endl;
384  diff_vec.clear();
385  //std::cout<<"same?, q= "<<hits_vec[itr-hits_vec.begin()]->Integral()<<std::endl;
386 
387  hit_energy=_hits[itr-hits_vec.begin()]->Integral();
388 
389  std::vector<sim::TrackIDE> trackides = bt_serv->HitToTrackIDEs(*itr);
390 
391  std::vector<sim::TrackIDE> eveides = bt_serv->HitToEveTrackIDEs(*itr);
392 
393  std::vector<sim::TrackIDE>::iterator idesitr = trackides.begin();
394 
395  while( idesitr != trackides.end() ){
396 
397  //int eveID = _particleList.EveId( (*idesitr).trackID );
398 
399  // std::cout<<"track id: " << (*idesitr).trackID<<" contributed " << (*idesitr).energyFrac<< " to the current hit and has eveID: " << eveID<<std::endl;
400 
401  // double energy=voxelData.Energy(i);
402  // std::cout<<"check4-3"<<std::endl;
403  // std::cout<<"energy= "<<energy<<std::endl;
404  // fEnergy->Fill(energy);
405 
406  vec_trackid.push_back((*idesitr).trackID);
407 
408  //for( unsigned int i=0; i<_particleList.size(); ++i )
409  // {
410  // // const sim::ParticleList* particleList = _particleList[i];
411  // //particleList = _particleList[i];
412  //
413  //
414  // // double energyTrackID=voxelData.Energy[trackID];
415  // // std::cout<<"ENERGY OF PRIMARY TRACKID= "<<energyTrackID<<std::endl;
416  const simb::MCParticle* particle = _particleList.at( (*idesitr).trackID);
417  //
418  int pdg = particle->PdgCode();
419  // std::cout<<"pdg= "<<pdg<<std::endl;
420  //
421  // double energy2=voxelData.Energy(i);
422  // // std::cout<<"energy2= "<<energy2<<std::endl;
423  //
424  // // std::cout<<"part eng= "<<particle->E()<<std::endl;
425  if(pdg==13 || pdg==-13){
426  _hit_13++;
427  _en_13+=hit_energy*((*idesitr).energyFrac);
428  // std::cout<<"_en_13= "<<_en_13<<std::endl;
429  }
430  // if(pdg==11){_hit_11++;
431  // _en_11+=energy2;
432  // // std::cout<<"in clus: _en_11="<<_en_11<<std::endl;
433  // }
434  // if(pdg==-11){_hit_m_11++;
435  // _en_m11+=energy2;}
436  // if(pdg==111){_hit_111++;
437  // _en_111+=energy2;}
438  // if(pdg==22){_hit_22++;
439  // _en_22+=energy2;}
440  // if(pdg==211){_hit_211++;
441  // _en_211+=energy2;}
442  // if(pdg==-211){_hit_m211++;
443  // _en_m211+=energy2;}
444  // if(pdg==2212){_hit_2212++;
445  // _en_2212+=energy2;}
446  // if(pdg==2112){_hit_2112++;
447  // _en_2112+=energy2;}
448  //
449  // //std::cout<<"True PDG= "<<pdg<<std::endl;
450  // vec_pdg.push_back(pdg);
451  // // std::cout<<"_en_11= "<<_en_11<<std::endl;
452  // //while particle is not a primary particle and going up in a chain of trackIDs is not going to change its pdg code, go up the chain.
453  // while ( (! _particleList.IsPrimary( trackID )) && (((_particleList.at(particle->Mother()))->PdgCode())==pdg))
454  // {
455  // trackID = particle->Mother();
456  // //std::cout<<"((NOt a PRIMARY ORIGINALLY!!! ) trackID= "<<trackID<<std::endl;
457  // particle = _particleList.at( trackID );
458  // pdg= particle->PdgCode();
459  // // std::cout<<"(NOt a PRIMARY ORIGINALLY!!! ) The PDG from HIT is: "<<pdg<<std::endl;
460  //
461  // }
462  //
463  // // std::cout<<"The PDG from HIT is: "<<pdg<<std::endl;
464  // //std::cout<<"after mother trackid= "<<trackID<<std::endl;
465  // vec_trackid_mother.push_back(trackID);
466  // if(energy>(7e-5)){ vec_trackid_mother_en.push_back(trackID);}
467  // }
468 
469 
470 
471  idesitr++;
472 
473  }
474 
476 
477  // int numberPrimaries = particleList->NumberOfPrimaries();
478  // std::cout<<"no of PRIMARIES: "<<numberPrimaries<<std::endl;
479  // for ( int i = 0; i != numberPrimaries; ++i )
480  // {
481  // const simb::MCParticle* primaryParticle = particleList->Primary(i);
482  // int trackID = primaryParticle->TrackId();
483  // std::cout<<"from PRIMARY trackID= "<<trackID<<std::endl;
484  // double energy=voxelData.Energy[trackID];
485  // std::cout<<"ENERGY OF PRIMARY TRACKID= "<<energy<<std::endl;
486  // }
487 
489 
490  // std::cout<<"check--3"<<std::endl;
491 
492  itr++;
493 
494  }//loop thru hits
495  //std::cout<<"check--4"<<std::endl;
496 
497  // std::cout<<"vec_pdg("<<vec_pdg.size()<<")= " ;
498  for(unsigned int i=0;i<vec_pdg.size();++i){
499 
500  // std::cout<<vec_pdg[i]<<" ";
501 
502  }
503  //std::cout<<std::endl;
504  // std::cout<<"vec_trackid("<<vec_trackid.size()<<")= ";
505  // for(unsigned int ii=0;ii<vec_trackid.size();++ii){
506 
507  // std::cout<<vec_trackid[ii]<<" ";
508 
509  // }
510 
511  // std::cout<<"vec_trackid_mother("<<vec_trackid_mother.size()<<")= ";
512  // for(unsigned int ii=0;ii<vec_trackid_mother.size();++ii){
513 
514  // std::cout<<vec_trackid_mother[ii]<<" ";
515 
516  // }
517 
518 
519  it=find(vec_pdg.begin(),vec_pdg.end(),13);
520  if(it!=vec_pdg.end()){
521  // std::cout<<"matched found at position="<<int(it-vec_pdg.begin())<<std::endl;
522  no_cl_for_muon++;
523  }
524  else{
525  // std::cout<<"no match!"<<std::endl;
526  }
527 
528  it2=find(vec_pdg.begin(),vec_pdg.end(),11);
529  if(it2!=vec_pdg.end()){
530  no_cl_for_electron++;
531  }
532 
533  it3=find(vec_pdg.begin(),vec_pdg.end(),-11);
534  if(it3!=vec_pdg.end()){
535  no_cl_for_positron++;
536  }
537 
538  it4=find(vec_pdg.begin(),vec_pdg.end(),111);
539  if(it4!=vec_pdg.end()){
540  no_cl_for_pion_111++;
541  }
542  it6=find(vec_pdg.begin(),vec_pdg.end(),211);
543  if(it6!=vec_pdg.end()){
544  no_cl_for_pion_211++;
545  }
546  it7=find(vec_pdg.begin(),vec_pdg.end(),-211);
547  if(it7!=vec_pdg.end()){
548  no_cl_for_pion_m211++;
549  }
550  it8=find(vec_pdg.begin(),vec_pdg.end(),2212);
551  if(it8!=vec_pdg.end()){
552  no_cl_for_proton++;
553  }
554 
555  // std::cout<<std::endl;
556  //std::cout<<"numberParticles= "<<numberParticles<<std::endl;
557  // std::cout<<"size of vec_pdg= "<<vec_pdg.size()<<std::endl;
558  sort( vec_pdg.begin(), vec_pdg.end() );
559  vec_pdg.erase( unique( vec_pdg.begin(), vec_pdg.end() ), vec_pdg.end() );
560  // std::cout<<" NO OF PARTICLES IN THIS CLUSTER IS: "<<vec_pdg.size()<<std::endl;
561  // std::cout<<"They are: ";
562  // for(unsigned int ii=0;ii<vec_pdg.size();++ii){
563  // std::cout<<vec_pdg[ii]<<" ";
564  // }
565  // std::cout<<std::endl;
566 
567  //same for vec_trackid:
568 
569  sort( vec_trackid.begin(), vec_trackid.end() );
570  vec_trackid.erase( unique( vec_trackid.begin(), vec_trackid.end() ), vec_trackid.end() );
571  // std::cout<<" NO OF DIFFERENT TRACKIDS IN THIS CLUSTER IS: "<<vec_trackid.size()<<std::endl;
572  // std::cout<<"They are: ";
573  for(unsigned int ii=0;ii<vec_trackid.size();++ii){
574  // std::cout<<vec_trackid[ii]<<" ";
575  all_trackids.push_back(vec_trackid[ii]);
576  //mytracklist.push_back(vec_trackid[i]);
577  }
578  // std::cout<<std::endl;
579  //...............................................................
580  //Also Make vec_trackid_mother unique:
581 
582  sort( vec_trackid_mother.begin(), vec_trackid_mother.end() );
583  vec_trackid_mother.erase( unique( vec_trackid_mother.begin(), vec_trackid_mother.end() ), vec_trackid_mother.end() );
584  // std::cout<<" NO OF DIFFERENT TRACKIDS_MOTHER IN THIS CLUSTER IS: "<<vec_trackid_mother.size()<<std::endl;
585  // std::cout<<"They are: ";
586  /*
587  for(unsigned int ii=0;ii<vec_trackid_mother.size();++ii){
588  std::cout<<vec_trackid_mother[ii]<<" ";
589 
590  }
591  */
592  // std::cout<<std::endl;
593 
594  //........................................................................
595 
596  //Also Make vec_trackid_mother_en unique:
597 
598  sort( vec_trackid_mother_en.begin(), vec_trackid_mother_en.end() );
599  vec_trackid_mother_en.erase( unique( vec_trackid_mother_en.begin(), vec_trackid_mother_en.end() ), vec_trackid_mother_en.end() );
600  // std::cout<<" NO OF DIFFERENT TRACKIDS_MOTHER_en IN THIS CLUSTER IS: "<<vec_trackid_mother_en.size()<<std::endl;
601  // std::cout<<"They are: ";
602  /*
603  for(unsigned int ii=0;ii<vec_trackid_mother_en.size();++ii){
604  std::cout<<vec_trackid_mother_en[ii]<<" ";
605 
606  }
607  */
608  // std::cout<<std::endl;
609 
610  //........................................................................
611 
612 
613  // Q: How many clusters it takes to contain a certain particle?
614  for(unsigned int ii=0;ii<mc_trackids.size();++ii){
615  it5=find(vec_trackid.begin(),vec_trackid.end(),mc_trackids[ii]);
616  if(it5!=vec_trackid.end()){
617  // std::cout<<"found match for: "<<mc_trackids[ii]<<" at position: "<<it5-vec_trackid.begin()<<std::endl;
618  ids.push_back(mc_trackids[ii]);//then make it unique
619  noCluster++;
620  }
621  }
622 
623  // std::cout<<"noCluster = "<< noCluster<<std::endl;
624 
625 
626 
627 
628  fNoParticles_pdg->Fill(vec_pdg.size());
629 
630  // std::cout<<"LOOK ->> vec_trackid.size()= "<<vec_trackid.size()<<std::endl;
631  fNoParticles_trackid->Fill(vec_trackid.size());
632 
633  fNoParticles_trackid_mother->Fill(vec_trackid_mother.size());
634 
635  no_of_clusters++;
636  // std::cout<<"MUON IS CONTAINED IN "<<no_cl_for_muon<<" CLUSTERS"<<std::endl;
637 
638  no_of_particles_in_cluster+=vec_pdg.size();
639  sum_vec_trackid+=vec_trackid_mother.size();
640  total_no_hits_in_clusters+=_hits.size();
641 
642  vec_pdg.clear();
643 
644  vec_trackid.clear();
645 
646  vec_trackid_mother.clear();
647 
648  vec_trackid_mother_en.clear();
649 
650 
651  }//end if cluster is in correct view
652  //clusterIter++;
653 
654  hits_vec.clear();
655  }//for each cluster
656 
657  // std::cout<<"sum_vec_trackid= "<<sum_vec_trackid<<std::endl;
658 
659  sort( all_trackids.begin(), all_trackids.end() );
660  all_trackids.erase( unique( all_trackids.begin(), all_trackids.end() ), all_trackids.end() );
661  // std::cout<<" NO OF DIFFERENT TRACKIDS IN THIS EVENT IS: "<<all_trackids.size()<<std::endl;
662  // std::cout<<"They are: ";
663  for(unsigned int ii=0;ii<all_trackids.size();++ii){
664  // std::cout<<all_trackids[ii]<<" ";
665 
666  }
667  // std::cout<<std::endl;
668 
669  //now I have a vector(all_trackids) that only contains unique trackids
670  //go to it and search for each trackid in every cluster
671 
672 
673 
674  sort(ids.begin(),ids.end() );
675  ids.erase( unique( ids.begin(), ids.end() ), ids.end() );
676  // std::cout<<"**NEW** NO OF USED TRACKIDS IN THIS EVENT IS: "<<ids.size()<<std::endl;
677  double no_of_clusters_per_track=noCluster/ids.size();
678  // std::cout<<"alright so i have no_of_clusters_per_track= "<<no_of_clusters_per_track<<std::endl;
679 
680 
681 
682  //Filling Histograms:
683 
684  if(no_cl_for_muon!=0){
685  fCl_for_Muon->Fill(no_cl_for_muon);}
686  // if(no_cl_for_electron!=0){
687  // fCl_for_Electron->Fill(no_cl_for_electron);}
688  // if(no_cl_for_positron!=0){
689  // fCl_for_Positron->Fill(no_cl_for_positron);}
690  // if(no_cl_for_pion_111!=0){
691  // fCl_for_Pion_111->Fill(no_cl_for_pion_111);}
692  // if(no_cl_for_pion_211!=0){
693  // fCl_for_Pion_211->Fill(no_cl_for_pion_211);}
694  // if(no_cl_for_pion_m211!=0){
695  // fCl_for_Pion_m211->Fill(no_cl_for_pion_m211);}
696  // if(no_cl_for_proton!=0){
697  // fCl_for_Proton->Fill(no_cl_for_proton);}
698 
699  fno_of_clusters_per_track->Fill(no_of_clusters_per_track);
700 
701  no_cl_for_muon=0;
702  no_cl_for_electron=0;
703  no_cl_for_positron=0;
704  no_cl_for_pion_111=0;
705  no_cl_for_pion_211=0;
706  no_cl_for_pion_m211=0;
707  no_cl_for_proton=0;
708  noCluster=0;
709  // std::cout<<"*************************************************"<<std::endl;
710  // std::cout<<"*************************************************"<<std::endl;
711  // std::cout<<" NEW PLANE "<<std::endl;
712  // std::cout<<"*************************************************"<<std::endl;
713  // std::cout<<"*************************************************"<<std::endl;
714  ids.clear();
715 
716  }//for each plane
717  //std::cout<<"no_of_particles_in_cluster= "<<no_of_particles_in_cluster<<std::endl;
718  //std::cout<<" no_of_clusters (with hits that are non-zero)= "<< no_of_clusters<<std::endl;
719  double result=no_of_particles_in_cluster/no_of_clusters;
720 
721  // std::cout<<"FINALLY NO OF PARTICLES PER CLUSTER IS: "<<result<<std::endl;
722  // std::cout<<"Sum of all hits in clusters= "<<total_no_hits_in_clusters<<std::endl;
723  double no_noise_hits=hits.size()-total_no_hits_in_clusters;
724  // std::cout<<"No of hits marked as noise= "<< no_noise_hits<<std::endl;
725  double percent_noise=double(no_noise_hits/hits.size())*100;
726  // std::cout<<"%%%%%%%%% NOISE IS: "<< percent_noise<<std::endl;
727 
728 
729 
730  double no_trackid_per_cl_per_event = sum_vec_trackid/no_of_clusters;
731  // std::cout<<"sum_vec_trackid= "<<sum_vec_trackid<<" no_of_clusters: "<<no_of_clusters<<" no_trackid_per_cl_per_event= "<<no_trackid_per_cl_per_event<<std::endl;
732  fNoParticles_trackid_per_event->Fill(no_trackid_per_cl_per_event);
733  fNoParticles_pdg_per_event->Fill(result);
734  fNoClustersInEvent->Fill(clusters.size());
735 
736  fPercentNoise->Fill(percent_noise);
737 
738  //????????here clean sum_vec_trackid
739 
741  //-----------------------------------------------------------------
742 
743 
744 
745 
746 
747  //-----------------------------------------------------------------------
748  for( unsigned int i = 0; i < mclist.size(); ++i ){
749  //art::Ptr<const simb::MCTruth> mc(mctruthListHandle,i);
750  art::Ptr<simb::MCTruth> mc(mclist[i]);
751  // simb::MCTruth mcp(*(mclist[i]));// We don't ever use this. EC, 6-Oct-2010.
752  // std::cout<<"Number of particles is : "<<mc->NParticles()<<std::endl;
753  // std::cout<<" mc->NParticles()="<< mc->NParticles()<<std::endl;
754  // TParticle part(mc->GetParticle(1));
755  // std::cout<<"part.Gte<PDG()->PdgCode()= "<<part.GetPDG()->PdgCode()<<std::endl;
756 
757  for(int ii = 0; ii < mc->NParticles(); ++ii){
758  simb::MCParticle part(mc->GetParticle(ii));
759  std::cout<<"FROM MC TRUTH,the particle's pdg code is: "<<part.PdgCode()<<std::endl;
760  std::cout<<"with energy= "<<part.E();
761  if(abs(part.PdgCode()) == 13){std::cout<<" I have a muon!!!"<<std::endl;
762  // std::cout<<"with energy= "<<part.Energy();
763  }
764  if(abs(part.PdgCode()) == 111){std::cout<<" I have a pi zero!!!"<<std::endl;}
765 
766  }
767 
768  }
769 
770  //-----------------------------------------------------------------------
771 
772 
773 
774 
775 
776  }
777 
778  // std::cout<<std::endl;
779  // std::cout<<" ********************************************************"<<std::endl;
780  // std::cout<<" ************* HITS ONLY *****************************"<<std::endl;
781  // std::cout<<" 88888888888888888888888888888888888888888888888888888888"<<std::endl;
782  //......................................................................
783  //......................................................................
784  //......................................................................
785  // *** ***************** *** *****************
786 
787  // Now I can load just hits (not clusters) and see how many of what kind are lost due to dbscan marking them as noise:
788 
789  //Load hits and copy most of the code from above:
790 
791  // std::cout<<"Take care of "<<hits.size()<<" hits"<<std::endl;
792  // std::cout<<"_en_11= "<<_en_11<<" _en_13= "<<_en_13<<std::endl;
793  double hit_13=0,hit_11=0,hit_m_11=0,hit_111=0,hit_211=0,hit_m211=0,hit_2212=0,hit_2112=0;
794  double en_13=0,en_11=0,en_m11=0,en_111=0,en_211=0,en_m211=0,en_2212=0,en_2112=0;
795  //int no_hits=0;
796  //unsigned int plane_k=0;
797  //double total_eng_hits_p0=0;
798  //double total_eng_hits_p1=0;
799 
800  // geo::View_t view_ind = geom->Plane(0).View();
801  // geo::View_t view_coll = geom->Plane(1).View();
802 
803  std::vector< art::Ptr<recob::Hit> >::iterator itr = hits.begin();
804  while(itr != hits.end()) {
805 
806  std::vector<sim::TrackIDE> trackides = bt_serv->HitToTrackIDEs(*itr);
807  std::vector<sim::TrackIDE> eveides = bt_serv->HitToEveTrackIDEs(*itr);
808 
809  std::vector<sim::TrackIDE>::iterator idesitr = trackides.begin();
810 
811  hit_energy=hits[itr-hits.begin()]->Integral();
812 
813  while( idesitr != trackides.end() ){
814 
815  //std::cout<<"0:TOTAL ENERGY FROM HITS for P=0 = "<<total_eng_hits_p0<<std::endl;
816  //std::cout<<"0:TOTAL ENERGY FROM HITS for P=1 = "<<total_eng_hits_p1<<std::endl;
817 
818 
819  const simb::MCParticle* particle = _particleList.at( (*idesitr).trackID);
820 
821  int pdg = particle->PdgCode();
822  //std::cout<<"pdg= "<<pdg<<std::endl;
823 
824 
825 
826 
827  diff_vec.clear();
828 
829 
830 
831 
832 
833  // double energy3=voxelData.Energy(i);
834  // //std::cout<<"plane= "<<plane_k<<std::endl;
835  // if(plane_k==0){total_eng_hits_p0+=energy3;}
836  // if(plane_k==1){total_eng_hits_p1+=energy3;}
837 
838 
839  if(pdg==13 || pdg==-13){hit_13++;
840  en_13+=hit_energy*((*idesitr).energyFrac);}
841  // if(pdg==11){hit_11++;
842  // en_11+=energy3;
843  // // std::cout<<"in hits: en_11="<<en_11<<std::endl;
844  // }
845  // if(pdg==-11){hit_m_11++;
846  // en_m11+=energy3;}
847  // if(pdg==111){hit_111++;
848  // en_111+=energy3;}
849  // if(pdg==22){hit_22++;
850  // en_22+=energy3;}
851  // if(pdg==211){hit_211++;
852  // en_211+=energy3;}
853  // if(pdg==-211){hit_m211++;
854  // en_m211+=energy3;}
855  // if(pdg==2212){hit_2212++;
856  // en_2212+=energy3;}
857  // if(pdg==2112){hit_2112++;
858  // en_2112+=energy3;}
859  // if(pdg !=22 && pdg!=111 && pdg!=-11 && pdg !=11 && pdg!=13 && pdg!=211 && pdg!=-211 && pdg!=2212 && pdg!=2112){
860  //
861  // std::cout<<"SOMETHING ELSE!!! PDG= "<<pdg<<std::endl;
862  // }
863 
864  idesitr++;
865 
866  } //trackIDs
867 
868  itr++;
869  } //hits
870  // std::cout<<"True PDG= "<<pdg<<std::endl;
871 
872 
873 
874  // std::cout<<"hit_13= "<<hit_13<<" "<<"hit_11= "<<hit_11<<" "<<"hit_m_11= "<<hit_m_11<<" "<<"hit_111= "<<hit_111<<" "<<"hit_22= "<<hit_22<<" ";
875 
876  // int sum=hit_13+hit_11+hit_m_11+hit_111+hit_22;
877  // std::cout<<"sum= "<<sum<<" no_hits= "<<no_hits<<" DIFF= "<<sum-no_hits<<std::endl;
878 
879 
880 
881 
882  // std::cout<<"PLANE_K= "<<plane_k<<std::endl;
883 
884 
885 
886  //std::cout<<"TOTAL ENERGY FROM HITS for P=0 = "<<total_eng_hits_p0<<std::endl;
887  //std::cout<<"TOTAL ENERGY FROM HITS for P=1 = "<<total_eng_hits_p1<<std::endl;
888 
889 
890 
891 
892  // std::cout<<"After hits,PLANE_K= "<<plane_k<<std::endl;
893  // }//plane
894  // }//non-zero hits
895  // std::cout<<"no_hits= "<<no_hits<<std::endl;
896  // std::cout<<"TOTAL # of muon hits= "<<hit_13<<std::endl;
897  // std::cout<<"TOTAL # of electron hits= "<<hit_11<<std::endl;
898  // std::cout<<"TOTAL # of positron hits= "<<hit_m_11<<std::endl;
899  // std::cout<<"TOTAL # of 111 hits= "<<hit_111<<std::endl;
900  // std::cout<<"TOTAL # of 211 hits= "<<hit_211<<std::endl;
901  // std::cout<<"TOTAL # of -211 hits= "<<hit_m211<<std::endl;
902  // std::cout<<"TOTAL # of 2212 hits= "<<hit_2212<<std::endl;
903  // std::cout<<"IN CLUSTERS # of muon hits= "<<_hit_13<<std::endl;
904  // std::cout<<"IN CLUSTERS # of electron hits= "<<_hit_11<<std::endl;
905  // std::cout<<"IN CLUSTERS # of positron hits= "<<_hit_m_11<<std::endl;
906  // std::cout<<"IN CLUSTERS # of 111 hits= "<<_hit_111<<std::endl;
907  // std::cout<<"IN CLUSTERS # of 211 hits= "<<_hit_211<<std::endl;
908  // std::cout<<"IN CLUSTERS # of -211 hits= "<<_hit_m211<<std::endl;
909  // std::cout<<"IN CLUSTERS # of 2212 hits= "<<_hit_2212<<std::endl;
910 
911  // std::cout<<"WE MISSED % of muon hits= "<<100-((_hit_13/hit_13)*100)<<"%"<<std::endl;
912 
913  // std::cout<<"WE MISSED % of electron hits= "<<100-((_hit_11/hit_11)*100)<<"%"<<std::endl;
914  // std::cout<<"WE MISSED % of positron hits= "<<100-((_hit_m_11/hit_m_11)*100)<<"%"<<std::endl;
915  // std::cout<<"WE MISSED % of 111 hits= "<<100-((_hit_111/hit_111)*100)<<"%"<<std::endl;
916  // std::cout<<"WE MISSED % of 211 hits= "<<100-((_hit_211/hit_211)*100)<<"%"<<std::endl;
917  // std::cout<<"WE MISSED % of -211 hits= "<<100-((_hit_m211/hit_m211)*100)<<"%"<<std::endl;
918  // std::cout<<"WE MISSED % of 2212 hits= "<<100-((_hit_2212/hit_2212)*100)<<"%"<<std::endl;
919 
920  if(hit_13!=0){ fPercent_lost_muon_hits->Fill(100-((_hit_13/hit_13)*100));}
921  if(hit_11!=0){ fPercent_lost_electron_hits->Fill(100-((_hit_11/hit_11)*100));}
922  if(hit_m_11!=0){ fPercent_lost_positron_hits->Fill(100-((_hit_m_11/hit_m_11)*100));}
923  if(hit_111!=0){ fPercent_lost_111_hits->Fill(100-((_hit_111/hit_111)*100));}
924 
925  if(hit_211!=0){ fPercent_lost_211_hits->Fill(100-((_hit_211/hit_211)*100));}
926 
927  if(hit_m211!=0){ fPercent_lost_m211_hits->Fill(100-((_hit_m211/hit_m211)*100));}
928  if(hit_2212!=0){ fPercent_lost_2212_hits->Fill(100-((_hit_2212/hit_2212)*100));}
929  if(hit_2112!=0){ fPercent_lost_2112_hits->Fill(100-((_hit_2112/hit_2112)*100));}
930 
931  // std::cout<<"*** _en_11= "<<_en_11<<" en_11= "<<en_11<<std::endl;
932  // std::cout<<"WE MISSED % of muon energy= "<<100-((_en_13/en_13)*100)<<"%"<<std::endl;
933 
934  // std::cout<<"WE MISSED % of electron energy= "<<100-((_en_11/en_11)*100)<<"%"<<std::endl;
935  // std::cout<<"WE MISSED % of positron energy= "<<100-((_en_m11/en_m11)*100)<<"%"<<std::endl;
936  // std::cout<<"WE MISSED % of 111 energy= "<<100-((_en_111/en_111)*100)<<"%"<<std::endl;
937  // std::cout<<"WE MISSED % of 211 energy= "<<100-((_en_211/en_211)*100)<<"%"<<std::endl;
938  // std::cout<<"WE MISSED % of -211 energy= "<<100-((_en_m211/en_m211)*100)<<"%"<<std::endl;
939  // std::cout<<"WE MISSED % of 2212 energy= "<<100-((_en_2212/en_2212)*100)<<"%"<<std::endl;
940  // std::cout<<"WE MISSED % of 2112 energy= "<<100-((_en_2112/en_2112)*100)<<"%"<<std::endl;
941  // if(en_13==0){std::cout<<"NO MU IN THIS EVENT (en) $$$$$$$$$$$$$$$$$"<<std::endl;}
942  //if(hit_13==0){std::cout<<"NO MU IN THIS EVENT (hit)$$$$$$$$$$$$$$$$$"<<std::endl;}
943  std::cout<<"****** mu E from clusters = "<<_en_13<<std::endl;
944  std::cout<<"****** mu E from hits = "<<en_13<<std::endl;
945 
946  if(en_13!=0){ fPercent_lost_muon_energy->Fill(100-((_en_13/en_13)*100));}
947  if(en_11!=0){ fPercent_lost_electron_energy->Fill(100-((_en_11/en_11)*100));}
948  if(en_m11!=0){ fPercent_lost_positron_energy->Fill(100-((_en_m11/en_m11)*100));
949  // std::cout<<"POSITRON E= "<<100-((_en_m11/en_m11)*100)<<std::endl;
950  }
951  if(en_111!=0){ fPercent_lost_111_energy->Fill(100-((_en_111/en_111)*100));}
952  if(en_211!=0){ fPercent_lost_211_energy->Fill(100-((_en_211/en_211)*100));}
953  if(en_m211!=0){ fPercent_lost_m211_energy->Fill(100-((_en_m211/en_m211)*100));}
954  if(en_2212!=0){ fPercent_lost_2212_energy->Fill(100-((_en_2212/en_2212)*100));}
955  if(en_2112!=0){ fPercent_lost_2112_energy->Fill(100-((_en_2112/en_2112)*100));}
956 
957 
958 
959 
960 
962  //-----------------------------------------------------------------------
963 
964  // FOR BRIAN:
965 
966 
967  //-------------------------------------------------------------------
968 
969  //std::cout<<"hello for Brian-->>>"<<std::endl;
970  //------------------------------------------------------------
971  // std::vector<const recob::Wire*> wirelist;
972  // try{
973  // evt.Reco().Get(fInputFolder.c_str(),wirelist);
974  // }
975  // catch(art::Exception e){
976  // std::cerr << "Error retrieving wire list, while looking for wires "
977  // << "in hit::FFTHitFinder::Ana(), "<< "directory : "
978  // << fInputFolder.c_str() << std::endl;
979  // return jobc::kFailed;
980  // }
981  // _electrons=0;
982  // electrons=0;
983  // int Total_Elec_p0=0;
984  // int Total_Elec_p1=0;
985  // int wire;
986  // int pl=0,t;
987  // unsigned int channel;
988  // //loop through all wires:
989 
990  // for(std::vector<const recob::Wire*>::iterator wireIter = wirelist.begin();
991  // wireIter != wirelist.end(); wireIter++) {
992 
993 
994 
995  // _rawdigit2 = (*wireIter)->RawDigit();
996  // sim::SimDigit* simdigit = dynamic_cast< sim::SimDigit*>(_rawdigit2);
997  // int numberOfElectrons = simdigit->NumberOfElectrons();
998 
999  // if(numberOfElectrons==0){std::cout<<" ZERO ELEC!!!"<<std::endl;}
1000  // //std::cout<<"# of elec: "<<numberOfElectrons<<" ";
1001  // if(simdigit==0){std::cout<<"simdigit=0 !!!!!!!!!!"<<std::endl;}
1002 
1003  // if(pl==0) {
1004  // Total_Elec_p0 += numberOfElectrons;}
1005  // if(pl==1) {
1006  // Total_Elec_p1 += numberOfElectrons;}
1007 
1008 
1009 
1010  // }//loop wires
1011 
1012 
1013  // std::cout<<"NO OF ELECTRONS, p0 = "<< Total_Elec_p0<<std::endl;
1014  // std::cout<<"NO OF ELECTRONS, p1 = "<< Total_Elec_p1<<std::endl;
1015 
1016  //now determine number of electrons for hits
1017 
1018 
1019  //----------------------------------------------------------------------------------------
1020  // Now, do the same thing for the found hits and compare the number of ionization electrons.
1021 
1022  //----------------------------------------------------------------------------------------
1023 
1024 
1025  // _electrons=0;
1026  // electrons=0;
1027  //
1028  // double sum=0;
1029  // double sum0=0;
1030  //
1031  //
1032  //
1033  // double no_ele_p0=0;
1034  // double no_ele_p1=0;
1035  // unsigned int plane=0;
1036  // double Tno_ele_p0=0;
1037  // double Tno_ele_p1=0;
1038 
1039 
1040  //for(unsigned int j = 0; j < hits.size(); ++j) // {
1041  //
1042  //
1043  // unsigned int channel = hits[j]->Wire()->RawDigit()->Channel();
1044  //
1045  // // std::cout<<"channel= "<<w_<<std::endl;
1046  // double XTime=hits[j]->PeakTime();
1047  //
1048  // // loop over the SimChannels to find this one
1049  // art::Ptr<sim::SimChannel> sc2;
1050  // for(unsigned int scs = 0; scs < simchans.size(); ++scs)
1051  // if(simchans[scs]->Channel() == channel) sc2 = simchans[scs];
1052  //
1053  // unsigned int numberOfElectrons = sc2->NumberOfElectrons();
1054  //
1055  // // if(numberOfElectrons==0){std::cout<<" ZERO ELEC!!!"<<std::endl;}
1056  // // std::cout<<"# of elec: "<<"for plane: "<<plane<<" is: "<<numberOfElectrons<<std::endl;
1057  // // std::cout<<"simdigit is: "<<simdigit<<std::endl;
1058  // for (size_t i = 0; i != numberOfElectrons; ++i )
1059  // {
1060  //
1061  // _electrons = sc2->GetElectrons(i);
1062  //
1063  // double ArrivalTime=(_electrons->ArrivalT())/200;
1064  //
1065  // double diff=XTime-ArrivalTime;
1066  // // std::cout<<"e's ArrivalT = "<<ArrivalTime<<" diff= "<<diff<<std::endl;
1067  //
1068  // if(plane==0) {
1069  // if((diff<22)&&(diff>13))
1070  // {
1071  //
1072  // electrons = sc2->GetElectrons(i);
1073  //
1074  // no_ele_p0= electrons-> NumElectrons();
1075  // // std::cout<<"p0,channel= "<<w_<<" e= "<<no_ele_p0<<" Hits "<<std::endl;
1076  // Tno_ele_p0+= no_ele_p0;
1077  // sum0=Tno_ele_p0;
1078  //
1079  // //std::cout<<"sum= "<<sum0<<std::endl;
1080  //
1081  // }
1082  // //std::cout<<" no_ele_p0= "<< no_ele_p0<<std::endl;
1083  // }
1084  //
1085  //
1086  //
1087  // if(plane==1) {
1088  // if((diff<36)&&(diff>27))
1089  // {
1090  // electrons = sc2->GetElectrons(i);
1091  // no_ele_p1= electrons-> NumElectrons();
1092  // //double _ArrivalTime=(electrons->ArrivalT())/200;
1093  // //double _diff=XTime-_ArrivalTime;
1094  // // std::cout<<"p1,channel= "<<w_<<" e= "<<no_ele_p1<<" Hits"<<std::endl;
1095  // Tno_ele_p1+= no_ele_p1;
1096  // sum=Tno_ele_p1;
1097  //
1098  // //std::cout<<"sum= "<<sum<<std::endl;
1099  // }
1100  // //std::cout<<" no_ele_p1= "<< no_ele_p1<<std::endl;
1101  // }
1102  //
1103  //
1104  //
1105  // }//numberofElectrons
1106  // // std::cout<<"TOTAL for p0 is: "<< Tno_ele_p0<<std::endl;
1107  // // std::cout<<"TOTAL for p1 is: "<< Tno_ele_p1<<std::endl;
1108  // sum=0;
1109  // sum0=0;
1110  // }//hits
1111 
1112 
1113  // std::cout<<"***TOTAL for p0 is: "<< Tno_ele_p0<<std::endl;
1114  // std::cout<<"***TOTAL for p1 is: "<< Tno_ele_p1<<std::endl;
1115 
1116 
1117  //-------------------first part FOR BRIAN done---------------------------
1118 
1119  //------------------------------------------------------------
1120  // no_ele_p0=0;
1121  // no_ele_p1=0;
1122  // double Tno_ele_p0_w=0;
1123  // double Tno_ele_p1_w=0;
1124  // sum=0;
1125  // sum0=0;
1126  // std::cout<<"hi"<<std::endl;
1127  // _electrons=0;
1128  // electrons=0;
1129  // // int Total_Elec_p0=0;
1130  // // int Total_Elec_p1=0;
1131  // unsigned int wire=0;
1132  // unsigned int pl=0;
1133  // unsigned int channel=0;
1134  //
1135  // //loop through all wires:
1136  //
1137  // for(art::PtrVector<recob::Wire>::const_iterator wireIter2 = wirelist.begin();
1138  // wireIter2 != wirelist.end(); wireIter2++) // {
1139  // //
1140  //
1141  //
1142  // // std::cout<<"channel: "<<wire<<std::endl;
1143  //
1144  // // loop over the SimChannels to find this one
1145  // art::Ptr<sim::SimChannel> sc;
1146  // for(unsigned int scs = 0; scs < simchans.size(); ++scs)
1147  // if(simchans[scs]->Channel() == channel) sc = simchans[scs];
1148  //
1149  // unsigned int numberOfElectrons = sc->NumberOfElectrons();
1150  //
1151  // // if(numberOfElectrons==0){std::cout<<" ZERO ELEC!!!"<<std::endl;}
1152  // //std::cout<<"# of elec: "<<numberOfElectrons<<" ";
1153  //
1154  // for (size_t i = 0; i!=numberOfElectrons; ++i )
1155  // {
1156  //
1157  // _electrons = sc->GetElectrons(i);
1158  //
1159  //
1160  // if(pl==0) {
1161  //
1162  //
1163  // electrons = sc->GetElectrons(i);
1164  //
1165  // no_ele_p0= electrons-> NumElectrons();
1166  // // std::cout<<"p0,channel= "<<wire<<" e= "<<no_ele_p0<<" Wires"<<std::endl;
1167  //
1168  // //std::cout<<" no_ele_p0= "<< no_ele_p0<<std::endl;
1169  // Tno_ele_p0_w+= no_ele_p0;
1170  // sum0=Tno_ele_p0;
1171  //
1172  // // std::cout<<"sum= "<<sum0<<std::endl;
1173  // }
1174  //
1175  //
1176  //
1177  // // std::cout<<"sum= "<< Tno_ele_p0<<std::endl;
1178  //
1179  //
1180  //
1181  // if(pl==1) {
1182  //
1183  //
1184  // electrons = sc->GetElectrons(i);
1185  // no_ele_p1= electrons-> NumElectrons();
1186  //
1187  // //std::cout<<"p1,channel= "<<wire<<" e= "<<no_ele_p1<<" Wires"<<std::endl;
1188  //
1189  // //std::cout<<" no_ele_p1= "<< no_ele_p1<<std::endl;
1190  // Tno_ele_p1_w+= no_ele_p1;
1191  // sum=Tno_ele_p1;
1192  //
1193  // //std::cout<<"sum= "<<sum<<std::endl;
1194  //
1195  // }
1196  //
1197  //
1198  // }//numberofElectrons
1199  // // std::cout<<"TOTAL for p0 is: "<< Tno_ele_p0_w<<std::endl;
1200  // // std::cout<<"TOTAL for p1 is: "<< Tno_ele_p1_w<<std::endl;
1201  //
1202  // no_ele_p0=0;
1203  // no_ele_p1=0;
1204  //
1205  //
1206  //
1207  //
1208  //
1209  // //--------------------------------------
1210  // // if(pl==0) {
1211  // // Total_Elec_p0 += numberOfElectrons;}
1212  // // if(pl==1) {
1213  // // Total_Elec_p1 += numberOfElectrons;}
1214  // // //------------------------------------
1215  //
1216  // sum=0;
1217  // sum0=0;
1218  // }//loop wires
1219 
1220  // std::cout<<"(wires)TOTAL for p0 is: "<< Tno_ele_p0_w<<std::endl;
1221  // std::cout<<"(wires)TOTAL for p1 is: "<< Tno_ele_p1_w<<std::endl;
1222  // fbrian_in->Fill(Tno_ele_p0_w, Tno_ele_p0 );
1223  // fbrian_coll->Fill(Tno_ele_p1_w, Tno_ele_p1 );
1224 
1225 
1226  }
1227 
1228 
1229 
1230 } //end namespace
1231 
1232 
1233 
1234 
1235 
1236 
1237 namespace cluster{
1238 
1239  DEFINE_ART_MODULE(DBclusterAna)
1240 
1241 }
1242 
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:216
cluster::DBclusterAna::DBclusterAna
DBclusterAna(fhicl::ParameterSet const &pset)
Definition: DBclusterAna_module.cc:129
geo::GeometryCore::Plane
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Definition: GeometryCore.h:2938
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
cheat::BackTrackerService::HitToTrackIDEs
const std::vector< sim::TrackIDE > HitToTrackIDEs(recob::Hit const &hit)
Definition: BackTrackerService_service.cc:145
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Hit.h
Declaration of signal hit object.
cheat::ParticleInventoryService::SetEveIdCalculator
void SetEveIdCalculator(sim::EveIdCalculator *ec)
Definition: ParticleInventoryService.h:56
sim::ParticleList::at
mapped_type at(const key_type &key)
Definition: ParticleList.h:330
std
STL namespace.
RawDigit.h
Definition of basic raw digits.
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:72
cluster::DBclusterAna::analyze
void analyze(const art::Event &evt)
read access to event
Definition: DBclusterAna_module.cc:200
art::Event::isRealData
bool isRealData() const
Definition: Event.h:83
cluster
Cluster finding and building.
Definition: SmallClusterFilter_module.cc:36
Event.h
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:171
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2792
hits
void hits()
Definition: readHits.C:15
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:42
art::PtrVector::push_back
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
part
TString part[npart]
Definition: Style.C:32
sim::ParticleList::begin
iterator begin()
Definition: ParticleList.h:305
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Cluster.h
Declaration of cluster object.
art::PtrVector::size
size_type size() const
Definition: PtrVector.h:308
simb::MCTruth::GetParticle
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:73
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
AssociationUtil.h
Utility object to perform functions of association.
PlaneGeo.h
Encapsulate the construction of a single detector plane.
EmEveIdCalculator.h
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
cheat::BackTrackerService::HitToEveTrackIDEs
const std::vector< sim::TrackIDE > HitToEveTrackIDEs(recob::Hit const &hit)
Definition: BackTrackerService_service.cc:166
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:117
Wire.h
Declaration of basic channel signal object.
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
ParticleList.h
Particle list in DetSim contains Monte Carlo particle information.
art::Event::run
RunNumber_t run() const
Definition: Event.h:77
art::Event::id
EventID id() const
Definition: Event.h:56
Geometry.h
art framework interface to geometry description