203 std::cout <<
"run : " << evt.
run() << std::endl;
205 std::cout <<
"event : " << evt.
id().
event() << std::endl;
211 std::cout<<
"**** DBclusterAna: Bailing. Don't call this module if you're not MC. "<<std::endl;
238 for (
size_t ii = 0; ii < rdListHandle->size(); ++ii){
249 std::vector<int> mc_trackids;
253 for (
auto i = _particleList.
begin(); i != _particleList.
end(); ++i ){
254 int trackID = (*i).first;
255 mc_trackids.push_back(trackID);
270 for (
size_t ii = 0; ii < mctruthListHandle->size(); ++ii){
282 std::vector< art::Ptr<recob::Hit> > hits_vec;
285 std::vector< art::Ptr<recob::Hit> >
hits;
290 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii){
295 std::cout<<
"in Efficiency, clusters.size()= "<<clusters.
size()<<std::endl;
300 for (
size_t ii = 0; ii < wireListHandle->size(); ++ii){
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;
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;
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;
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;
346 double total_Q_cluster_hits=0;
355 if(clusters.
size()!=0 && hits.size()!=0){
356 for(
unsigned int plane=0;plane<geom->
Nplanes();++plane){
359 for(
size_t j = 0; j < clusters.
size(); ++j){
362 if( clusters[j]->View() == view){
365 std::vector< art::Ptr<recob::Hit> > _hits = fmh.at(j);
370 for(
size_t p = 0; p<_hits.size(); ++p){
372 hits_vec.push_back(_hits_ptr);
374 total_Q_cluster_hits += _hits[p]->Integral();
378 std::vector< art::Ptr<recob::Hit> >
::iterator itr = hits_vec.begin();
382 while(itr != hits_vec.end()) {
387 hit_energy=_hits[itr-hits_vec.begin()]->Integral();
389 std::vector<sim::TrackIDE> trackides = bt_serv->
HitToTrackIDEs(*itr);
395 while( idesitr != trackides.end() ){
406 vec_trackid.push_back((*idesitr).trackID);
425 if(pdg==13 || pdg==-13){
427 _en_13+=hit_energy*((*idesitr).energyFrac);
498 for(
unsigned int i=0;i<vec_pdg.size();++i){
519 it=find(vec_pdg.begin(),vec_pdg.end(),13);
520 if(it!=vec_pdg.end()){
528 it2=find(vec_pdg.begin(),vec_pdg.end(),11);
529 if(it2!=vec_pdg.end()){
530 no_cl_for_electron++;
533 it3=find(vec_pdg.begin(),vec_pdg.end(),-11);
534 if(it3!=vec_pdg.end()){
535 no_cl_for_positron++;
538 it4=find(vec_pdg.begin(),vec_pdg.end(),111);
539 if(it4!=vec_pdg.end()){
540 no_cl_for_pion_111++;
542 it6=find(vec_pdg.begin(),vec_pdg.end(),211);
543 if(it6!=vec_pdg.end()){
544 no_cl_for_pion_211++;
546 it7=find(vec_pdg.begin(),vec_pdg.end(),-211);
547 if(it7!=vec_pdg.end()){
548 no_cl_for_pion_m211++;
550 it8=find(vec_pdg.begin(),vec_pdg.end(),2212);
551 if(it8!=vec_pdg.end()){
558 sort( vec_pdg.begin(), vec_pdg.end() );
559 vec_pdg.erase( unique( vec_pdg.begin(), vec_pdg.end() ), vec_pdg.end() );
569 sort( vec_trackid.begin(), vec_trackid.end() );
570 vec_trackid.erase( unique( vec_trackid.begin(), vec_trackid.end() ), vec_trackid.end() );
573 for(
unsigned int ii=0;ii<vec_trackid.size();++ii){
575 all_trackids.push_back(vec_trackid[ii]);
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() );
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() );
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()){
618 ids.push_back(mc_trackids[ii]);
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();
646 vec_trackid_mother.clear();
648 vec_trackid_mother_en.clear();
659 sort( all_trackids.begin(), all_trackids.end() );
660 all_trackids.erase( unique( all_trackids.begin(), all_trackids.end() ), all_trackids.end() );
663 for(
unsigned int ii=0;ii<all_trackids.size();++ii){
674 sort(ids.begin(),ids.end() );
675 ids.erase( unique( ids.begin(), ids.end() ), ids.end() );
677 double no_of_clusters_per_track=noCluster/ids.size();
684 if(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;
719 double result=no_of_particles_in_cluster/no_of_clusters;
723 double no_noise_hits=hits.size()-total_no_hits_in_clusters;
725 double percent_noise=double(no_noise_hits/hits.size())*100;
730 double no_trackid_per_cl_per_event = sum_vec_trackid/no_of_clusters;
748 for(
unsigned int i = 0; i < mclist.
size(); ++i ){
757 for(
int ii = 0; ii < mc->NParticles(); ++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;
764 if(abs(
part.PdgCode()) == 111){std::cout<<
" I have a pi zero!!!"<<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;
803 std::vector< art::Ptr<recob::Hit> >
::iterator itr = hits.begin();
804 while(itr != hits.end()) {
806 std::vector<sim::TrackIDE> trackides = bt_serv->
HitToTrackIDEs(*itr);
811 hit_energy=hits[itr-hits.begin()]->Integral();
813 while( idesitr != trackides.end() ){
839 if(pdg==13 || pdg==-13){hit_13++;
840 en_13+=hit_energy*((*idesitr).energyFrac);}
943 std::cout<<
"****** mu E from clusters = "<<_en_13<<std::endl;
944 std::cout<<
"****** mu E from hits = "<<en_13<<std::endl;
std::string fCalDataModuleLabel
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
TH1F * fPercent_lost_2212_energy
TH1F * fPercent_lost_muon_energy
TH1F * fNoParticles_trackid
const std::vector< sim::TrackIDE > HitToTrackIDEs(recob::Hit const &hit)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
TH1F * fno_of_clusters_per_track
void SetEveIdCalculator(sim::EveIdCalculator *ec)
TH1F * fPercent_lost_111_hits
mapped_type at(const key_type &key)
std::string fGenieGenModuleLabel
TH1F * fNoClustersInEvent
TH1F * fPercent_lost_electron_energy
std::string fDigitModuleLabel
TH1F * fPercent_lost_111_energy
View_t View() const
Which coordinate does this plane measure.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
TH1F * fPercent_lost_2112_energy
TH1F * fPercent_lost_m211_hits
void push_back(Ptr< U > const &p)
TH1F * fPercent_lost_positron_hits
TH1F * fPercent_lost_m211_energy
TH1F * fNoParticles_trackid_per_event
TH1F * fPercent_lost_electron_hits
TH1F * fPercent_lost_muon_hits
const std::vector< sim::TrackIDE > HitToEveTrackIDEs(recob::Hit const &hit)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
TH1F * fNoParticles_pdg_per_event
TH1F * fPercent_lost_positron_energy
EventNumber_t event() const
TH1F * fNoParticles_trackid_mother
const sim::ParticleList & ParticleList()
std::string fHitsModuleLabel
TH1F * fPercent_lost_211_hits
TH1F * fPercent_lost_2112_hits
TH1F * fPercent_lost_211_energy
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TH1F * fPercent_lost_2212_hits
std::string fClusterFinderModuleLabel