| @@ -0,0 +1,162 @@ | ||
| #ifndef HITVIEWER_CXX | ||
| #define HITVIEWER_CXX | ||
|
|
||
| #include "HitViewer.h" | ||
|
|
||
| namespace larlite { | ||
|
|
||
|
|
||
| //******************************** | ||
| HitViewer::HitViewer(): ana_base(), _hHits_U(), _hHits_V(), _hHits_Y() | ||
| //******************************** | ||
| { | ||
| //Class Name | ||
| _name = "HitViewer"; | ||
| //set initialization for pointers | ||
| _hHits_U = 0; | ||
| _hHits_V = 0; | ||
| _hHits_Y = 0; | ||
|
|
||
| _w2cm = larutil::GeometryUtilities::GetME()->WireToCm(); | ||
| _t2cm = larutil::GeometryUtilities::GetME()->TimeToCm(); | ||
|
|
||
| _evtNum = 0; | ||
|
|
||
| _minQ = 0; | ||
|
|
||
| } | ||
|
|
||
| //******************************** | ||
| bool HitViewer::initialize() | ||
| //******************************** | ||
| { | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
|
|
||
| //********************************************** | ||
| bool HitViewer::analyze(storage_manager* storage) | ||
| //********************************************** | ||
| { | ||
|
|
||
|
|
||
|
|
||
| //Get Hits | ||
| auto hits = storage->get_data<event_hit>("gaushit"); | ||
|
|
||
| //clean up histograms if they already exist (from previous event) | ||
| if (_hHits_U) {delete _hHits_U; _hHits_U = 0;}; | ||
| if (_hHits_V) {delete _hHits_V; _hHits_V = 0;}; | ||
| if (_hHits_Y) {delete _hHits_Y; _hHits_Y = 0;}; | ||
|
|
||
| //Define axis ranges | ||
| std::vector<double> chmax, chmin, wiremax, wiremin, timemax, timemin; | ||
| //Find axis boundary | ||
| GetAxisRange(chmax, chmin, wiremax, wiremin, timemax, timemin, hits); | ||
| //proceed only if values actually reset | ||
| if ( wiremax[0] <= -1 ) | ||
| { | ||
| print(msg::kWARNING,__FUNCTION__, | ||
| "Did not find any reconstructed hits in view 0. Skipping this event..."); | ||
| return true; | ||
| } | ||
|
|
||
| //if all ok, plot wire vs. time for hits | ||
| _hHits_U = Prepare2DHisto(Form("Event %i - Hit Charge [ Area in ADCs ] U-Plane",_evtNum), | ||
| wiremin[0]*_w2cm, wiremax[0]*_w2cm, timemin[0]*_t2cm, timemax[0]*_t2cm); | ||
| _hHits_V = Prepare2DHisto(Form("Event %i - Hit Charge [ Area in ADCs ] V-Plane",_evtNum), | ||
| wiremin[1]*_w2cm, wiremax[1]*_w2cm, timemin[1]*_t2cm, timemax[1]*_t2cm); | ||
| _hHits_Y = Prepare2DHisto(Form("Event %i - Hit Charge [ Area in ADCs ] Y-Plane",_evtNum), | ||
| wiremin[2]*_w2cm, wiremax[2]*_w2cm, timemin[2]*_t2cm, timemax[2]*_t2cm); | ||
|
|
||
| //loop over hits | ||
| for (size_t i=0; i<hits->size(); i++) | ||
| { | ||
| const hit *this_hit = (&(hits->at(i))); | ||
| // Place an optional cut on the minimum charge of hits we want to display (settable value) | ||
| if (this_hit->Integral() < _minQ) | ||
| continue; | ||
| //place in right plane | ||
| if ( this_hit->View()==0 ) | ||
| _hHits_U->Fill( this_hit->WireID().Wire*_w2cm, this_hit->PeakTime()*_t2cm, this_hit->Integral() ); | ||
| if ( this_hit->View()==1 ) | ||
| _hHits_V->Fill( this_hit->WireID().Wire*_w2cm, this_hit->PeakTime()*_t2cm, this_hit->Integral() ); | ||
| if ( this_hit->View()==2 ) | ||
| _hHits_Y->Fill( this_hit->WireID().Wire*_w2cm, this_hit->PeakTime()*_t2cm, this_hit->Integral() ); | ||
|
|
||
| }//end loop over hits | ||
|
|
||
| _evtNum +=1; | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| //**************************************************************** | ||
| TH2I* HitViewer::Prepare2DHisto(std::string name, | ||
| double wiremin, double wiremax, | ||
| double timemin, double timemax) | ||
| //**************************************************************** | ||
| { | ||
|
|
||
| TH2I* h=0; | ||
| if(h) delete h; | ||
|
|
||
| h = new TH2I("2DViewer", name.c_str(), | ||
| 200, wiremin-10, wiremax+10, | ||
| 200, timemin-10, timemax+10); | ||
|
|
||
| h->SetXTitle("Wire [cm] "); | ||
| h->SetYTitle("Time [cm]"); | ||
|
|
||
| return h; | ||
| } | ||
|
|
||
| bool HitViewer::finalize() { | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
|
|
||
| //********************************************** | ||
| void HitViewer::GetAxisRange(std::vector<Double_t> &chmax, | ||
| std::vector<Double_t> &chmin, | ||
| std::vector<Double_t> &wiremax, | ||
| std::vector<Double_t> &wiremin, | ||
| std::vector<Double_t> &timemax, | ||
| std::vector<Double_t> &timemin, | ||
| const event_hit* hits) const | ||
| { | ||
|
|
||
| // Make sure input vector is of size wire plane with initial value -1 (if not yet set) | ||
| chmax.resize((larlite::geo::kW+1),-1); | ||
| wiremax.resize((larlite::geo::kW+1),-1); | ||
| timemax.resize((larlite::geo::kW+1),-1); | ||
| chmin.resize((larlite::geo::kW+1),-1); | ||
| wiremin.resize((larlite::geo::kW+1),-1); | ||
| timemin.resize((larlite::geo::kW+1),-1); | ||
|
|
||
| for (size_t i=0; i < hits->size(); i++){ | ||
|
|
||
| const hit *h = (&(hits->at(i))); | ||
|
|
||
| larlite::geo::View_t view = h->View(); | ||
| Double_t wire = (Double_t)(h->WireID().Wire); | ||
| Double_t ch = (Double_t)(h->Channel()); | ||
| Double_t tstart = h->StartTick(); | ||
| Double_t tend = h->EndTick(); | ||
|
|
||
| if( wiremax[view] < 0 || wiremax[view] < wire ) wiremax[view] = wire; | ||
| if( chmax[view] < 0 || chmax[view] < ch ) chmax[view] = ch; | ||
| if( timemax[view] < 0 || timemax[view] < tend ) timemax[view] = tend; | ||
|
|
||
| if( wiremin[view] < 0 || wiremin[view] > wire ) wiremin[view] = wire; | ||
| if( chmin[view] < 0 || chmin[view] > ch ) chmin[view] = ch; | ||
| if( timemin[view] < 0 || timemin[view] > tstart ) timemin[view] = tstart; | ||
| } | ||
|
|
||
| return; | ||
| } | ||
|
|
||
| } | ||
| #endif |
| @@ -0,0 +1,116 @@ | ||
| /** | ||
| * \file HitViewer.h | ||
| * | ||
| * \ingroup Analysis | ||
| * | ||
| * \brief Class def header for a class HitViewer | ||
| * | ||
| * @author David Caratelli | ||
| */ | ||
|
|
||
| /** \addtogroup Analysis | ||
| @{*/ | ||
|
|
||
| #ifndef HITVIEWER_H | ||
| #define HITVIEWER_H | ||
|
|
||
| #include "Analysis/ana_base.h" | ||
| #include "LArUtil/GeometryUtilities.h" | ||
| #include "DataFormat/rawdigit.h" | ||
| #include "DataFormat/hit.h" | ||
| #include <TH2I.h> | ||
| #include <TGraph.h> | ||
| #include <TCanvas.h> | ||
| #include <TPad.h> | ||
| #include <string.h> | ||
|
|
||
| namespace larlite { | ||
| /** | ||
| \class HitViewer | ||
| User custom analysis class made by SHELL_USER_NAME | ||
| */ | ||
| class HitViewer : public ana_base{ | ||
|
|
||
| public: | ||
|
|
||
| /// Default constructor | ||
| HitViewer(); | ||
|
|
||
| /// Default destructor | ||
| virtual ~HitViewer(){}; | ||
|
|
||
| /** IMPLEMENT in HitViewer.cc! | ||
| Initialization method to be called before the analysis event loop. | ||
| */ | ||
| virtual bool initialize(); | ||
|
|
||
| /** IMPLEMENT in HitViewer.cc! | ||
| Analyze a data event-by-event | ||
| */ | ||
| virtual bool analyze(storage_manager* storage); | ||
|
|
||
| /** IMPLEMENT in HitViewer.cc! | ||
| Finalize method to be called after all events processed. | ||
| */ | ||
| virtual bool finalize(); | ||
|
|
||
| /// set minimum hit charge for it to be visible | ||
| void setMinQ(double q) { _minQ = q; } | ||
|
|
||
| /// A utility function to (re)create Th3D histogram of a specified boundary & name | ||
| TH2I* Prepare2DHisto(std::string name, | ||
| double wiremin, double wiremax, | ||
| double timemin, double timemax); | ||
|
|
||
| TGraph* PrepareGraph(); | ||
|
|
||
| /// Getter for hit TH2I histo, weighted by charge | ||
| const TH2I* GetHisto_Hits (int view) const { | ||
| if(view==0) | ||
| return _hHits_U; | ||
| else if(view==1) | ||
| return _hHits_V; | ||
| else if(view==2) | ||
| return _hHits_Y; | ||
| else { | ||
| std::cout<<"*******************you screwed something up. view should be 0 1 or 2!"<<std::endl; | ||
| std::cout<<"returning _hHits_U because i don't know what else to return"<<std::endl; | ||
| return _hHits_U; | ||
| } | ||
| }; | ||
|
|
||
| protected: | ||
|
|
||
| /// Counter for event Number | ||
| int _evtNum; | ||
|
|
||
| /// Main canvas | ||
| TCanvas* _c1; | ||
| /// Main pad | ||
| TPad* _p1; | ||
|
|
||
| /// Hit histograms to sit next to cluster ones, | ||
| TH2I* _hHits_U; | ||
| TH2I* _hHits_V; | ||
| TH2I* _hHits_Y; | ||
|
|
||
| geo::View_t iview; | ||
| std::vector<hit> ihit_v; | ||
|
|
||
| double _w2cm; | ||
| double _t2cm; | ||
|
|
||
| /// Set minimum hit Integarl for hit to be viewable | ||
| double _minQ; | ||
|
|
||
| void GetAxisRange(std::vector<Double_t> &chmax,std::vector<Double_t> &chmin, | ||
| std::vector<Double_t> &wiremax, std::vector<Double_t> &wiremin, | ||
| std::vector<Double_t> &timemax, std::vector<Double_t> &timemin, | ||
| const event_hit* hits) const; | ||
|
|
||
| }; | ||
| } | ||
| #endif | ||
|
|
||
| /** @} */ // end of doxygen group |
| @@ -0,0 +1,13 @@ | ||
| #ifndef POSITION_CXX | ||
| #define POSITION_CXX | ||
|
|
||
| #include "Position.h" | ||
|
|
||
| namespace radius{ | ||
|
|
||
| //Point::Point() | ||
| //{ this->clear(); } | ||
|
|
||
| } | ||
|
|
||
| #endif |
| @@ -0,0 +1,35 @@ | ||
| #ifndef POSITION_H | ||
| #define POSITION_H | ||
|
|
||
| #include <cmath> | ||
| #include <vector> | ||
| #include <iostream> | ||
|
|
||
| namespace radius{ | ||
|
|
||
| class Point : public std::vector<double> { | ||
|
|
||
| public: | ||
|
|
||
| Point() {}; | ||
|
|
||
| Point(const std::vector<double> &v) : std::vector<double>(v) | ||
| {} | ||
|
|
||
| virtual ~Point(){}; | ||
|
|
||
| bool operator< (const Point& p1) | ||
| { | ||
| for (size_t i=0; i < 3; i++){ | ||
| if ((*this)[i] < p1[i]) return true; | ||
| if ((*this)[i] > p1[i]) return false; | ||
| } | ||
| return false; | ||
| } | ||
|
|
||
| private: | ||
|
|
||
| }; | ||
| } | ||
|
|
||
| #endif |
| @@ -0,0 +1,198 @@ | ||
| #ifndef LARLITE_SHOWERQUALITYSIMPLE_CXX | ||
| #define LARLITE_SHOWERQUALITYSIMPLE_CXX | ||
|
|
||
| #include "ShowerQualitySimple.h" | ||
| #include "DataFormat/shower.h" | ||
| #include "DataFormat/mcshower.h" | ||
| #include <TTree.h> | ||
|
|
||
| namespace larlite { | ||
|
|
||
| ShowerQualitySimple::ShowerQualitySimple() | ||
| : _mc_tree(nullptr) | ||
| , _reco_tree(nullptr) | ||
| { | ||
| _fout = 0; | ||
| _name = "ShowerQualitySimple"; | ||
| _producer_name = "showerreco"; | ||
| _minE = 100; | ||
| } | ||
|
|
||
| bool ShowerQualitySimple::initialize() { | ||
|
|
||
| _evt = 0; | ||
|
|
||
| if (_mc_tree) { delete _mc_tree; } | ||
| _mc_tree = new TTree("_mc_tree","MC Shower Info Tree"); | ||
| _mc_tree->Branch("_evt",&_evt,"evt/I"); | ||
| _mc_tree->Branch("_pdg_mc",&_pdg_mc,"pdg_mc/I"); | ||
| _mc_tree->Branch("_e_det_mc",&_e_det_mc,"e_det_mc/D"); | ||
| _mc_tree->Branch("_x_det_mc",&_x_det_mc,"x_det_mc/D"); | ||
| _mc_tree->Branch("_y_det_mc",&_y_det_mc,"y_det_mc/D"); | ||
| _mc_tree->Branch("_z_det_mc",&_z_det_mc,"z_det_mc/D"); | ||
| _mc_tree->Branch("_px_det_mc",&_px_det_mc,"px_det_mc/D"); | ||
| _mc_tree->Branch("_py_det_mc",&_py_det_mc,"py_det_mc/D"); | ||
| _mc_tree->Branch("_pz_det_mc",&_pz_det_mc,"pz_det_mc/D"); | ||
| _mc_tree->Branch("_e_sim_mc",&_e_sim_mc,"e_sim_mc/D"); | ||
| _mc_tree->Branch("_x_sim_mc",&_x_sim_mc,"x_sim_mc/D"); | ||
| _mc_tree->Branch("_y_sim_mc",&_y_sim_mc,"y_sim_mc/D"); | ||
| _mc_tree->Branch("_z_sim_mc",&_z_sim_mc,"z_sim_mc/D"); | ||
| _mc_tree->Branch("_px_sim_mc",&_px_sim_mc,"px_sim_mc/D"); | ||
| _mc_tree->Branch("_py_sim_mc",&_py_sim_mc,"py_sim_mc/D"); | ||
| _mc_tree->Branch("_pz_sim_mc",&_pz_sim_mc,"pz_sim_mc/D"); | ||
| _mc_tree->Branch("_containment",&_containment,"containment/D"); | ||
|
|
||
| if (_reco_tree) { delete _reco_tree; } | ||
| _reco_tree = new TTree("_reco_tree","RECO Shower Info Tree"); | ||
| _reco_tree->Branch("_evt",&_evt,"evt/I"); | ||
| _reco_tree->Branch("_best_mcs",&_best_mcs,"best_mcs/I"); | ||
| _reco_tree->Branch("_dedx",&_dedx,"dedx/D"); | ||
| _reco_tree->Branch("_dist",&_dist,"dist/D"); | ||
| _reco_tree->Branch("_dir",&_dir,"dir/D"); | ||
| _reco_tree->Branch("_e_reco",&_e_reco,"e_reco/D"); | ||
| _reco_tree->Branch("_x_reco",&_x_reco,"x_reco/D"); | ||
| _reco_tree->Branch("_y_reco",&_y_reco,"y_reco/D"); | ||
| _reco_tree->Branch("_z_reco",&_z_reco,"z_reco/D"); | ||
| _reco_tree->Branch("_px_reco",&_px_reco,"px_reco/D"); | ||
| _reco_tree->Branch("_py_reco",&_py_reco,"py_reco/D"); | ||
| _reco_tree->Branch("_pz_reco",&_pz_reco,"pz_reco/D"); | ||
| _reco_tree->Branch("_pdg_best",&_pdg_best,"pdg_best/I"); | ||
| _reco_tree->Branch("_e_det_best",&_e_det_best,"e_det_best/D"); | ||
| _reco_tree->Branch("_x_det_best",&_x_det_best,"x_det_best/D"); | ||
| _reco_tree->Branch("_y_det_best",&_y_det_best,"y_det_best/D"); | ||
| _reco_tree->Branch("_z_det_best",&_z_det_best,"z_det_best/D"); | ||
| _reco_tree->Branch("_px_det_best",&_px_det_best,"px_det_best/D"); | ||
| _reco_tree->Branch("_py_det_best",&_py_det_best,"py_det_best/D"); | ||
| _reco_tree->Branch("_pz_det_best",&_pz_det_best,"pz_det_best/D"); | ||
| _reco_tree->Branch("_e_sim_best",&_e_sim_best,"e_sim_best/D"); | ||
| _reco_tree->Branch("_x_sim_best",&_x_sim_best,"x_sim_best/D"); | ||
| _reco_tree->Branch("_y_sim_best",&_y_sim_best,"y_sim_best/D"); | ||
| _reco_tree->Branch("_z_sim_best",&_z_sim_best,"z_sim_best/D"); | ||
| _reco_tree->Branch("_px_sim_best",&_px_sim_best,"px_sim_best/D"); | ||
| _reco_tree->Branch("_py_sim_best",&_py_sim_best,"py_sim_best/D"); | ||
| _reco_tree->Branch("_pz_sim_best",&_pz_sim_best,"pz_sim_best/D"); | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| bool ShowerQualitySimple::analyze(storage_manager* storage) { | ||
|
|
||
| _evt += 1; | ||
|
|
||
| // Retrieve mcshower data product | ||
| auto ev_mcs = storage->get_data<event_mcshower>("mcreco"); | ||
|
|
||
| if(!ev_mcs || !(ev_mcs->size())) { | ||
| print(msg::kERROR,__FUNCTION__,"MCShower data product not found!"); | ||
| return false; | ||
| } | ||
|
|
||
| // Get shower | ||
| auto ev_shr = storage->get_data<event_shower>(_producer_name); | ||
| if(!ev_shr) { | ||
| print(msg::kERROR,__FUNCTION__, | ||
| Form("Did not find shower produced by \"%s\"",_producer_name.c_str()) | ||
| ); | ||
| return false; | ||
| } | ||
|
|
||
| // Loop through mcshowers and recoshowers | ||
| for (size_t i=0; i < ev_mcs->size(); i++){ | ||
| if (ev_mcs->at(i).DetProfile().E() < _minE) { continue; } | ||
| auto const& mcs = ev_mcs->at(i); | ||
| auto const& det = mcs.DetProfile(); | ||
| auto const& sim = mcs.Start(); | ||
| _pdg_mc = mcs.PdgCode(); | ||
| _e_det_mc = det.E(); | ||
| _x_det_mc = det.X(); | ||
| _y_det_mc = det.Y(); | ||
| _z_det_mc = det.Z(); | ||
| _px_det_mc = det.Px(); | ||
| _py_det_mc = det.Py(); | ||
| _pz_det_mc = det.Pz(); | ||
| _e_sim_mc = sim.E(); | ||
| _x_sim_mc = sim.X(); | ||
| _y_sim_mc = sim.Y(); | ||
| _z_sim_mc = sim.Z(); | ||
| _px_sim_mc = sim.Px(); | ||
| _py_sim_mc = sim.Py(); | ||
| _pz_sim_mc = sim.Pz(); | ||
| _containment = _e_det_mc/_e_sim_mc; | ||
| _mc_tree->Fill(); | ||
| } | ||
|
|
||
| // loop through reco showers | ||
| for (size_t j=0; j < ev_shr->size(); j++){ | ||
|
|
||
| auto const& shr = ev_shr->at(j); | ||
| auto bestplane = shr.best_plane(); | ||
| _dedx = shr.dEdx()[bestplane]; | ||
| _e_reco = shr.Energy()[bestplane]; | ||
| _x_reco = shr.ShowerStart()[0]; | ||
| _y_reco = shr.ShowerStart()[1]; | ||
| _z_reco = shr.ShowerStart()[2]; | ||
| _px_reco = shr.Direction()[0]; | ||
| _py_reco = shr.Direction()[1]; | ||
| _pz_reco = shr.Direction()[2]; | ||
|
|
||
| // loop through all mc again and compare -> try and find the best one | ||
| _best_mcs = -1; | ||
| double best_dist = 1036; | ||
| for (size_t i=0; i < ev_mcs->size(); i++){ | ||
| if (ev_mcs->at(i).DetProfile().E() < _minE) { continue; } | ||
| auto const& mcs = ev_mcs->at(i); | ||
| auto const& det = mcs.DetProfile(); | ||
| auto const& sim = mcs.Start(); | ||
| double dist = sqrt ( (det.X()-_x_reco)*(det.X()-_x_reco) + | ||
| (det.X()-_x_reco)*(det.X()-_x_reco) + | ||
| (det.X()-_x_reco)*(det.X()-_x_reco) ); | ||
| if (dist < best_dist){ | ||
| best_dist = dist; | ||
| _best_mcs = i; | ||
| } | ||
| } | ||
| if (_best_mcs >= 0){ | ||
| auto const& mcs = ev_mcs->at(_best_mcs); | ||
| auto const& det = mcs.DetProfile(); | ||
| auto const& sim = mcs.Start(); | ||
| _pdg_best = mcs.PdgCode(); | ||
| _e_det_best = det.E(); | ||
| _x_det_best = det.X(); | ||
| _y_det_best = det.Y(); | ||
| _z_det_best = det.Z(); | ||
| _px_det_best = det.Px(); | ||
| _py_det_best = det.Py(); | ||
| _pz_det_best = det.Pz(); | ||
| _e_sim_best = sim.E(); | ||
| _x_sim_best = sim.X(); | ||
| _y_sim_best = sim.Y(); | ||
| _z_sim_best = sim.Z(); | ||
| _px_sim_best = sim.Px(); | ||
| _py_sim_best = sim.Py(); | ||
| _pz_sim_best = sim.Pz(); | ||
| _dist = sqrt ( (det.X()-_x_reco)*(det.X()-_x_reco) + | ||
| (det.X()-_x_reco)*(det.X()-_x_reco) + | ||
| (det.X()-_x_reco)*(det.X()-_x_reco) ); | ||
| _dir = ( ( (det.Px()*_px_reco) + (det.Py()*_py_reco) + (det.Pz()*_pz_reco) ) / | ||
| sqrt ( (_px_reco*_px_reco) + (_py_reco*_py_reco) + (_pz_reco*_pz_reco) ) ); | ||
| }// if we found a best shower | ||
| _reco_tree->Fill(); | ||
|
|
||
| } | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| bool ShowerQualitySimple::finalize() { | ||
|
|
||
| if (_fout){ | ||
| if (_mc_tree) { _mc_tree->Write(); } | ||
| if (_reco_tree) { _reco_tree->Write(); } | ||
| } | ||
|
|
||
|
|
||
| return true; | ||
| } | ||
|
|
||
| } | ||
| #endif |
| @@ -0,0 +1,114 @@ | ||
| /** | ||
| * \file ShowerQualitySimple.h | ||
| * | ||
| * \ingroup ShowerReco3D | ||
| * | ||
| * \brief Class def header for a class ShowerQualitySimple | ||
| * | ||
| * @author david caratelli | ||
| */ | ||
|
|
||
| /** \addtogroup ShowerReco3D | ||
| @{*/ | ||
|
|
||
| #ifndef LARLITE_SHOWERQUALITYSIMPLE_H | ||
| #define LARLITE_SHOWERQUALITYSIMPLE_H | ||
|
|
||
| #include "Analysis/ana_base.h" | ||
|
|
||
| namespace larlite { | ||
| /** | ||
| \class ShowerQualitySimple | ||
| User custom analysis class made by SHELL_USER_NAME | ||
| */ | ||
| class ShowerQualitySimple : public ana_base{ | ||
|
|
||
| public: | ||
|
|
||
| /// Default constructor | ||
| ShowerQualitySimple(); | ||
|
|
||
| /// Default destructor | ||
| virtual ~ShowerQualitySimple(){} | ||
|
|
||
| virtual bool initialize(); | ||
|
|
||
| virtual bool analyze(storage_manager* storage); | ||
|
|
||
| virtual bool finalize(); | ||
|
|
||
| /// Set producer name | ||
| void setProducerName(std::string name) { _producer_name = name; } | ||
|
|
||
| /// Set minimum energy to use | ||
| void setMinE(double e) { _minE = e; } | ||
|
|
||
| protected: | ||
|
|
||
| std::string _producer_name; | ||
|
|
||
| int _evt; | ||
|
|
||
| double _minE; | ||
|
|
||
| TTree* _mc_tree; | ||
| int _pdg_mc; | ||
| double _e_det_mc; | ||
| double _x_det_mc; | ||
| double _y_det_mc; | ||
| double _z_det_mc; | ||
| double _px_det_mc; | ||
| double _py_det_mc; | ||
| double _pz_det_mc; | ||
| double _e_sim_mc; | ||
| double _x_sim_mc; | ||
| double _y_sim_mc; | ||
| double _z_sim_mc; | ||
| double _px_sim_mc; | ||
| double _py_sim_mc; | ||
| double _pz_sim_mc; | ||
| double _containment; | ||
|
|
||
|
|
||
| TTree* _reco_tree; | ||
| double _dedx; | ||
| double _e_reco; | ||
| double _x_reco; | ||
| double _y_reco; | ||
| double _z_reco; | ||
| double _px_reco; | ||
| double _py_reco; | ||
| double _pz_reco; | ||
| int _pdg_best; | ||
| int _best_mcs; | ||
| double _e_det_best; | ||
| double _x_det_best; | ||
| double _y_det_best; | ||
| double _z_det_best; | ||
| double _px_det_best; | ||
| double _py_det_best; | ||
| double _pz_det_best; | ||
| double _e_sim_best; | ||
| double _x_sim_best; | ||
| double _y_sim_best; | ||
| double _z_sim_best; | ||
| double _px_sim_best; | ||
| double _py_sim_best; | ||
| double _pz_sim_best; | ||
| double _dir; | ||
| double _dist; | ||
|
|
||
| }; | ||
| } | ||
| #endif | ||
|
|
||
| //************************************************************************** | ||
| // | ||
| // For Analysis framework documentation, read Manual.pdf here: | ||
| // | ||
| // http://microboone-docdb.fnal.gov:8080/cgi-bin/ShowDocument?docid=3183 | ||
| // | ||
| //************************************************************************** | ||
|
|
||
| /** @} */ // end of doxygen group |
| @@ -0,0 +1,187 @@ | ||
| #ifndef LARLITE_TRACKHITREMOVER_CXX | ||
| #define LARLITE_TRACKHITREMOVER_CXX | ||
|
|
||
| #include "TrackHitRemover.h" | ||
|
|
||
| namespace larlite { | ||
|
|
||
| TrackHitRemover::TrackHitRemover(){ | ||
| _name = "TrkHitRemover"; | ||
| _fout = 0; | ||
| _trackProducer = ""; | ||
| _showerProducer = ""; | ||
| _outHitProducer = "outhits"; | ||
| _verbose = false; | ||
| _getShowerHits = false; | ||
| } | ||
|
|
||
| bool TrackHitRemover::initialize() { | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| bool TrackHitRemover::analyze(storage_manager* storage) { | ||
|
|
||
| // clear vectors to be filled with hit indices: | ||
| _trk_hit_indices.clear(); | ||
| _shr_hit_indices.clear(); | ||
| _out_hits.clear(); | ||
|
|
||
|
|
||
| // ************************************ | ||
| // Now get hits associated with showers | ||
| if (_getShowerHits){ | ||
| std::cout << "looking for shower hits" << std::endl; | ||
| // Get shower | ||
| auto ev_shower = storage->get_data<event_shower>(_showerProducer); | ||
| if(!ev_shower) { | ||
| print(msg::kERROR,__FUNCTION__,Form("Did not find shower produced by \"%s\"",_showerProducer.c_str())); | ||
| return false; | ||
| } | ||
|
|
||
| std::cout << "there are " << ev_shower->size() << " showers" << std::endl; | ||
|
|
||
| // get associated clusters | ||
| event_cluster* ev_cluster = nullptr; | ||
| auto const& ass_cluster_v = storage->find_one_ass(ev_shower->id(),ev_cluster,ev_shower->name()); | ||
|
|
||
| if (!ev_cluster){ | ||
| print(msg::kERROR,__FUNCTION__, Form("No associated cluster found to a shower produced by \"%s\"", _showerProducer.c_str())); | ||
| return false; | ||
| } | ||
| else if(ev_cluster->size()<1) { | ||
| print(msg::kERROR,__FUNCTION__,Form("There are 0 clusters in this event! Skipping......")); | ||
| return false; | ||
| } | ||
|
|
||
| // get associated hits | ||
| event_hit* ev_hit_shr = nullptr; | ||
| auto const& ass_hit_v = storage->find_one_ass(ev_cluster->id(),ev_hit_shr,ev_cluster->name()); | ||
|
|
||
| if (!ev_hit_shr){ | ||
| print(msg::kERROR,__FUNCTION__, Form("No associated hit found to a shower produced by \"%s\"", ev_cluster->name().c_str())); | ||
| return false; | ||
| } | ||
| else if(ev_hit_shr->size()<1) { | ||
| print(msg::kERROR,__FUNCTION__,Form("There are 0 hits in this event! Skipping......")); | ||
| return false; | ||
| } | ||
|
|
||
| _event_hit_shr = *ev_hit_shr; | ||
|
|
||
| // fill vector of shower-hit index | ||
| std::cout << "number of associated clusters: " << ass_cluster_v.size() << std::endl; | ||
| for (size_t i=0; i < ass_cluster_v.size(); i++){ | ||
| std::vector<unsigned int> shrHits; | ||
| for (size_t j=0; j < ass_cluster_v[i].size(); j++){ | ||
| for (size_t k=0; k < ass_hit_v[ass_cluster_v[i][j]].size(); k++) | ||
| shrHits.push_back(ass_hit_v[ass_cluster_v[i][j]][k]); | ||
| } | ||
| _shr_hit_indices.push_back(shrHits); | ||
| } | ||
|
|
||
| }// if we should get hits associated with showers | ||
| // ********************************************** | ||
|
|
||
| // If we did not specify a track producer | ||
| // just move on! | ||
| if (_trackProducer == "") | ||
| return true; | ||
|
|
||
| // Get Tracks that we want to remove | ||
| auto ev_track = storage->get_data<event_track>(_trackProducer); | ||
| if(!ev_track) { | ||
| print(msg::kERROR,__FUNCTION__,Form("Did not find shower produced by \"%s\"",_trackProducer.c_str())); | ||
| return false; | ||
| } | ||
|
|
||
| // get associated hits | ||
| event_hit* ev_hit_trk = nullptr; | ||
| auto const& ass_hit_v = storage->find_one_ass(ev_track->id(),ev_hit_trk,ev_track->name()); | ||
|
|
||
| _event_hit_trk = *ev_hit_trk; | ||
| _ass_v = ass_hit_v; | ||
|
|
||
| if (!ev_hit_trk){ | ||
| print(msg::kERROR,__FUNCTION__, Form("No associated hit found to a track produced by \"%s\"", _trackProducer.c_str())); | ||
| return false; | ||
| } | ||
| else if(ev_hit_trk->size()<1) { | ||
| print(msg::kERROR,__FUNCTION__,Form("There are 0 hits in this event! Skipping......")); | ||
| return false; | ||
| } | ||
|
|
||
|
|
||
| // loop over associated hit vector to find hits associated with tracks | ||
| int num_ass_hits = 0; | ||
| for (size_t t=0; t < ass_hit_v.size(); t++){ | ||
| if (_verbose) { std::cout << "track " << t << " has " << ass_hit_v[t].size() << " hits associated with it" << std::endl; } | ||
| num_ass_hits += ass_hit_v[t].size(); | ||
| } | ||
|
|
||
| // Get all hits of this type | ||
| event_hit* all_hit = storage->get_data<event_hit>(ev_hit_trk->name()); | ||
|
|
||
| if (_verbose){ | ||
| std::cout << "There are " << ass_hit_v.size() << " tracks in this event" << std::endl; | ||
| std::cout << "Hits associated with track of type " << _trackProducer.c_str() << ": " << num_ass_hits << std::endl; | ||
| std::cout << "Hits of type with track of type " << ev_hit_trk->name().c_str() << ": " << all_hit->size() << std::endl; | ||
| std::cout << std::endl; | ||
| } | ||
|
|
||
| // vector containing all indices of all hits to be removed | ||
| std::vector<unsigned int> hits_to_remove_indices; | ||
|
|
||
| // For the hits in each track, make a Cluster and calculate its parameters using CPAN | ||
| for (size_t t=0; t < ass_hit_v.size(); t++){ | ||
| std::vector<const larlite::hit*> hits; | ||
| for (size_t h=0; h < ass_hit_v[t].size(); h++) | ||
| hits.push_back(&(ev_hit_trk->at(ass_hit_v[t][h]))); | ||
| ::cluster::ClusterParamsAlg cluster(hits); | ||
| cluster.FillParams(true,true,true,true,true,false); | ||
| // get values of quantities to be used | ||
| auto const& mod_hit_dens = cluster.GetParams().modified_hit_density; | ||
| auto const& multi_hit_wires = cluster.GetParams().multi_hit_wires; | ||
| auto const& principal = TMath::Log(1-cluster.GetParams().eigenvalue_principal); | ||
| auto const& opening_angle = cluster.GetParams().opening_angle; | ||
|
|
||
| // if ( (mod_hit_dens < 1.4) || (multi_hit_wires < 3.5) || (principal < -6) || (hits.size() > 200) ){ | ||
| if ( (hits.size() > 200) || ( opening_angle < 0.5 && hits.size() > 50) ) { | ||
| // then we have a track | ||
| std::vector<unsigned int> trkHitIndices; | ||
| for (size_t h=0; h < ass_hit_v[t].size(); h++){ | ||
| trkHitIndices.push_back(ass_hit_v[t][h]); // add hits to list of hits to remove | ||
| hits_to_remove_indices.push_back(ass_hit_v[t][h]); | ||
| } | ||
| _trk_hit_indices.push_back(trkHitIndices); | ||
| }// if cluster is track-like | ||
| } | ||
|
|
||
| // iterator to search list of removed hits | ||
| std::vector<unsigned int>::iterator it; | ||
| // create new event_hit object where to store output hits | ||
| auto ev_hits_out = storage->get_data<event_hit>(_outHitProducer); | ||
| for (size_t n=0; n < ev_hit_trk->size(); n++){ | ||
| // find this hit | ||
| it = std::find(hits_to_remove_indices.begin(), hits_to_remove_indices.end(), n); | ||
| // if this hit is not found in the list to be removed | ||
| if (it == hits_to_remove_indices.end()){ | ||
| _out_hits.push_back(n); | ||
| ev_hits_out->push_back(ev_hit_trk->at(n)); | ||
| } | ||
| }// for all hits | ||
|
|
||
| std::cout << "input hits: " << ev_hit_trk->size() << std::endl; | ||
| std::cout << "output hits: " << ev_hits_out->size() << std::endl; | ||
|
|
||
|
|
||
| return true; | ||
| } | ||
|
|
||
| bool TrackHitRemover::finalize() { | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| } | ||
| #endif |
| @@ -0,0 +1,130 @@ | ||
| /** | ||
| * \file TrackHitRemover.h | ||
| * | ||
| * \ingroup Playground | ||
| * | ||
| * \brief Class def header for a class TrackHitRemover | ||
| * | ||
| * @author david caratelli | ||
| */ | ||
|
|
||
| /** \addtogroup Playground | ||
| @{*/ | ||
|
|
||
| #ifndef LARLITE_TRACKHITREMOVER_H | ||
| #define LARLITE_TRACKHITREMOVER_H | ||
|
|
||
| #include "Analysis/ana_base.h" | ||
| #include "DataFormat/track.h" | ||
| #include "DataFormat/shower.h" | ||
| #include "DataFormat/cluster.h" | ||
| #include "DataFormat/hit.h" | ||
| #include "ClusterParamsAlg.h" | ||
| #include <algorithm> // std::find | ||
|
|
||
| namespace larlite { | ||
| /** | ||
| \class TrackHitRemover | ||
| User custom analysis class made by SHELL_USER_NAME | ||
| */ | ||
| class TrackHitRemover : public ana_base{ | ||
|
|
||
| public: | ||
|
|
||
| /// Default constructor | ||
| TrackHitRemover(); | ||
|
|
||
| /// Default destructor | ||
| virtual ~TrackHitRemover(){} | ||
|
|
||
| virtual bool initialize(); | ||
|
|
||
| virtual bool analyze(storage_manager* storage); | ||
|
|
||
| virtual bool finalize(); | ||
|
|
||
| /// Set producer name for tracks | ||
| void setTrackProducer(std::string s) { _trackProducer = s; } | ||
|
|
||
| /// Set producer name for showers | ||
| void setShowerProducer(std::string s) { _showerProducer = s; } | ||
|
|
||
| /// Set producer name for output hits | ||
| void setOutHitProducer(std::string s) { _outHitProducer = s; } | ||
|
|
||
| /// Hit Getter (for Track Hits) | ||
| ::larlite::event_hit getTrkHits() { return _event_hit_trk; } | ||
|
|
||
| /// Hit Getter (for Shower Hits) | ||
| ::larlite::event_hit getShrHits() { return _event_hit_shr; } | ||
|
|
||
| /// Ass Getter | ||
| std::vector<std::vector<unsigned int> > getAss() { return _ass_v; } | ||
|
|
||
| /// getter for vector of track-like hit indices | ||
| std::vector<std::vector<unsigned int> > getTrkHitIndices() { return _trk_hit_indices; } | ||
|
|
||
| /// getter for vector of shower-like hit indices | ||
| std::vector<std::vector<unsigned int> > getShrHitIndices() { return _shr_hit_indices; } | ||
|
|
||
| /// getter for vector of output hit indices | ||
| std::vector<unsigned int> getOutHits() { return _out_hits; } | ||
|
|
||
| /// Set Verbosity | ||
| void setVerbose(bool on) { _verbose = on; } | ||
|
|
||
| /// Set whether to save hits associated with showers | ||
| void setGetShowerHits(bool on) { _getShowerHits = on; } | ||
|
|
||
| protected: | ||
|
|
||
| // verbosity boolean | ||
| bool _verbose; | ||
|
|
||
| // boolean to decide if to get hits ass. with showers | ||
| bool _getShowerHits; | ||
|
|
||
| // track producer name | ||
| std::string _trackProducer; | ||
|
|
||
| // shower producer name | ||
| std::string _showerProducer; | ||
|
|
||
| // output hit producer name | ||
| std::string _outHitProducer; | ||
|
|
||
| // hit vector (for track hits) | ||
| ::larlite::event_hit _event_hit_trk; | ||
|
|
||
| // hit vector (for shower hits) | ||
| ::larlite::event_hit _event_hit_shr; | ||
|
|
||
| // output hit index vector | ||
| std::vector<unsigned int> _out_hits; | ||
|
|
||
| // associated hit vector | ||
| std::vector<std::vector<unsigned int> > _ass_v; | ||
|
|
||
| /// vector of track-like hit indices | ||
| std::vector<std::vector<unsigned int> > _trk_hit_indices; | ||
|
|
||
| /// vector of shower-like hit indices | ||
| std::vector<std::vector<unsigned int> > _shr_hit_indices; | ||
|
|
||
| // ClusterParamsAlg instance | ||
| ::cluster::ClusterParamsAlg _clusterAlg; | ||
|
|
||
| }; | ||
| } | ||
| #endif | ||
|
|
||
| //************************************************************************** | ||
| // | ||
| // For Analysis framework documentation, read Manual.pdf here: | ||
| // | ||
| // http://microboone-docdb.fnal.gov:8080/cgi-bin/ShowDocument?docid=3183 | ||
| // | ||
| //************************************************************************** | ||
|
|
||
| /** @} */ // end of doxygen group |
| @@ -0,0 +1,256 @@ | ||
| #ifndef LARLITE_VERTEXFINDER_CXX | ||
| #define LARLITE_VERTEXFINDER_CXX | ||
|
|
||
| #include "VertexFinder.h" | ||
| #include "DataFormat/hit.h" | ||
| #include "DataFormat/mctruth.h" | ||
| #include "LArUtil/GeometryUtilities.h" | ||
| #include "LArUtil/Geometry.h" | ||
| #include <math.h> | ||
|
|
||
| namespace larlite { | ||
|
|
||
| bool VertexFinder::initialize() { | ||
|
|
||
| _wire2cm = larutil::GeometryUtilities::GetME()->WireToCm(); | ||
| _time2cm = larutil::GeometryUtilities::GetME()->TimeToCm(); | ||
| _detwidth = 2*larutil::Geometry::GetME()->DetHalfWidth(); | ||
|
|
||
| if (!_vtx_tree) { _vtx_tree = new TTree("_vtx_tree","Vertex Tree"); } | ||
| _vtx_tree->Branch("_distY",&_distY,"distY/D"); | ||
| _vtx_tree->Branch("_distV",&_distV,"distV/D"); | ||
| _vtx_tree->Branch("_distU",&_distU,"distU/D"); | ||
| _vtx_tree->Branch("_distYZ",&_distYZ,"distYZ/D"); | ||
| _vtx_tree->Branch("_distX",&_distX,"distX/D"); | ||
| _vtx_tree->Branch("_dist3D",&_dist3D,"dist3D/D"); | ||
| _vtx_tree->Branch("_area",&_area,"area/D"); | ||
| _vtx_tree->Branch("_xAvg",&_xAvg,"xAvg/D"); | ||
| _vtx_tree->Branch("_yAvg",&_yAvg,"yAvg/D"); | ||
| _vtx_tree->Branch("_zAvg",&_zAvg,"zAvg/D"); | ||
| _vtx_tree->Branch("_nu_x",&_nu_x,"nu_x/D"); | ||
| _vtx_tree->Branch("_nu_y",&_nu_y,"nu_y/D"); | ||
| _vtx_tree->Branch("_nu_z",&_nu_z,"nu_z/D"); | ||
| _vtx_tree->Branch("_QvsTheta","std::vector<double>",&_QvsTheta); | ||
| _vtx_tree->Branch("_nProtons",&_nProtons,"nProtons/I"); | ||
| _vtx_tree->Branch("_evt",&_evt,"evt/I"); | ||
|
|
||
| _evt = 0; | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| bool VertexFinder::analyze(storage_manager* storage) { | ||
|
|
||
| auto evt_hits = storage->get_data<event_hit>("gaushit"); | ||
|
|
||
| if (!evt_hits){ | ||
| std::cout << "No hits!" << std::endl; | ||
| return false; | ||
| } | ||
|
|
||
| // get neutrino vtx info | ||
| auto evt_truth = storage->get_data<event_mctruth>("generator"); | ||
| if (!evt_truth or !evt_truth->size()){ | ||
| std::cout << "No mctruth!" << std::endl; | ||
| return false; | ||
| } | ||
| auto const& truth = evt_truth->at(0); | ||
| auto const& nu = truth.GetNeutrino(); | ||
| auto const& firstStep = nu.Nu().Trajectory().back(); | ||
| _nu_z = firstStep.Z(); | ||
| _nu_x = firstStep.X()+_detwidth; | ||
| _nu_y = firstStep.Y(); | ||
| // count the number of protons (that start from < 1 cm of neutrino) | ||
| auto const& parts = truth.GetParticles(); | ||
| _nProtons = 0; | ||
| for (auto p : parts){ | ||
| if ( (p.PdgCode() == 2212) and ((p.Trajectory().front().E()-p.Mass()) > 0.02) ){ | ||
| auto const& ps = p.Trajectory().front(); | ||
| double d = sqrt( (ps.X()+_detwidth-_nu_x)*(ps.X()+_detwidth-_nu_x) + | ||
| (ps.Y()-_nu_y)*(ps.Y()-_nu_y) + | ||
| (ps.Z()-_nu_z)*(ps.Z()-_nu_z) ); | ||
| if (d < 1) | ||
| _nProtons += 1; | ||
| }// if proton > 20 MeV | ||
| }// for all parts | ||
|
|
||
|
|
||
| if (_verbose){ std::cout << "Event: " << _evt << std::endl; } | ||
|
|
||
| // prepare 3-element vector to store best hit on each plane | ||
| std::vector<::larlite::hit> bestH_v; | ||
|
|
||
| for (int plane=0; plane < 3; plane++){ | ||
|
|
||
| // get best hit on a plane | ||
| size_t bestHitQ = GetBestHit(evt_hits,plane); | ||
|
|
||
| if (bestHitQ == 9999) | ||
| return false; | ||
|
|
||
| auto const& bestH = evt_hits->at(bestHitQ); | ||
| bestH_v.push_back(bestH); | ||
| double t = bestH.PeakTime()*_time2cm; | ||
| double w = bestH.WireID().Wire*_wire2cm; | ||
|
|
||
| double dist = sqrt( (_nu_z-w)*(_nu_z-w) + (_nu_x-t)*(_nu_x-t) ); | ||
| if (plane == 0) _distU = dist; | ||
| if (plane == 1) _distV = dist; | ||
| if (plane == 2) _distY = dist; | ||
|
|
||
| if (_verbose){ | ||
|
|
||
| std::cout << "Plane: " << plane << std::endl; | ||
| std::cout << "Best Hit : [" << w << ", " << t << "]" << std::endl; | ||
| std::cout << "Neutrino pos: [" << _nu_z << ", " << _nu_x << "]" << std::endl; | ||
| std::cout << "dist: " << dist << std::endl; | ||
| } | ||
|
|
||
| if (plane == 2){ | ||
| // now loop through hit-list and fill vector with angle to start point direction | ||
| _QvsTheta.clear(); | ||
| _QvsTheta = std::vector<double>(360,0.); | ||
| for (size_t h=0; h < evt_hits->size(); h++){ | ||
| auto const& hit = evt_hits->at(h); | ||
| double t2 = hit.PeakTime()*_time2cm; | ||
| double w2 = hit.WireID().Wire*_wire2cm; | ||
| if ( ((w2-w) == 0) and ((t2-t) == 0) ) | ||
| continue; | ||
| double c = (t2-t)/sqrt( (w2-w)*(w2-w) + (t2-t)*(t2-t) ); | ||
| double angle = (180/3.14)*acos(c); | ||
| if ( (w2-w) < 0){ | ||
| if (c < 0) | ||
| angle += 90; | ||
| else | ||
| angle += 270; | ||
| } | ||
| // std::cout << angle << std::endl; | ||
| _QvsTheta[int(angle)] += hit.Integral(); | ||
| } | ||
| }// if plane == 2 | ||
| }// loop over planes | ||
|
|
||
| // Now reconstruct a 3D vertex from the 3 vertices on each plane | ||
| double yUV, zUV; | ||
| double yUY, zUY; | ||
| double yVY, zVY; | ||
| larutil::Geometry::GetME()->IntersectionPoint( bestH_v[0].WireID().Wire, bestH_v[1].WireID().Wire, | ||
| 0, 1, yUV, zUV); | ||
| larutil::Geometry::GetME()->IntersectionPoint( bestH_v[1].WireID().Wire, bestH_v[2].WireID().Wire, | ||
| 1, 2, yVY, zVY); | ||
| larutil::Geometry::GetME()->IntersectionPoint( bestH_v[0].WireID().Wire, bestH_v[2].WireID().Wire, | ||
| 0, 2, yUY, zUY); | ||
|
|
||
| // average the points | ||
| _yAvg = (yUV+yUY+yVY)/3.; | ||
| _zAvg = (zUV+zUY+zVY)/3.; | ||
|
|
||
| // area from the 3 points (measure of uncertainty) | ||
| _area = ( zUV*(yUY-zVY) + zUY*(yVY-yUV) + zVY*(yUV-yUY) ) / 2.; | ||
| if (_area < 0) { _area *= -1; } | ||
|
|
||
| _distYZ = sqrt( (_yAvg-_nu_y)*(_yAvg-_nu_y) + (_zAvg-_nu_z)*(_zAvg-_nu_z) ); | ||
|
|
||
| // best time average | ||
| _xAvg = (bestH_v[0].PeakTime()+bestH_v[1].PeakTime()+bestH_v[2].PeakTime())*_time2cm/3.; | ||
| _distX = _nu_x-_xAvg; | ||
| if (_distX < 0) _distX *= -1; | ||
|
|
||
| _dist3D = sqrt ( _distX*_distX + _distYZ*_distYZ ); | ||
|
|
||
| _vtx_tree->Fill(); | ||
|
|
||
| _evt += 1; | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| bool VertexFinder::finalize() { | ||
|
|
||
| if (_fout && _vtx_tree) | ||
| _vtx_tree->Write(); | ||
|
|
||
| return true; | ||
| } | ||
|
|
||
| void VertexFinder::MakeHitMap(const event_hit* hitlist, int plane){ | ||
|
|
||
| _hitMap.clear(); | ||
| // temporary pair | ||
| std::pair<int,int> tmpPair; | ||
|
|
||
| for (size_t h=0; h < hitlist->size(); h++){ | ||
| auto const& hit = hitlist->at(h); | ||
| // skip if not of plane we want | ||
| if (hit.View() != plane) | ||
| continue; | ||
| double t = hit.PeakTime()*_time2cm; | ||
| double w = hit.WireID().Wire*_wire2cm; | ||
| // map is (i,j) -> hit list | ||
| // i : ith bin in wire of some width | ||
| // j : jth bin in time of some width | ||
| int i = int(w/_cellSize); | ||
| int j = int(t/_cellSize); | ||
| tmpPair = std::make_pair(i,j); | ||
| // does this entry exist in the map? | ||
| // if yes -> append to vector | ||
| // if no create new vector and add to map | ||
| if (_hitMap.find(tmpPair) == _hitMap.end()){ | ||
| std::vector<size_t> aaa = {h}; | ||
| _hitMap[tmpPair] = aaa; | ||
| } | ||
| else | ||
| _hitMap[tmpPair].push_back(h); | ||
| }// for all hits | ||
|
|
||
| return; | ||
| } | ||
|
|
||
| size_t VertexFinder::GetBestHit(const event_hit* evt_hits, int plane){ | ||
|
|
||
|
|
||
| MakeHitMap(evt_hits,plane); | ||
|
|
||
| // index for best hit | ||
| size_t bestHitQ = 9999; | ||
| size_t bestHitN = 9999; | ||
| // charge at best hit | ||
| double bestCharge = -1; | ||
| double bestNum = 0; | ||
|
|
||
| std::map<std::pair<int,int>, std::vector<size_t> >::iterator it; | ||
| for (it = _hitMap.begin(); it != _hitMap.end(); it++){ | ||
| auto const& hitlist = it->second; | ||
| for (size_t h1=0; h1 < hitlist.size(); h1++){ | ||
| auto const& hit1 = evt_hits->at(hitlist[h1]); | ||
| double charge = hit1.Integral(); | ||
| double num = 0; | ||
| for (size_t h2=0; h2 < hitlist.size(); h2++){ | ||
| if (h1 == h2) continue; | ||
| auto const& hit2 = evt_hits->at(hitlist[h2]); | ||
| double t1 = hit1.PeakTime()*_time2cm; | ||
| double w1 = hit1.WireID().Wire*_wire2cm; | ||
| double t2 = hit2.PeakTime()*_time2cm; | ||
| double w2 = hit2.WireID().Wire*_wire2cm; | ||
| double dist = sqrt( (t1-t2)*(t1-t2) + (w1-w2)*(w1-w2) ); | ||
| if (dist < _radius){ | ||
| charge += hit2.Integral(); | ||
| num += 1; | ||
| } | ||
| }// for second loop | ||
| if (charge > bestCharge){ | ||
| bestCharge = charge; | ||
| bestHitQ = hitlist[h1]; | ||
| } | ||
| if (num > bestNum){ | ||
| bestNum = num; | ||
| bestHitN = hitlist[h1]; | ||
| } | ||
| }// for first loop | ||
| }// loop over map entries | ||
|
|
||
| return bestHitQ; | ||
| } | ||
|
|
||
| } | ||
| #endif |
| @@ -0,0 +1,93 @@ | ||
| /** | ||
| * \file VertexFinder.h | ||
| * | ||
| * \ingroup Playground | ||
| * | ||
| * \brief Class def header for a class VertexFinder | ||
| * | ||
| * @author david caratelli | ||
| */ | ||
|
|
||
| /** \addtogroup Playground | ||
| @{*/ | ||
|
|
||
| #ifndef LARLITE_VERTEXFINDER_H | ||
| #define LARLITE_VERTEXFINDER_H | ||
|
|
||
| #include "Analysis/ana_base.h" | ||
| #include <utility> | ||
|
|
||
| namespace larlite { | ||
| /** | ||
| \class VertexFinder | ||
| User custom analysis class made by SHELL_USER_NAME | ||
| */ | ||
| class VertexFinder : public ana_base{ | ||
|
|
||
| public: | ||
|
|
||
| /// Default constructor | ||
| VertexFinder(){ _name="VertexFinder"; _fout=0; _verbose=false; } | ||
|
|
||
| /// Default destructor | ||
| virtual ~VertexFinder(){} | ||
|
|
||
| virtual bool initialize(); | ||
|
|
||
| virtual bool analyze(storage_manager* storage); | ||
|
|
||
| virtual bool finalize(); | ||
|
|
||
| /// Set the size of each cell | ||
| void setCellSize(double d) { _cellSize = d; } | ||
| /// Set the radius around which to search for hits | ||
| void setRadius(double d) { _radius = d; } | ||
| /// Set which plane to select hits from | ||
| void setPlane(int pl) { _plane = pl; } | ||
| /// Verbosity setter | ||
| void setVerbose(bool on) { _verbose = on; } | ||
|
|
||
| protected: | ||
|
|
||
| /// verbosity flag | ||
| bool _verbose; | ||
|
|
||
| /// conversion factors for hits | ||
| double _wire2cm, _time2cm; | ||
|
|
||
| /// map connecting coordinate index (i,j) to [h1,h2,h3] (hit index list) | ||
| std::map<std::pair<int,int>, std::vector<size_t> > _hitMap; | ||
|
|
||
| /// size of each cell [cm] | ||
| double _cellSize; | ||
|
|
||
| /// radius to count charge around [cm] | ||
| double _radius; | ||
|
|
||
| /// plane to select hits from | ||
| int _plane; | ||
|
|
||
| /// Map making function | ||
| void MakeHitMap(const event_hit* hitlist, int plane); | ||
|
|
||
| /// find the best hit on a plane | ||
| size_t GetBestHit(const event_hit* evt_hits, int plane); | ||
|
|
||
| TTree* _vtx_tree; | ||
| double _distU, _distV, _distY; | ||
| double _distYZ, _distX, _dist3D; | ||
| double _nu_x, _nu_y, _nu_z; | ||
| double _xAvg, _yAvg, _zAvg; | ||
| std::vector<double> _QvsTheta; | ||
| int _nProtons; | ||
| int _evt; | ||
| double _area; // area spanned by triangle from projection of 3 points on YV plane | ||
|
|
||
| // detector width | ||
| double _detwidth; | ||
|
|
||
| }; | ||
| } | ||
| #endif | ||
|
|
| @@ -0,0 +1,91 @@ | ||
| import sys,os | ||
|
|
||
| if len(sys.argv) < 2: | ||
| msg = '\n' | ||
| msg += "Usage 1: %s $INPUT_ROOT_FILE\n" % sys.argv[0] | ||
| msg += '\n' | ||
| sys.stderr.write(msg) | ||
| sys.exit(1) | ||
|
|
||
| from ROOT import gSystem | ||
| #gSystem.Load("libSimpleWFAna") | ||
| from ROOT import larlite as fmwk | ||
|
|
||
| # Create ana_processor instance | ||
| my_proc = fmwk.ana_processor() | ||
|
|
||
| # Set input root file | ||
| run=int(sys.argv[1]) | ||
| DATA_DIR="/home/david/uBooNE/data/WireBias/" | ||
| temp_files=[x for x in os.listdir(DATA_DIR) if x.find('LArLite-%03d-' % run)==0] | ||
| files={} | ||
| for f in temp_files: | ||
| subrun=f.replace('LArLite-%03d-' % run,'') | ||
| subrun=int(subrun.replace('.root','')) | ||
| files[subrun]=f | ||
| files = {} | ||
| files[0] = 'LArTF_Test_519_larlight.root' | ||
| subrun_start = files.keys()[0] | ||
| subrun_end = files.keys()[-1] | ||
| if len(sys.argv)==4: | ||
| subrun_start = int(sys.argv[2]) | ||
| subrun_end = int(sys.argv[3]) | ||
| my_proc.set_ana_output_file("corr_run%03d_subrun%03d_%03d.root" % (run,subrun_start,subrun_end)) | ||
| else: | ||
| my_proc.set_ana_output_file("corr_run%03d.root" % run); | ||
|
|
||
| for x in files.keys(): | ||
| my_proc.add_input_file('%s/%s' % (DATA_DIR,files[x])) | ||
|
|
||
| #my_proc.set_input_rootdir("scanner") | ||
| # Specify IO mode | ||
| my_proc.set_io_mode(fmwk.storage_manager.kREAD) | ||
|
|
||
| # Specify output root file name | ||
| my_proc.set_ana_output_file("noise_zigzag_run%03d.root" % run); | ||
|
|
||
| corr_unit = fmwk.ZigZag() | ||
| #corr_unit.setVerbose(True) | ||
| corr_unit.setZigZagMinLength(4) | ||
|
|
||
| contents = open('mapping_v01.txt','r') | ||
| crate_slot_map={} | ||
| for line in contents: | ||
|
|
||
| if len(line.split()) < 4: continue | ||
| (larch,crate,slot,ch) = [int(x) for x in line.split()] | ||
| corr_unit.SetMap(larch,crate,slot,ch) | ||
| if not crate in crate_slot_map.keys(): | ||
| crate_slot_map[crate]={} | ||
| if not slot in crate_slot_map[crate].keys(): | ||
| crate_slot_map[crate][slot]={} | ||
| if not ch in crate_slot_map[crate][slot].keys(): | ||
| crate_slot_map[crate][slot][ch]=larch | ||
|
|
||
|
|
||
| #reference = open('%s/noise_ch.txt' % os.environ['HOME'],'r').read().split('\n') | ||
| reference = open("zigzag_small.txt",'r') | ||
| for line in reference: | ||
| chans = line.split(",") | ||
| for l in chans: | ||
| print l | ||
| corr_unit.add(int(l)) | ||
|
|
||
| my_proc.add_process(corr_unit) | ||
| #my_proc.add_process(fmwk.SimpleWFAna()); | ||
|
|
||
| #my_proc.set_verbosity(fmwk.MSG.DEBUG) | ||
| print "Finished configuring ana_processor. Start event loop!" | ||
|
|
||
| # Let's run it. | ||
| my_proc.run(0,10) | ||
|
|
||
| # done! | ||
| print "Finished running ana_processor event loop!" | ||
|
|
||
| sys.exit(0) | ||
|
|
| @@ -0,0 +1,40 @@ | ||
| # | ||
| # Example PyROOT script to run analysis module, ana_base. | ||
| # The usage is same for inherited analysis class instance. | ||
| # | ||
|
|
||
| # Load libraries | ||
| import os, ROOT, sys | ||
| from ROOT import * | ||
| from ROOT import gSystem | ||
| import time | ||
|
|
||
| gSystem.Load("libLArLite_Base") | ||
| gSystem.Load("libLArLite_Analysis") | ||
| gSystem.Load("libLArLite_LArUtil") | ||
| gSystem.Load("libBasicTool_GeoAlgo") | ||
|
|
||
| # Now import ana_processor & your class. For this example, ana_base. | ||
| from ROOT import larlite as fmwk | ||
|
|
||
| # Create ana_processor instance | ||
| my_proc=fmwk.ana_processor() | ||
|
|
||
| # Specify IO mode | ||
| my_proc.set_io_mode(fmwk.storage_manager.kREAD) | ||
| #my_proc.set_io_mode(storage_manager.WRITE) | ||
| #my_proc.set_io_mode(storage_manager.BOTH) | ||
|
|
||
| for x in xrange(len(sys.argv)-1): | ||
| my_proc.add_input_file(sys.argv[x+1]) | ||
|
|
||
| my_proc.set_ana_output_file("hitenergy.root") | ||
|
|
||
|
|
||
| ana = fmwk.HitEnergy() | ||
| ana.useShowers(True) | ||
|
|
||
| my_proc.add_process(ana) | ||
|
|
||
| my_proc.run() | ||
| sys.exit(0); |
| @@ -0,0 +1,137 @@ | ||
| # | ||
| # Example PyROOT script to run analysis module, ana_base. | ||
| # The usage is same for inherited analysis class instance. | ||
| # | ||
|
|
||
| # Load libraries | ||
| import os, ROOT, sys | ||
| from ROOT import * | ||
| from ROOT import gSystem | ||
| import time | ||
| import matplotlib.pyplot as plt | ||
| import numpy as np | ||
|
|
||
| plt.ion() | ||
|
|
||
| # function that given hit charge returns marker size | ||
| def markerSize(Q): | ||
| return (Q/10.) | ||
|
|
||
| gSystem.Load("libLArLite_Base") | ||
| gSystem.Load("libLArLite_Analysis") | ||
| gSystem.Load("libLArLite_LArUtil") | ||
| gSystem.Load("libBasicTool_GeoAlgo") | ||
|
|
||
| # Now import ana_processor & your class. For this example, ana_base. | ||
| from ROOT import larlite as fmwk | ||
|
|
||
| w2cm = larutil.GeometryUtilities.GetME().WireToCm() | ||
| t2cm = larutil.GeometryUtilities.GetME().TimeToCm() | ||
|
|
||
| # Create ana_processor instance | ||
| my_proc=fmwk.ana_processor() | ||
|
|
||
| # Specify IO mode | ||
| my_proc.set_io_mode(fmwk.storage_manager.kBOTH) | ||
|
|
||
| for x in xrange(len(sys.argv)-1): | ||
| my_proc.add_input_file(sys.argv[x+1]) | ||
|
|
||
| my_proc.set_ana_output_file("hitremover.root") | ||
| my_proc.set_output_file("out_hits.root") | ||
|
|
||
| # Create analysis class instance. For this example, ana_base. | ||
| # To show how one can run multiple analysis modules at once, | ||
| # we make multiple ana_base instance. | ||
|
|
||
| my_ana = fmwk.TrackHitRemover() | ||
| outhitproducer = 'outhits' | ||
| my_ana.setOutHitProducer(outhitproducer) | ||
| #my_ana.setTrackProducer("cctrack") | ||
| #my_ana.setTrackProducer("pandoraCosmicKHit") | ||
| #my_ana.setTrackProducer("pandoraNuKHit") | ||
| my_ana.setTrackProducer("trackkalmanhit") | ||
| #my_ana.setTrackProducer("stitchkalmanhit") | ||
| my_proc.add_process(my_ana) | ||
| my_proc.set_data_to_write(fmwk.data.kHit,outhitproducer) | ||
| #my_proc.run() | ||
| #sys.exit() | ||
|
|
||
| # track colors | ||
| trk_c = ['r','g','c','m','b'] | ||
|
|
||
| fig, (axU, axV, axY) = plt.subplots(nrows=3,figsize=(15,12)) | ||
|
|
||
| while my_proc.process_event(): | ||
|
|
||
| hits = my_ana.getTrkHits() | ||
| out_hit_indices = my_ana.getOutHits() | ||
| print "there are ", hits.size(), " hits" | ||
| hit_t = [[],[],[]] | ||
| hit_w = [[],[],[]] | ||
| hit_s = [[],[],[]] | ||
| axU.clear() | ||
| axV.clear() | ||
| axY.clear() | ||
| fig.gca() | ||
| # loop over all hits and color them blue | ||
| for h in xrange(hits.size()): | ||
| #for h in out_hit_indices: | ||
| pl = hits.at(h).View() | ||
| hit_w[pl].append(hits.at(h).WireID().Wire*w2cm) | ||
| hit_t[pl].append(hits.at(h).PeakTime()*t2cm) | ||
| hit_s[pl].append(markerSize(hits.at(h).Integral())) | ||
| axU.scatter(hit_w[0],hit_t[0],s=hit_s[0],color='b') | ||
| axV.scatter(hit_w[1],hit_t[1],s=hit_s[1],color='b') | ||
| axY.scatter(hit_w[2],hit_t[2],s=hit_s[2],color='b') | ||
|
|
||
| ''' | ||
| # now loop over hits with association to track and color them red | ||
| ass = my_ana.getAss() | ||
| ctr = 0 | ||
| for t in xrange(ass.size()): | ||
| # remove tracks with less than X hits | ||
| if (ass.at(t).size() < 0): | ||
| continue | ||
| # prepare vector for new cluster | ||
| hit_t_ass = [] | ||
| hit_w_ass = [] | ||
| hit_s_ass = [] | ||
| for h2 in xrange(ass.at(t).size()): | ||
| # don't consider Y plane hits | ||
| if (hits.at(ass.at(t).at(h2)).View() != 2): | ||
| continue | ||
| hit_w_ass.append(hits.at(ass.at(t).at(h2)).WireID().Wire*w2cm) | ||
| hit_t_ass.append(hits.at(ass.at(t).at(h2)).PeakTime()*t2cm) | ||
| hit_s_ass.append(markerSize(hits.at(ass.at(t).at(h2)).Integral())) | ||
| plt.scatter(hit_w_ass,hit_t_ass,s=hit_s_ass,color=trk_c[ctr%5]) | ||
| ctr += 1 | ||
| ''' | ||
|
|
||
| #''' | ||
| # now get hits associated with track-like clusters and add them | ||
| trk_hits = my_ana.getTrkHitIndices() | ||
|
|
||
| ctr = 0 | ||
| for t in xrange(trk_hits.size()): | ||
| trk_hit_t = [[],[],[]] | ||
| trk_hit_w = [[],[],[]] | ||
| trk_hit_s = [[],[],[]] | ||
| for th in xrange(trk_hits[t].size()): | ||
| # only Y-plane hits for now | ||
| pl = hits.at(trk_hits[t][th]).View() | ||
| trk_hit_w[pl].append(hits.at(trk_hits[t][th]).WireID().Wire*w2cm) | ||
| trk_hit_t[pl].append(hits.at(trk_hits[t][th]).PeakTime()*t2cm) | ||
| trk_hit_s[pl].append(markerSize(hits.at(trk_hits[t][th]).Integral())) | ||
|
|
||
| axU.scatter(trk_hit_w[0],trk_hit_t[0],s=trk_hit_s[0],color=trk_c[ctr%5]) | ||
| axV.scatter(trk_hit_w[1],trk_hit_t[1],s=trk_hit_s[1],color=trk_c[ctr%5]) | ||
| axY.scatter(trk_hit_w[2],trk_hit_t[2],s=trk_hit_s[2],color=trk_c[ctr%5]) | ||
| ctr += 1 | ||
| #''' | ||
|
|
||
| plt.show() | ||
| usrinput = raw_input("Hit Enter: next evt || q: exit viewer\n") | ||
| if ( usrinput == "q" ): | ||
| sys.exit(0) | ||
|
|
| @@ -0,0 +1,45 @@ | ||
| import sys | ||
|
|
||
| if len(sys.argv) < 2: | ||
| msg = '\n' | ||
| msg += "Usage 1: %s $INPUT_ROOT_FILE\n" % sys.argv[0] | ||
| msg += '\n' | ||
| sys.stderr.write(msg) | ||
| sys.exit(1) | ||
|
|
||
| from ROOT import larlite as fmwk | ||
|
|
||
|
|
||
| # Create ana_processor instance | ||
| my_proc = fmwk.ana_processor() | ||
|
|
||
| # Set input root files | ||
| for x in xrange(len(sys.argv)): | ||
| if x < 1: | ||
| continue | ||
| my_proc.add_input_file(sys.argv[x]) | ||
|
|
||
| # Specify IO mode | ||
| my_proc.set_io_mode(fmwk.storage_manager.kREAD) | ||
|
|
||
| # Specify analysis output root file name | ||
| my_proc.set_ana_output_file("simple_shr_qual.root"); | ||
|
|
||
| # Specify data output root file name | ||
| my_proc.set_output_file("outfile.root") | ||
|
|
||
| # Create analysis unit | ||
| module = fmwk.ShowerQualitySimple() | ||
| module.setProducerName("showerreco") | ||
| module.setMinE(100) | ||
|
|
||
| my_proc.add_process(module) | ||
|
|
||
| print "Finished configuring ana_processor. Start event loop!" | ||
|
|
||
| my_proc.run() | ||
|
|
||
| sys.exit(0) | ||
|
|
| @@ -36,5 +36,5 @@ | ||
|
|
||
| my_proc.add_process(ana) | ||
|
|
||
| my_proc.run() | ||
| sys.exit(0); | ||
| @@ -0,0 +1,108 @@ | ||
| # | ||
| # Example PyROOT script to run analysis module, ana_base. | ||
| # The usage is same for inherited analysis class instance. | ||
| # | ||
|
|
||
| # Load libraries | ||
| import os, ROOT, sys | ||
| from ROOT import * | ||
| from ROOT import gSystem | ||
| import time | ||
| import matplotlib.pyplot as plt | ||
| import numpy as np | ||
|
|
||
| plt.ion() | ||
|
|
||
| # function that given hit charge returns marker size | ||
| def markerSize(Q): | ||
| return (Q/10.) | ||
|
|
||
| gSystem.Load("libLArLite_Base") | ||
| gSystem.Load("libLArLite_Analysis") | ||
| gSystem.Load("libLArLite_LArUtil") | ||
| gSystem.Load("libBasicTool_GeoAlgo") | ||
|
|
||
| # Now import ana_processor & your class. For this example, ana_base. | ||
| from ROOT import larlite as fmwk | ||
|
|
||
| w2cm = larutil.GeometryUtilities.GetME().WireToCm() | ||
| t2cm = larutil.GeometryUtilities.GetME().TimeToCm() | ||
|
|
||
| # Create ana_processor instance | ||
| my_proc=fmwk.ana_processor() | ||
|
|
||
| # Specify IO mode | ||
| my_proc.set_io_mode(fmwk.storage_manager.kBOTH) | ||
|
|
||
| for x in xrange(len(sys.argv)-1): | ||
| my_proc.add_input_file(sys.argv[x+1]) | ||
|
|
||
| my_proc.set_ana_output_file("hitremover.root") | ||
| my_proc.set_output_file("out_hits.root") | ||
|
|
||
| # Create analysis class instance. For this example, ana_base. | ||
| # To show how one can run multiple analysis modules at once, | ||
| # we make multiple ana_base instance. | ||
|
|
||
| my_ana = fmwk.TrackHitRemover() | ||
| outhitproducer = 'outhits' | ||
| my_ana.setOutHitProducer(outhitproducer) | ||
| my_ana.setGetShowerHits(True) | ||
| my_ana.setShowerProducer("showerreco") | ||
| my_proc.add_process(my_ana) | ||
| my_proc.set_data_to_write(fmwk.data.kHit,outhitproducer) | ||
| #my_proc.run() | ||
| #sys.exit() | ||
|
|
||
| # track colors | ||
| trk_c = ['r','g','c','m','b'] | ||
|
|
||
| fig, (axU, axV, axY) = plt.subplots(nrows=3,figsize=(15,12)) | ||
|
|
||
| while my_proc.process_event(): | ||
|
|
||
| hits = my_ana.getShrHits() | ||
| out_hit_indices = my_ana.getOutHits() | ||
| hit_t = [[],[],[]] | ||
| hit_w = [[],[],[]] | ||
| hit_s = [[],[],[]] | ||
| axU.clear() | ||
| axV.clear() | ||
| axY.clear() | ||
| fig.gca() | ||
| # loop over all hits and color them blue | ||
| for h in xrange(hits.size()): | ||
| #for h in out_hit_indices: | ||
| pl = hits.at(h).View() | ||
| hit_w[pl].append(hits.at(h).WireID().Wire*w2cm) | ||
| hit_t[pl].append(hits.at(h).PeakTime()*t2cm) | ||
| hit_s[pl].append(markerSize(hits.at(h).Integral())) | ||
| axU.scatter(hit_w[0],hit_t[0],s=hit_s[0],color='b') | ||
| axV.scatter(hit_w[1],hit_t[1],s=hit_s[1],color='b') | ||
| axY.scatter(hit_w[2],hit_t[2],s=hit_s[2],color='b') | ||
|
|
||
| # now get hits associated with track-like clusters and add them | ||
| shr_hits = my_ana.getShrHitIndices() | ||
|
|
||
| ctr = 0 | ||
| for t in xrange(shr_hits.size()): | ||
| trk_hit_t = [[],[],[]] | ||
| trk_hit_w = [[],[],[]] | ||
| trk_hit_s = [[],[],[]] | ||
| for th in xrange(shr_hits[t].size()): | ||
| # determine plane | ||
| pl = hits.at(shr_hits[t][th]).View() | ||
| trk_hit_w[pl].append(hits.at(shr_hits[t][th]).WireID().Wire*w2cm) | ||
| trk_hit_t[pl].append(hits.at(shr_hits[t][th]).PeakTime()*t2cm) | ||
| trk_hit_s[pl].append(markerSize(hits.at(shr_hits[t][th]).Integral())) | ||
|
|
||
| axU.scatter(trk_hit_w[0],trk_hit_t[0],s=trk_hit_s[0],color=trk_c[ctr%5]) | ||
| axV.scatter(trk_hit_w[1],trk_hit_t[1],s=trk_hit_s[1],color=trk_c[ctr%5]) | ||
| axY.scatter(trk_hit_w[2],trk_hit_t[2],s=trk_hit_s[2],color=trk_c[ctr%5]) | ||
| ctr += 1 | ||
|
|
||
| plt.show() | ||
| usrinput = raw_input("Hit Enter: next evt || q: exit viewer\n") | ||
| if ( usrinput == "q" ): | ||
| sys.exit(0) | ||
|
|
| @@ -0,0 +1,90 @@ | ||
| # Load libraries | ||
| import sys, os | ||
| from ROOT import * | ||
| from ROOT import gSystem | ||
| from ROOT import larlite as fmwk | ||
|
|
||
| # vertex finder | ||
| # Create ana_processor instance | ||
| vtxproc=fmwk.ana_processor() | ||
|
|
||
| # Specify IO mode | ||
| vtxproc.set_io_mode(fmwk.storage_manager.kREAD) | ||
|
|
||
| for x in xrange(len(sys.argv)-1): | ||
| vtxproc.add_input_file(sys.argv[x+1]) | ||
|
|
||
| vtxproc.set_ana_output_file("vtxfinder.root") | ||
|
|
||
| # Create analysis class instance. For this example, ana_base. | ||
| # To show how one can run multiple analysis modules at once, | ||
| # we make multiple ana_base instance. | ||
|
|
||
| vtxfinder = fmwk.VertexFinder() | ||
| vtxfinder.setCellSize(2.) | ||
| vtxfinder.setRadius(1.) | ||
| vtxfinder.setPlane(2) | ||
|
|
||
| vtxproc.add_process(vtxfinder) | ||
|
|
||
| vtxproc.run() | ||
|
|
||
| # Create ana_processor instance | ||
| my_proc=fmwk.ana_processor() | ||
|
|
||
| # Specify IO mode | ||
| #my_proc.set_io_mode(fmwk.storage_manager.kBOTH) | ||
| #my_proc.set_io_mode(fmwk.storage_manager.kWRITE) | ||
| my_proc.set_io_mode(fmwk.storage_manager.kREAD) | ||
|
|
||
|
|
||
| for x in xrange(len(sys.argv)-1): | ||
| my_proc.add_input_file(sys.argv[x+1]) | ||
|
|
||
|
|
||
| #is there a way to disable ana_proc from creating an output file at all? | ||
| my_proc.set_ana_output_file("") | ||
|
|
||
| my_ana = fmwk.HitViewer() | ||
| my_ana.setMinQ(200) | ||
|
|
||
| my_proc.add_process(my_ana) | ||
|
|
||
|
|
||
| gStyle.SetOptStat(0) | ||
| gStyle.SetTitleFontSize(0.1) | ||
| gStyle.SetTitleOffset(0.4,"X") | ||
| gStyle.SetTitleSize(0.08,"X") | ||
| gStyle.SetTitleOffset(0.4,"Y") | ||
| gStyle.SetTitleSize(0.08,"Y") | ||
| gStyle.SetLabelSize(0.08,"X") | ||
| gStyle.SetLabelSize(0.08,"Y") | ||
| gStyle.SetLabelSize(0.08,"Z") | ||
| gStyle.SetOptLogz(1) | ||
|
|
||
| c=TCanvas("c","Wire v. Time Hit Viewer",900,600) | ||
| c.Divide(1,3) | ||
|
|
||
| while my_proc.process_event(): | ||
|
|
||
| #vtxproc.process_event() | ||
|
|
||
| currentview = 0; | ||
| #First fill the 6 pads on the main canvas with stuff | ||
| for pad in xrange(1,4,1): | ||
|
|
||
| c.cd(pad) | ||
| h = my_ana.GetHisto_Hits(int(currentview)) | ||
| h.SetTitleFont(62,'X') | ||
| h.SetTitleFont(62,'Y') | ||
| h.Draw("COLZ") | ||
|
|
||
| currentview = currentview + 1 | ||
| c.Update() | ||
|
|
||
|
|
||
| usrinput = raw_input("Hit Enter: next evt || q: exit viewer\n") | ||
| if ( usrinput == "q" ): | ||
| sys.exit(0) | ||
|
|
||
| # done! |
| @@ -0,0 +1,82 @@ | ||
| # | ||
| # Example PyROOT script to run analysis module, ana_base. | ||
| # The usage is same for inherited analysis class instance. | ||
| # | ||
|
|
||
| # Load libraries | ||
| import os, ROOT, sys | ||
| from ROOT import * | ||
| from ROOT import gSystem | ||
| import time | ||
|
|
||
| gSystem.Load("libLArLite_Base") | ||
| gSystem.Load("libLArLite_Analysis") | ||
| gSystem.Load("libLArLite_LArUtil") | ||
| gSystem.Load("libBasicTool_GeoAlgo") | ||
|
|
||
| # Now import ana_processor & your class. For this example, ana_base. | ||
| from ROOT import larlite as fmwk | ||
|
|
||
| # Create ana_processor instance | ||
| my_proc=fmwk.ana_processor() | ||
|
|
||
| # Specify IO mode | ||
| my_proc.set_io_mode(fmwk.storage_manager.kREAD) | ||
|
|
||
| for x in xrange(len(sys.argv)-1): | ||
| my_proc.add_input_file(sys.argv[x+1]) | ||
|
|
||
| my_proc.set_ana_output_file("vtxfinder.root") | ||
|
|
||
| # Create analysis class instance. For this example, ana_base. | ||
| # To show how one can run multiple analysis modules at once, | ||
| # we make multiple ana_base instance. | ||
|
|
||
| vtxfinder = fmwk.VertexFinder() | ||
| vtxfinder.setCellSize(4.) | ||
| vtxfinder.setRadius(2.) | ||
| vtxfinder.setPlane(2) | ||
| vtxfinder.setVerbose(False) | ||
| my_ana = fmwk.HitViewer() | ||
| my_proc.add_process(vtxfinder) | ||
| #my_proc.add_process(my_ana) | ||
| my_proc.run() | ||
| sys.exit() | ||
| ''' | ||
| gStyle.SetOptStat(0) | ||
| gStyle.SetTitleFontSize(0.1) | ||
| gStyle.SetTitleOffset(0.4,"X") | ||
| gStyle.SetTitleSize(0.08,"X") | ||
| gStyle.SetTitleOffset(0.4,"Y") | ||
| gStyle.SetTitleSize(0.08,"Y") | ||
| gStyle.SetLabelSize(0.08,"X") | ||
| gStyle.SetLabelSize(0.08,"Y") | ||
| gStyle.SetLabelSize(0.08,"Z") | ||
| gStyle.SetOptLogz(1) | ||
| c=TCanvas("c","Wire v. Time Hit Viewer",900,600) | ||
| c.Divide(1,3) | ||
| while my_proc.process_event(): | ||
| #vtxproc.process_event() | ||
| currentview = 0; | ||
| #First fill the 6 pads on the main canvas with stuff | ||
| for pad in xrange(1,4,1): | ||
| c.cd(pad) | ||
| h = my_ana.GetHisto_Hits(int(currentview)) | ||
| h.SetTitleFont(62,'X') | ||
| h.SetTitleFont(62,'Y') | ||
| h.Draw("COLZ") | ||
| currentview = currentview + 1 | ||
| c.Update() | ||
| usrinput = raw_input("Hit Enter: next evt || q: exit viewer\n") | ||
| if ( usrinput == "q" ): | ||
| sys.exit(0) | ||
| ''' |