diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 2cf3326..570e572 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -15,16 +15,15 @@ defaults: shell: bash jobs: + build_and_test: + name: Build and test runs-on: ubuntu-latest container: image: codecr.jlab.org/hallb/clas12/container-forge/base_root:latest - steps: - - name: checkout repository uses: actions/checkout@v5 - - name: set clas12root env vars run: | CLAS12ROOT=$GITHUB_WORKSPACE @@ -35,14 +34,49 @@ jobs: echo PATH=$PATH | tee -a $GITHUB_ENV echo LD_LIBRARY_PATH=$LD_LIBRARY_PATH | tee -a $GITHUB_ENV echo ROOT_INCLUDE_PATH=$ROOT_INCLUDE_PATH | tee -a $GITHUB_ENV - - name: print versions run: | cmake --version root --version - - name: build clas12root run: ./installC12Root - - name: run example on data file run: ./bin/clas12root -q ./tests/read_file.C + + generate_documentation: + name: Generate documentation + runs-on: ubuntu-latest + container: + image: codecr.jlab.org/hallb/clas12/container-forge/base_root:latest + options: --user root + steps: + - uses: actions/checkout@v5 + - name: install dependencies + run: | + pacman -Sy --noconfirm + pacman -S --noconfirm doxygen graphviz + - name: doxygen + run: docs/doxygen/generate.sh + - name: rename + run: mv share/doc/clas12root/html pages + - uses: actions/upload-pages-artifact@v4 + with: + retention-days: 3 + path: pages/ + + deploy_webpages: + if: ${{ github.ref == 'refs/heads/master' }} + name: Deploy documentation webpage + needs: + - generate_documentation + permissions: + pages: write + id-token: write + environment: + name: github-pages + url: ${{ steps.deployment.outputs.page_url }} + runs-on: ubuntu-latest + steps: + - name: deployment + id: deployment + uses: actions/deploy-pages@v4 diff --git a/.gitignore b/.gitignore index ba54cd5..37af0fc 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ bin build lib man +/share # dependency installations ############## @@ -54,6 +55,7 @@ ehthumbs.db Thumbs.db .nfs* *~ +.cache #Dropbox .dropbox.attr diff --git a/.ignore b/.ignore new file mode 100644 index 0000000..f1b711c --- /dev/null +++ b/.ignore @@ -0,0 +1,8 @@ +# ignore configuration file for some IDE tools + +# submodule source code +lz4/** +clas12-qadb/** +ccdb/** +hipo_src/** +rcdb/** diff --git a/CMakeLists.txt b/CMakeLists.txt index 8f0e2c0..78cb09f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,7 +13,7 @@ endif() message(STATUS "Using cmake version ${CMAKE_VERSION}" ) # Define the project name, version, and language -project(clas12root VERSION 1.9 LANGUAGES CXX) +project(clas12root VERSION 1.9.0 LANGUAGES CXX) # ============================================================================== # Options diff --git a/Clas12Banks/clas12reader.cpp b/Clas12Banks/clas12reader.cpp index 43502cb..222b69c 100644 --- a/Clas12Banks/clas12reader.cpp +++ b/Clas12Banks/clas12reader.cpp @@ -56,6 +56,7 @@ namespace clas12 { addBank(additionalBank); _applyQA=other._applyQA; + _readEventUserAction=other._readEventUserAction; } void clas12reader::initReader(){ @@ -461,8 +462,9 @@ namespace clas12 { _event.getStructure(*ibank.get()); } - - return true; + // now that we have read all the banks, call user's custom read action; + // return its return value, since it's the last thing `readEvent()` does + return _readEventUserAction(this); } //////////////////////////////////////////////////////// ///initialise next event from the reader @@ -643,19 +645,35 @@ namespace clas12 { return true; } - //////////////////////////////////////////////////////// - ///Filter and return detParticles by given PID - std::vector clas12reader::getByID(int id){ - return container_filter(_detParticles, [id](region_part_ptr dr) - {return dr->getPid()==id;}); + //////////////////////////////////////////////////////// + // Filter and return detParticles + + std::vector clas12reader::getDetParticles(bool const& applyBankFilter) const + { + return applyBankFilter ? + container_filter(_detParticles, [](region_part_ptr dr) { return dr->isAllowed(); }) : + _detParticles; + } + + std::vector clas12reader::getByID(int id, bool const& applyBankFilter) const + { + return applyBankFilter ? + container_filter(_detParticles, [&id](region_part_ptr dr) { return dr->getPid()==id && dr->isAllowed(); }) : + container_filter(_detParticles, [&id](region_part_ptr dr) { return dr->getPid()==id; }); } - std::vector clas12reader::getByRegion(int ir){ - return container_filter(_detParticles, [ir](region_part_ptr dr) - {return dr->getRegion()==ir;}); + + std::vector clas12reader::getByRegion(int ir, bool const& applyBankFilter) const + { + return applyBankFilter ? + container_filter(_detParticles, [&ir](region_part_ptr dr) { return dr->getRegion()==ir && dr->isAllowed(); }) : + container_filter(_detParticles, [&ir](region_part_ptr dr) { return dr->getRegion()==ir; }); } - std::vector clas12reader::getByCharge(int ch){ - return container_filter(_detParticles, [ch](region_part_ptr dr) - {return dr->par()->getCharge()==ch;}); + + std::vector clas12reader::getByCharge(int ch, bool const& applyBankFilter) const + { + return applyBankFilter ? + container_filter(_detParticles, [&ch](region_part_ptr dr) { return dr->par()->getCharge()==ch && dr->isAllowed(); }) : + container_filter(_detParticles, [&ch](region_part_ptr dr) { return dr->par()->getCharge()==ch; }); } //////////////////////////////////////////////////////////////// diff --git a/Clas12Banks/clas12reader.h b/Clas12Banks/clas12reader.h index bb83fba..676cd02 100644 --- a/Clas12Banks/clas12reader.h +++ b/Clas12Banks/clas12reader.h @@ -44,6 +44,7 @@ #include #include #include +#include #ifdef RCDB_MYSQL #include "rcdb_reader.h" @@ -58,6 +59,18 @@ namespace clas12 { using std::endl; using std::cerr; + /// @brief Clas12root HIPO file reader + /// @details + /// @par Accessors + /// - Use the following methods to access `region_particle` objects: + /// - `getDetParticles` + /// - `getDetParticlesPtr` + /// - @link getDetParticles(bool const&) const `getDetParticles(bool)` @endlink + /// - `getByID` + /// - `getByRegion` + /// - `getByCharge` + /// - To access `hipo::bank` objects, @ref bank_acc "see the list of bank accessors". + /// @see `clas12root::HipoChain` to read a set of HIPO files class clas12reader { @@ -114,19 +127,89 @@ namespace clas12 { _rbands.push_back(std::move(reg)); } - + // bank accessors + // - the method name convention is `getBankName` where `BankName` is the bank name without the colons (`::`) + // - the `class` tag is needed disambiguate between classes and `clas12reader` methods with the same name + /// @bank_accessor{REC::Particle} + /// @see getParticleBank + class particle& getRECParticle() const { return *_bparts; } + /// @bank_accessor{REC::FTParticle} + /// @see getParticleBank + class ftbparticle& getRECFTParticle() const { return *_bftbparts; } + /// @bank_accessor{HEL::online} + class helonline& getHELonline() const { return *_bhelonline; } + /// @bank_accessor{HEL::flip} + class helflip& getHELflip() const { return *_bhelflip; } + /// @bank_accessor{RUN::config} + class runconfig& getRUNconfig() const { return *_brunconfig; } + /// @bank_accessor{REC::Event} + class event& getRECEvent() const { return *_bevent; } + /// @bank_accessor{REC::FTEvent} + class ftbevent& getRECFTEvent() const { return *_bftbevent; } + /// @bank_accessor{RAW::vtp} + class vtp& getRAWvtp() const { return *_bvtp; } + /// @bank_accessor{REC::VertDoca} + class vertdoca& getRECVertDoca() const { return *_bvertdoca; } + /// @bank_accessor{REC::Calorimeter} + class calorimeter& getRECCalorimeter() const { return *_bcal; } + /// @bank_accessor{REC::Scintillator} + class scintillator& getRECScintillator() const { return *_bscint; } + /// @bank_accessor{REC::Track} + class tracker& getRECTrack() const { return *_btrck; } + /// @bank_accessor{REC::CovMat} + class covmatrix& getRECCovMat() const { return *_bcovmat; } + /// @bank_accessor{REC::UTrack} + class utracker& getRECUTrack() const { return *_butrck; } + /// @bank_accessor{REC::Traj} + class traj& getRECTraj() const { return *_btraj; } + /// @bank_accessor{REC::Cherenkov} + class cherenkov& getRECCherenkov() const { return *_bcher; } + /// @bank_accessor{RICH::Particle} + class rich& getRICHParticle() const { return *_brich; } + /// @bank_accessor{REC::ForwardTagger} + class forwardtagger& getRECForwardTagger() const { return *_bft; } + /// @bank_accessor{MC::Lund} + class mcparticle& getMCLund() const { return *_bmcparts; } + /// @bank_accessor{MC::Event} + class mcevent& getMCEvent() const { return *_bmcevent; } + + // bank pointer accessors + /// @bank_ptr_accessor{HEL::online|getHELonline} helonline_ptr helonline() const{return _bhelonline.get();}; + /// @bank_ptr_accessor{HEL::flip|getHELflip} helflip_ptr helflip() const{return _bhelflip.get();}; + /// @bank_ptr_accessor{RUN::config|getRUNconfig} runconfig_ptr runconfig() const{return _brunconfig.get();}; + /// @bank_ptr_accessor{REC::Event|getRECEvent} event_ptr event() const{return _bevent.get();}; + /// @bank_ptr_accessor{REC::FTEvent|getRECFTEvent} ftbevent_ptr ftbevent() const{return _bftbevent.get();}; + /// @bank_ptr_accessor{RAW::vtp|getRAWvtp} vtp_ptr vtp() const{return _bvtp.get();}; + /// @bank_ptr_accessor{REC::VertDoca|getRECVertDoca} vertdoca_ptr vertdoca() const{return _bvertdoca.get();}; - - + /// @bank_ptr_accessor{MC::Lund|getMCLund} mcpar_ptr mcparts() const{return _bmcparts.get();}; + /// @bank_ptr_accessor{MC::Event|getMCEvent} mcevt_ptr mcevent() const{return _bmcevent.get();}; + // additional bank accessors + + /// @brief get the particle bank, depending on your configuration + /// @returns a reference to `REC::Particle` by default, or a reference to `RECFT::Particle` if you are using + /// FT-based PID (you have called `useFTBased`) and the `RECFT::Particle` bank is not empty + /// (the empty check ignores `hipo::bank::rowlist` filtering, _e.g._, applied by an Iguana filter algorithm) + /// @note this method returns a base-class reference, `hipo::bank&`, while many other bank reference accessors + /// return derived-class references, such as `getRECParticle` + /// @see specific particle-bank accessors: + /// - getRECParticle + /// - getRECFTParticle + hipo::bank& getParticleBank() const { + if(_useFTBased && getRECFTParticle().getRows()>0) + return getRECFTParticle(); + else + return getRECParticle(); + } //support for generic non-DST banks uint addBank(const std::string& name){ @@ -170,12 +253,66 @@ namespace clas12 { } ///////////////////////////////// - + // region particle accessors + + /// @brief Get the number of `region_particle`s + /// @note This method does not take into account any `hipo::bank` filters (_e.g._, from iguana) + /// @returns The number of `region_particle`s + int getNParticles() const noexcept{return _detParticles.size();} + + /// @brief Get the _full_ list of `region_particle`s + /// @note This method does not take into account any `hipo::bank` filters (_e.g._, from iguana) + /// @returns A reference to the _full_ list of `region_particle`s + /// @see Use @link getDetParticles(bool const&) const `getDetParticles(bool)` @endlink if you need a `hipo::bank` filter applied std::vector& getDetParticles(){return _detParticles;} + + /// @brief Get the _full_ list of `region_particle`s + /// @note This method does not take into account any `hipo::bank` filters (_e.g._, from iguana) + /// @returns A pointer to the _full_ list of `region_particle`s + /// @see Use @link getDetParticles(bool const&) const `getDetParticles(bool)` @endlink if you need a `hipo::bank` filter applied std::vector* getDetParticlesPtr(){return &_detParticles;} - std::vector getByID(int id); - std::vector getByRegion(int ir); - std::vector getByCharge(int ch); + + /// @brief Get the list of `region_particle`s + /// @param applyBankFilter If `true`, apply any `hipo::bank` filters (_e.g._, from iguana) + /// @returns A _copy_ of the list of `region_particle`s, with or without applying a `hipo::bank` filter + /// @see Use getDetParticles() or getDetParticlesPtr() if you don't need the filter, and you want to avoid copying + /// @see Additional methods: + /// - getByID + /// - getByRegion + /// - getByCharge + std::vector getDetParticles(bool const& applyBankFilter) const; + + /// @brief Get the list of `region_particle`s, filtered by PDG + /// @returns A _copy_ of the list of corresponding `region_particle`s + /// @param id The PDG + /// @param applyBankFilter If `true`, apply any `hipo::bank` filters (_e.g._, from iguana) + /// @see Additional methods: + /// - getByRegion + /// - getByCharge + /// - getDetParticles + std::vector getByID(int id, bool const& applyBankFilter=false) const; + + /// @brief Get the list of `region_particle`s, filtered by region + /// @returns A _copy_ of the list of corresponding `region_particle`s + /// @param ir The region + /// @param applyBankFilter If `true`, apply any `hipo::bank` filters (_e.g._, from iguana) + /// @see Additional methods: + /// - getByID + /// - getByCharge + /// - getDetParticles + std::vector getByRegion(int ir, bool const& applyBankFilter=false) const; + + /// @brief Get the list of `region_particle`s, filtered by charge + /// @returns A _copy_ of the list of corresponding `region_particle`s + /// @param ch The charge + /// @param applyBankFilter If `true`, apply any `hipo::bank` filters (_e.g._, from iguana) + /// @see Additional methods: + /// - getByID + /// - getByRegion + /// - getDetParticles + std::vector getByCharge(int ch, bool const& applyBankFilter=false) const; + + ///////////////////////////////// const std::vector& preCheckPids(); const std::vector& preCheckPidsOrCharge(); @@ -193,9 +330,14 @@ namespace clas12 { bool passPidSelect(); + /// @brief use FT-based PID, rather than the default + /// @see usingFTBased void useFTBased(){_useFTBased=true;} + + /// @returns whether or not FT-based PID is used + /// @see useFTBased + bool const usingFTBased() const { return _useFTBased; } - int getNParticles() const noexcept{return _detParticles.size();} const std::vector &getPids() const noexcept{return _pids;} bool checkTriggerBit(uint k){ @@ -257,6 +399,15 @@ namespace clas12 { _verbose=level; _reader.setVerbose(level); } + + /// @brief Set a "read action", a custom lambda function that is executed for every event within `readEvent()`. + /// @details `readEvent()` is called by methods like `clas12reader::next()` and `clas12root::HipoChain::Next()` + /// @param readEventUserAction lambda function, where its argument is a pointer to an instance of this `clas12reader` class, + /// and its `bool` return value controls whether to proceed with the event or not + void SetReadAction(std::function readEventUserAction) { + _readEventUserAction = readEventUserAction; + } + protected: @@ -345,7 +496,10 @@ namespace clas12 { bool _useFTBased{false}; bool _isOpen{false}; std::vector _addBankNames; - + + /// user-definable lambda, called in `readEvent()`; the default is a no-op returning `true` + std::function _readEventUserAction = [](clas12reader* r) {return true;}; + ///////////////////////////////DB private: int _runNo{0}; diff --git a/Clas12Banks/region_particle.h b/Clas12Banks/region_particle.h index f0882c1..7e92be0 100644 --- a/Clas12Banks/region_particle.h +++ b/Clas12Banks/region_particle.h @@ -61,6 +61,9 @@ namespace clas12 { /// i.e. how the detector banks relate to that region virtual bool sort(){ _pentry=_parts->getEntry(); + //check if allowed by particle-bank filters + _allowed = bankAllowsRow(_pentry, _parts); + _allowed_ftb = bankAllowsRow(_pentry, _ftbparts); //check for covariance matrix if(_covmat)_pcmat=_covmat->getIndex(_pentry); if(_mcpart)_pmc=_mcpart->match_to(_pentry); @@ -158,9 +161,28 @@ namespace clas12 { double getMCPhiDiff() const {return getPhi()-mc()->getPhi();} double getMCPDiff() const {return getP()-mc()->getP();} + /// @brief Whether or not this `region_particle` is allowed + /// @returns true if this `region_particle` is "allowed", _e.g._, by a HIPO bank filter, + /// which can be applied by an iguana algorithm + bool const& isAllowed() const { + if(_ftbparts==nullptr) return _allowed; + return _useFTBPid*_ftbparts->getRows() ? _allowed_ftb : _allowed; + } + //if(_parts->getCharge()) protected: + /// @returns `true` if `bank` contains row `row` in its `hipo::bank::rowlist` + /// @param row the row to look up + /// @param bank the bank to check + bool const bankAllowsRow(int const& row, hipo::bank const* bank) const { + // NOTE: `hipo::bank::getRowList()` returns the list of rows which are _allowed_ by `hipo::bank::filter()`; for example, + // `hipo::bank::filter` is used by iguana fiducial cut algorithms to select particles from `REC::Particle` + if(bank!=nullptr && bank->getRows()>0) + return std::find(bank->getRowList().begin(), bank->getRowList().end(), row) != bank->getRowList().end(); + return false; + } + par_ptr _parts={nullptr}; ftbpar_ptr _ftbparts={nullptr}; covmat_ptr _covmat={nullptr}; @@ -183,6 +205,10 @@ namespace clas12 { short _pcmat=-1; short _region=-1; short _useFTBPid=0; + + // filter status (from `hipo::bank::filter`) + bool _allowed{true}; // allowed by `_parts`, i.e., `REC::Particle` + bool _allowed_ftb{true}; // allowed by `_ftbparts`, i.e., `RECFT::Particle` }; //pointer "typedef" using region_part_ptr=clas12::region_particle*; diff --git a/Clas12Root/HipoChain.h b/Clas12Root/HipoChain.h index 6a5c751..1d46386 100644 --- a/Clas12Root/HipoChain.h +++ b/Clas12Root/HipoChain.h @@ -15,6 +15,8 @@ namespace clas12root { + /// @brief Clas12root reader of multiple HIPO files + /// @see `clas12::clas12reader` class HipoChain : public TNamed { diff --git a/README.md b/README.md index 0c12560..54f0620 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,7 @@ Data Analysis Tools for hipo4 data format. +[:star: **Click here for the Source Code Documentation (doxygen)** :star:](https://jeffersonlab.github.io/clas12root) Examples are given for running in interactive ROOT sessions and ROOT-Jupyter notebooks. @@ -578,63 +579,9 @@ Note the event number is just its position in the file, not the DST RUN::Config: ## Ex 11 Iguana interface To run iguana routines in clas12root you should first set the environment to point to an -installed version of iguana, by setting the IGUANA variable to the INSTALLATION directory. - -There are 2 methods of using Iguana. For speed and simplicity both use just the iguana -action function and will not operate on the underlying banks structures. -In the first method, which can use any iguan algorithm, you just directly create the iguana -algorithm and run it yourself. See Ex11_Iguana_example_01_bank_rows.C - -The second method is though a higher level interface which simplifies the usage, but is not -garuanteed to be able to use all algorithms. In these cases you may revert to the first method -to apply any additional algorithms. The interface for this is in the $CLAS12ROOT/iguana/ directory -and an example is available at Ex11_Iguana_MomentumCorrection.C. - -These example use a HipoChain and you will need to set file paths yourself etc. -Using clas12root:region_particles means that information for -alogithms is readily available and so we can just pass these particle objects -into the action functions. - -Highlighting the iguana parts for method 2 : - - //clas12root-iguana interface - clas12root::Iguana ig; - - //choose some algorithms to apply - ig.GetTransformers().Use("clas12::MomentumCorrection"); - ig.GetFilters().Use("clas12::ZVertexFilter"); - ig.GetCreators().Use("physics::InclusiveKinematics"); - - ig.SetOptionAll("log", "debug"); - ig.Start(); - - - ... - // get particles by type - // note we applied a filter to ensure size of all ==1 - auto electron=c12->getByID(11)[0]; - auto pip=c12->getByID(211)[0]; - auto pim=c12->getByID(-211)[0]; - - ///Now do some iguana analysis!! - - //filter on z-vertices of the particles - //note I can pass a vector of clas12root particles - if( !(ig.GetFilters().doZVertexFilter({electron,pip,pim})) ) { - continue; - } - - - //correct momentum and get 4-vectors - //I pass a vector of clas12root particles and LorentzVectors - ig.GetTransformers().doMomentumCorrections({electron,pip,pim},{&p4el,&p4pip,&p4pim}); - - //calculate inclusive kinematics - //use original bank info - auto kine = ig.GetCreators().doInclusiveKinematics(electron); - //use momentum corrected momentum - auto corrkine = ig.GetCreators().doInclusiveKinematics(p4el); - - ... - +installed version of iguana, by setting the `IGUANA` variable to the INSTALLATION directory. +Clas12root does _not_ depend on Iguana, but running `clas12root` will load its libraries for you if +you have set the `IGUANA` environment variable. +For usage of Iguana with `clas12root`, see the examples, such as +- [`Ex11_Iguana.C`](/RunRoot/Ex11_Iguana.C) diff --git a/RunRoot/Ex11_Iguana.C b/RunRoot/Ex11_Iguana.C new file mode 100644 index 0000000..8dbd36e --- /dev/null +++ b/RunRoot/Ex11_Iguana.C @@ -0,0 +1,362 @@ +/** + * Example demonstrating how to read a HIPO file with clas12root, + * and use Iguana algorithms to process its banks + * + * NOTE: this requires Iguana version 1.0.0 or newer + * + * USAGE: + * ``` + * clas12root RunRoot/Ex11_Iguana.C+ --in=hipo_file_1.hipo --in=hipo_file_2.hipo + * ``` + */ + +// ROOT +#include +#include +#include +#include +#include +#include + +// clas12root +#include + +// include iguana algorithms +#include +#include +#include +#include +#include + +void Ex11_Iguana() { + + // parse arguments, which should be HIPO filename(s) prefixed with `--in=`; add them to a `HipoChain` + clas12root::HipoChain chain; + for(int i=2; iArgc(); i++) { + TString inputFile = gApplication->Argv(2); + inputFile(TRegexp("^--in=")) = ""; + std::cout << "reading file " << inputFile << std::endl; + chain.Add(inputFile); + chain.SetReaderTags({0}); // read tag-0 events only + } + if(chain.GetNFiles() == 0) { + std::cerr << " *** please provide HIPO file name(s)" << std::endl; + exit(1); + } + + ////////////////////////////////////////////////////////////////////////////////// + + // histograms + int const n_bins = 100; + auto* Q2_vs_x = new TH2D("Q2_vs_x", "Q^{2} vs. x;x;Q^{2} [GeV^{2}]", n_bins, 0, 1, n_bins, 0, 12); + auto* Q2_vs_W = new TH2D("Q2_vs_W", "Q^{2} vs. W;W [GeV];Q^{2} [GeV^{2}]", n_bins, 0, 5, n_bins, 0, 12); + auto* y_dist = new TH1D("y_dist", "y distribution;y", n_bins, 0, 1); + auto* vz_dist = new TH1D("vz_dist", "electron v_{z};v_{z} [cm]", n_bins, -30, 30); + auto* deltaP_vs_P = new TH2D("deltaP_vs_P", "electron momentum correction;p_{meas} [GeV];p_{corr}-p_{meas} [GeV]", 10, 0, 12, n_bins, -0.2, 0.2); + // histogram styles + for(auto hist : {y_dist, vz_dist}) { + hist->SetLineColor(kAzure); + hist->SetFillColor(kAzure); + } + + ////////////////////////////////////////////////////////////////////////////////// + + // create iguana algorithms + iguana::clas12::ZVertexFilter algo_vz_filter; // filter the z-vertex (a filter algorithm) + iguana::clas12::SectorFinder algo_sec_finder; // get the sector for each particle (a creator algorithm) + iguana::clas12::rga::FiducialFilterPass2 algo_fidu_filter; // fiducial cuts (a filter algorithm) + iguana::clas12::rga::MomentumCorrection algo_mom_corr; // momentum corrections (a transformer algorithm) + iguana::physics::InclusiveKinematics algo_inc_kin; // calculate inclusive kinematics (a creator algorithm) + + /* + * configure iguana algorithms here (see iguana documentation) + */ + // algo_sec_finder.SetOption("log", "trace"); // for example, more detailed logging + + // start iguana algorithms, after configuring them + algo_vz_filter.Start(); + algo_sec_finder.Start(); + algo_fidu_filter.Start(); + algo_mom_corr.Start(); + algo_inc_kin.Start(); + + // create bank objects, which creator-type algorithms will populate + // - by creating them here, this macro is the owner of these bank objects, but + // iguana algorithms' Run functions will handle them by reference + auto iguana_bank_sec = algo_sec_finder.GetCreatedBank(); + auto iguana_bank_inc_kin = algo_inc_kin.GetCreatedBank(); + + // let's also create a "cache" of particle momentum magnitudes, since we want to check the effects + // of momentum corrections + std::vector p_measured; + + ////////////////////////////////////////////////////////////////////////////////// + + // define a lambda function that processes HIPO banks, in particular, with iguana + // + // - this function will be executed by `clas12reader` just after each event's `hipo::bank` + // objects are filled, and before `clas12reader` further processes the banks; this happens + // within functions such as: + // - `clas12reader::next()` + // - `HipoChain::Next()` + // + // - behavior from iguana algorithm types: + // - filter algorithms typically filter certain bank rows, such as a z-vertex + // filter will choose rows of `REC::Particle`; accessors such as `clas12reader::getDetParticles` + // and `clas12reader::getByID` can be used to access the filtered (or unfiltered) particles + // - transformer algorithms modify values within banks, for example, a momentum corrections algorithm + // corrects momenta of particles in the `REC::Particle` bank; then clas12root objects which are based + // on `REC::Particle`, such as `region_particle`, will simply have the corrected momenta + // - creator algorithms create new banks that are not known by clas12root; they are `hipo::bank` + // objects and you may use them as any other bank directly + // + // - this example macro owns algorithm instances and created bank objects, therefore the lambda + // function must capture them by reference + // + // - the lambda function's parameter `cr` is a pointer, which is a `clas12reader` instance, whether you are + // using `clas12reader` directly (see `RunRoot/Ex1_CLAS12Reader.C`) or + // using `HipoChain` (see `RunRoot/Ex1_CLAS12ReaderChain.C`) + // - when you call functions like `clas12reader::next()` or `HipoChain::Next()`, the + // corresponding `clas12reader` instance will be passed into this lambda function + // - use `cr` to access DST banks from the `clas12reader`, using `cr->getBankName()` methods, + // where "BankName" is the name of the bank without the colons (`::`); for example, use + // `cr->getRECFTParticle` to access the `RECFT::Particle` bank + // + // - the lambda function's return value type must be `bool`: + // - return `true` if you want the event to proceed + // - return `false` if you want to skip the event + // + auto iguana_action = [ + // capture the algorithm instances; use the ampersand (`&`) to capture them by reference + &algo_vz_filter, + &algo_sec_finder, + &algo_fidu_filter, + &algo_mom_corr, + &algo_inc_kin, + // capture the iguana-created banks, again by reference + &iguana_bank_sec, + &iguana_bank_inc_kin, + // also capture our `p_measured` cache, since we want to populate it + &p_measured + ](clas12::clas12reader* cr) + { + // before anything, let's "cache" our original momentum magnitudes: loop over `REC::Particle` + // and store |p| for each particle + p_measured.clear(); + for(auto const& row : cr->getRECParticle().getRowList()) + p_measured.push_back(std::hypot( + cr->getRECParticle().getFloat("px", row), + cr->getRECParticle().getFloat("py", row), + cr->getRECParticle().getFloat("pz", row))); + + // call Iguana `Run` functions + // - think carefully about the ordering; for example, do you correct momenta before or after + // applying a filter which depends on momenta? + // - these `Run` functions take `hipo::bank` objects as parameters (technically, lvalue references) + // - some bank objects may be updated, depending on the algorithm type + // - filters tend to filter certain bank's rows + // - transformers tend to mutate certain values + // - creators create new banks, i.e., the `iguana_bank_*` objects will be populated with data + // - DST banks are read from the `clas12reader` instance, `cr` + // - created banks must exist so they can be filled; this is why we created them beforehand + // - `Run` functions return boolean, which we can use to skip an event at any time, i.e., use + // the pattern `if(!algo.Run(...)) return false;` + + // start with the z-vertex filter; it returns false if no electrons pass the filter + if(!algo_vz_filter.Run( + cr->getRECParticle(), // REC::Particle + cr->getRUNconfig() // RUN::config + )) return false; + + // next, let's apply the fiducial cuts; we'll skip FT cuts for this example, but take a look + // at the algorithm's documentation for other `Run` functions that use FT data + if(!algo_fidu_filter.Run( + cr->getRECParticle(), // REC::Particle + cr->getRUNconfig(), // RUN::config + cr->getRECCalorimeter(), // REC::Calorimeter + cr->getRECTraj() // REC::Traj + )) return false; + + // momentum corrections require sector information, so call the sector finder first; + // note that the momentum correction algorithm will alter the `px,py,pz` of certain + // particles, which is why we filled `p_measured` first above + if(!algo_sec_finder.Run( + cr->getRECParticle(), // REC::Particle + cr->getRECTrack(), // REC::Track + cr->getRECCalorimeter(), // REC::Calorimeter + cr->getRECScintillator(), // REC::Scintillator + iguana_bank_sec // result of iguana::clas12::SectorFinder ----. + )) return false; // | + // | + if(!algo_mom_corr.Run( // | + cr->getRECParticle(), // REC::Particle | + iguana_bank_sec, // result of iguana::clas12::SectorFinder <---' + cr->getRUNconfig() // RUN::config + )) return false; + + // finally, calculate inclusive kinematics + if(!algo_inc_kin.Run( + cr->getRECParticle(), // REC::Particle + cr->getRUNconfig(), // RUN::config + iguana_bank_inc_kin // result of iguana::physics::InclusiveKinematics + )) return false; + + // all algorithms are done, return true to keep the event + // (otherwise return false if you have some reason to not keep the event) + return true; + }; + + // attach the iguana-running lambda function to the `clas12reader` event reader + // - `chain.Next()` will call `iguana_action` internally, just after the HIPO banks have been + // read, but before any bank-dependent objects are created, such as `region_particle` + chain.GetC12Reader()->SetReadAction(iguana_action); + + ////////////////////////////////////////////////////////////////////////////////// + + // now get a pointer to the internal `clas12reader` + // this will point to the correct place when file changes + auto const& c12 = chain.C12ref(); + + // loop over all events + for(int numEvents=0; chain.Next(); numEvents++) { + + // let's be verbose for the first few events, to demonstrate what Iguana did + if(numEvents < 30) { + std::cout << "===== EVENT " << c12->runconfig()->getEvent() << "===========\n"; + + // print banks + // ----------- + std::cout << "------- FULL PARTICLE BANK -------\n"; + c12->getRECParticle().show(true); // use `true`, otherwise only filter-allowed rows are printed + std::cout << "----- FILTERED PARTICLE BANK -----\n"; + c12->getRECParticle().show(); + std::cout << "---------- IGUANA BANKS ----------\n"; + iguana_bank_sec.show(); + iguana_bank_inc_kin.show(); + std::cout << "----------------------------------\n"; + + // Accessing bank rows and filtering + // --------------------------------- + // if you want to loop over filtered rows, use `hipo::bank::getRowList()`; otherwise, + // just use the usual `for` loop from `0` up to `bank->getRows()` + std::cout << "REC::Particle filter-allowed rows:"; + for(const auto& row : c12->getRECParticle().getRowList()) + std::cout << " " << row; + std::cout << " (out of " << c12->getRECParticle().getRows() << " rows total)\n"; + + // Get particles by type + // --------------------- + // NOTE: to make sure that only the particles which passed Iguana filters, + // set the additional boolean argument to `true`, otherwise you will get ALL + // the particles; this may done in any `region_particle` list accessor, such as: + // - `clas12reader::getByID` + // - `clas12reader::getByCharge` + // - `clas12reader::getByRegion` + // - `clas12reader::getDetParticles` + std::cout << "electrons allowed by Iguana:"; + for(auto const* electron : c12->getByID(11, true)) // loops over electrons which passed the filter(s) + std::cout << " " << electron->getIndex(); + std::cout << "\n"; + + // alternatively, you could use `region_particle::isAllowed()` to filter as needed + std::cout << "electrons filtered out by Iguana:"; + for(auto const* electron : c12->getByID(11)) { // loops over ALL electrons + if( ! electron->isAllowed()) { // selects electrons which were filtered OUT + std::cout << " " << electron->getIndex(); + } + } + std::cout << "\n"; + } + else if(numEvents % 10000 == 0) + std::cout << "read " << numEvents << " events\n"; + + // from here, refer to other examples on how to proceed; + // in this example, we'll fill the histograms + + // inclusive kinematics; the created bank has only 1 row + Q2_vs_x->Fill(iguana_bank_inc_kin.getDouble("x", 0), iguana_bank_inc_kin.getDouble("Q2", 0)); + Q2_vs_W->Fill(iguana_bank_inc_kin.getDouble("W", 0), iguana_bank_inc_kin.getDouble("Q2", 0)); + y_dist->Fill(iguana_bank_inc_kin.getDouble("y", 0)); + + // scattered electron pindex + auto pindex_ele = iguana_bank_inc_kin.getShort("pindex", 0); + + // loop over electrons, and choose the scattered electron + // note: it would be faster to use `c12->getRECParticle()` and choose row `pindex_ele`, but here + // we want to demonstrate how to use `c12->getByID` to access the same information + for(auto const* electron : c12->getByID(11, true)) { // loops over electrons which passed the filter(s) + if(electron->getIndex() == pindex_ele) { // choose only the scattered electron + // electron vertex + auto vz = electron->par()->getVz(); + vz_dist->Fill(vz); + // electron momentum, which was corrected by Iguana's momentum corrections; + // notice we use the `p_measured` cache here, to plot the momentum correction amount + auto p_corrected = electron->getP(); + deltaP_vs_P->Fill(p_corrected, p_corrected - p_measured.at(pindex_ele)); + // note: alternatively, we could have just gotten these values directly + // from the `REC::Particle` bank; let's cross check the value to prove it: + auto vz_from_bank = c12->getRECParticle().getFloat("vz", pindex_ele); + auto p_from_bank = std::hypot( + c12->getRECParticle().getFloat("px", pindex_ele), + c12->getRECParticle().getFloat("py", pindex_ele), + c12->getRECParticle().getFloat("pz", pindex_ele)); + if(std::abs(vz - vz_from_bank) > 1e-4) + throw std::runtime_error("vz cross check failed: " + std::to_string(vz) + " vs. " + std::to_string(vz_from_bank)); + if(std::abs(p_corrected - p_from_bank) > 1e-4) + throw std::runtime_error("|p| cross check failed: " + std::to_string(p_corrected) + " vs. " + std::to_string(p_from_bank)); + // exit the loop over electrons + break; + } + } + + } // end event loop + + ////////////////////////////////////////////////////////////////////////////////// + + // stop the iguana algorithms, now that the event loop is done + algo_vz_filter.Stop(); + algo_sec_finder.Stop(); + algo_fidu_filter.Stop(); + algo_mom_corr.Stop(); + algo_inc_kin.Stop(); + + // draw the plots + int const n_rows = 2; + int const n_cols = 3; + auto canv = new TCanvas("canv", "canv", n_cols * 800, n_rows * 600); + canv->Divide(n_cols, n_rows); + for(int pad_num = 1; pad_num <= n_rows * n_cols; pad_num++) { + auto pad = canv->GetPad(pad_num); + pad->cd(); + pad->SetGrid(1, 1); + pad->SetLeftMargin(0.12); + pad->SetRightMargin(0.12); + pad->SetBottomMargin(0.12); + switch(pad_num) { + case 1: + pad->SetLogz(); + Q2_vs_x->Draw("colz"); + break; + case 2: + pad->SetLogz(); + Q2_vs_W->Draw("colz"); + break; + case 3: + y_dist->Draw(); + break; + case 4: + pad->SetLogy(); + vz_dist->Draw(); + break; + case 5: + pad->SetLogz(); + deltaP_vs_P->Draw("colz"); + auto prof = deltaP_vs_P->ProfileX("_pfx", 1, -1, "s"); + prof->SetLineColor(kBlack); + prof->SetLineWidth(5); + prof->Draw("same"); + break; + } + } + +} diff --git a/RunRoot/Ex11_Iguana_MomentumCorrection.C b/RunRoot/Ex11_Iguana_MomentumCorrection.C deleted file mode 100644 index 260a2b7..0000000 --- a/RunRoot/Ex11_Iguana_MomentumCorrection.C +++ /dev/null @@ -1,141 +0,0 @@ -#include "clas12reader.h" -#include "region_particle.h" -#include "HipoChain.h" -#include "Iguana.h" -#include -#include -#include -#include -#include -#include -#include -#include - -void Ex11_Iguana_MomentumCorrection(){ - - using FourVector = ROOT::Math::PxPyPzMVector; - // clas12databases::SetRCDBRootConnection("rcdb.root"); - gBenchmark->Start("iguana"); - //create hipo chain for looping over events - clas12root::HipoChain chain; - chain.Add("/home/dglazier/Jlab/clas12/data/hipo/DVPipPimP_006733.hipo"); - chain.SetReaderTags({0}); //create clas12reader with just tag 0 events - //add some filters - auto config_c12=chain.GetC12Reader(); - config_c12->addAtLeastPid(11,1); //exactly 1 electron - config_c12->addAtLeastPid(211,1); //exactly 1 pi+ - config_c12->addAtLeastPid(-211,1); //exactly 1 pi- - config_c12->ignoreBank("REC::CovMat"); - - - // create the chosen algorithms - clas12root::Iguana ig; - ig.SetClas12(chain.C12ptr());//connect to clas12reader - ig.GetTransformers().Use("clas12::MomentumCorrection"); - ig.GetTransformers().Use("clas12::FTEnergyCorrection"); - ig.GetFilters().Use("clas12::ZVertexFilter"); - ig.GetCreators().Use("physics::InclusiveKinematics"); - - ig.SetOptionAll("log", "debug"); - ig.Start(); - - - //decalre 4-vector objects - auto db=TDatabasePDG::Instance(); - FourVector p4beam(0,0,10.6,10.6); - FourVector p4target(0,0,0,db->GetParticle(2212)->Mass()); - FourVector p4el(0,0,0,db->GetParticle(11)->Mass()); - FourVector p4pip(0,0,0,db->GetParticle(211)->Mass()); - FourVector p4pim(0,0,0,db->GetParticle(-211)->Mass()); - TLorentzVector tp4el(0,0,0,db->GetParticle(11)->Mass()); - - TH1D hQ2{"Q2","Q^{2}",100,0,10}; - TH1D hx{"x","x",100,0,1}; - TH1D hy{"y","y",100,0,1}; - TH1D hW{"W","W",100,0,6}; - TH1D hcQ2{"cQ2","corrected Q^{2}",100,0,10}; - TH1D hcx{"cx","corrected x",100,0,1}; - TH1D hcy{"cy","corrected y",100,0,1}; - TH1D hcW{"cW","corrected W",100,0,6}; - - - //get run clas12reader - auto& c12=chain.C12ref(); - while ( chain.Next() ){ - ig.PrepareEvent(c12->getRunNumber()); - - //auto ebeam = c12->rcdb()->current().beam_energy/1000; - //p4beam.SetXYZT(0,0,ebeam,ebeam); //approx. mass =0 - //filter on z-vertices of the particles - //filter on fiducial cuts - //filter on PhotonGBT - ig.GetFilters().doAllFilters(); - - - //Check if particles still exist after filters - //and if they do assign them - auto getFirstParticle=[&c12](short pdg)->clas12::region_particle*{ - auto particles=c12->getByID(pdg); - if(particles.size()>0){ - return particles[0]; - } - return nullptr; - }; - - auto electron=getFirstParticle(11); - if(electron==nullptr) continue; - auto pip=getFirstParticle(211); - if(pip==nullptr) continue; - auto pim=getFirstParticle(-211); - if(pim==nullptr) continue; - - //correct momentum and get 4-vectors - ig.GetTransformers().doAllCorrections({electron,pip,pim},{&p4el,&p4pip,&p4pim}); - - //calculate inclusive kinematics - auto kine = ig.GetCreators().doInclusiveKinematics(electron); - auto corrkine = ig.GetCreators().doInclusiveKinematics(p4el); - - hQ2.Fill(kine.Q2); - hx.Fill(kine.x); - hy.Fill(kine.y); - hW.Fill(kine.W); - - hcQ2.Fill(corrkine.Q2); - hcx.Fill(corrkine.x); - hcy.Fill(corrkine.y); - hcW.Fill(corrkine.W); - - //if you unfortunately prefer TLorentzVector, there is a transform function to help. - // ig.GoodToBad(p4el,tp4el); - // cout<<"check good to bad "<Stop("iguana"); - gBenchmark->Print("iguana"); - - auto* can = new TCanvas(); - can->Divide(2,2); - - can->cd(1); - hQ2.DrawCopy(); - hcQ2.SetLineColor(2); - hcQ2.DrawCopy("same"); - - can->cd(2); - hW.DrawCopy(); - hcW.SetLineColor(2); - hcW.DrawCopy("same"); - - can->cd(3); - hx.DrawCopy(); - hcx.SetLineColor(2); - hcx.DrawCopy("same"); - - can->cd(4); - hy.DrawCopy(); - hcy.SetLineColor(2); - hcy.DrawCopy("same"); - - std::cout<<"analysed "< -#include - -void Ex11_Iguana_example_01_bank_rows(){ - // clas12databases::SetRCDBRootConnection("rcdb.root"); - - clas12root::HipoChain chain; - chain.Add("/home/dglazier/Jlab/clas12/data/hipo/DVPipPimP_006733.hipo"); - chain.SetReaderTags({0}); //create clas12reader with just tag 0 events - auto config_c12=chain.GetC12Reader(); - config_c12->addExactPid(11,1); //exactly 1 electron - config_c12->addExactPid(211,1); //exactly 1 pi+ - config_c12->addExactPid(-211,1); //exactly 1 pi- - - // create the algorithms - iguana::clas12::MomentumCorrection algo_momentum_correction; - // set log levels - algo_momentum_correction.SetOption("log", "debug"); - // start the algorithms - algo_momentum_correction.Start(); - - //get run clas12reader - auto& c12=chain.C12ref(); - int counter =0; - - //decalre 4-vector objects - auto db=TDatabasePDG::Instance(); - TLorentzVector p4beam(0,0,10.6,10.6); - TLorentzVector p4target(0,0,0,db->GetParticle(2212)->Mass()); - TLorentzVector p4el(0,0,0,db->GetParticle(11)->Mass()); - TLorentzVector p4pip(0,0,0,db->GetParticle(211)->Mass()); - TLorentzVector p4pim(0,0,0,db->GetParticle(-211)->Mass()); - - //define mometnum correction function taking a region_particle and a lorentzvector - auto momCorrection = [&algo_momentum_correction,&c12](region_particle *p, TLorentzVector& p4){ - - auto [px, py, pz] = algo_momentum_correction.Transform(p->getPx(),p->getPy(),p->getPz(),p->getSector(),p->getPid(),c12->runconfig()->getTorus()); - - p4.SetXYZM(px,py,pz,p4.M()); - }; - - //loop over all events - - while (chain.Next()){ - //auto ebeam = c12->rcdb()->current().beam_energy/1000; - // p4beam.SetXYZT(0,0,ebeam,ebeam); //approx. mass =0 - // get particles by type - // note we applied a filter to ensure size of all ==1 - auto electron=c12->getByID(11)[0]; - auto pip=c12->getByID(211)[0]; - auto pim=c12->getByID(-211)[0]; - - momCorrection(electron,p4el); - momCorrection(pip,p4pip); - momCorrection(pim,p4pim); - cout<<"el px "<getPx()<<" to "<getPy()<<" to "<getPz()<<" to "<Getenv("IGUANA"); gInterpreter->AddIncludePath(IGUANA+"/include"); - - gROOT->ProcessLine("#include "); - gROOT->ProcessLine("#include "); + gROOT->ProcessLine("#include "); - + gSystem->Load("$IGUANA/lib/libIguanaServices.so"); gSystem->Load("$IGUANA/lib/libIguanaAlgorithms.so"); - gSystem->Load("$IGUANA/lib/libIguanaValidators.so"); - - //clas12root iguana stuff - TString CLAS12ROOT = gSystem->Getenv("CLAS12ROOT"); - gInterpreter->AddIncludePath(Form("%s/iguana",CLAS12ROOT.Data())); - gROOT->ProcessLine(Form("#include \"%s/iguana/IguanaAlgo.h\" ",CLAS12ROOT.Data())); - gROOT->ProcessLine(Form("#include \"%s/iguana/Transformers.h\" ",CLAS12ROOT.Data())); - gROOT->ProcessLine(Form("#include \"%s/iguana/Filters.h\" ",CLAS12ROOT.Data())); - gROOT->ProcessLine(Form("#include \"%s/iguana/Creators.h\" ",CLAS12ROOT.Data())); - gROOT->ProcessLine(Form("#include \"%s/iguana/Iguana.h\" ",CLAS12ROOT.Data())); - } diff --git a/Doxyfile b/docs/doxygen/Doxyfile similarity index 98% rename from Doxyfile rename to docs/doxygen/Doxyfile index 4982ccb..12cebcd 100644 --- a/Doxyfile +++ b/docs/doxygen/Doxyfile @@ -58,7 +58,7 @@ PROJECT_LOGO = # entered, it will be relative to the location where doxygen was started. If # left blank the current directory will be used. -OUTPUT_DIRECTORY = "../clas12rootDoxy" +OUTPUT_DIRECTORY = share/doc/clas12root # If the CREATE_SUBDIRS tag is set to YES then doxygen will create 4096 sub- # directories (in 2 levels) under the output directory of each output format and @@ -239,6 +239,8 @@ TAB_SIZE = 4 # newlines. ALIASES = +ALIASES += bank_accessor{1}="@brief get the `\1` bank^^@xrefitem bank_acc \"Click here for the full list of bank accessors\" \"Bank Accessors\" `\1`^^@return a reference to the `\1` bank^^" +ALIASES += bank_ptr_accessor{2|}="@brief get a pointer to the `\1` bank^^@return a pointer to the `\1` bank^^@see Bank reference accessor `\2`^^" # This tag can be used to specify a number of word-keyword mappings (TCL only). # A mapping has the form "name=value". For example adding "class=itcl::class" @@ -578,7 +580,7 @@ SORT_MEMBER_DOCS = YES # this will also influence the order of the classes in the class list. # The default value is: NO. -SORT_BRIEF_DOCS = NO +SORT_BRIEF_DOCS = YES # If the SORT_MEMBERS_CTORS_1ST tag is set to YES then doxygen will sort the # (brief and detailed) documentation of class members so that constructors and @@ -590,7 +592,7 @@ SORT_BRIEF_DOCS = NO # detailed member documentation. # The default value is: NO. -SORT_MEMBERS_CTORS_1ST = NO +SORT_MEMBERS_CTORS_1ST = YES # If the SORT_GROUP_NAMES tag is set to YES then doxygen will sort the hierarchy # of group names into alphabetical order. If set to NO the group names will @@ -703,7 +705,7 @@ FILE_VERSION_FILTER = # DoxygenLayout.xml, doxygen will parse it automatically even if the LAYOUT_FILE # tag is left empty. -LAYOUT_FILE = +LAYOUT_FILE = docs/doxygen/DoxygenLayout.xml # The CITE_BIB_FILES tag can be used to specify one or more bib files containing # the reference definitions. This must be a list of .bib files. The .bib @@ -790,7 +792,10 @@ WARN_LOGFILE = # spaces. See also FILE_PATTERNS and EXTENSION_MAPPING # Note: If this tag is empty the current directory is searched. -INPUT = "hipo4" "Clas12Banks" "Clas12Root" +INPUT = hipo4 \ + Clas12Banks \ + Clas12Root \ + docs/doxygen/mainpage.md # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses @@ -982,7 +987,7 @@ FILTER_SOURCE_PATTERNS = # (index.html). This can be useful if you have a project on for instance GitHub # and want to reuse the introduction page also for the doxygen output. -USE_MDFILE_AS_MAINPAGE = +USE_MDFILE_AS_MAINPAGE = docs/doxygen/mainpage.md #--------------------------------------------------------------------------- # Configuration options related to source browsing @@ -1204,6 +1209,19 @@ HTML_EXTRA_STYLESHEET = HTML_EXTRA_FILES = +# The HTML_COLORSTYLE tag can be used to specify if the generated HTML output +# should be rendered with a dark or light theme. +# Possible values are: LIGHT always generate light mode output, DARK always +# generate dark mode output, AUTO_LIGHT automatically set the mode according to +# the user preference, use light mode if no preference is set (the default), +# AUTO_DARK automatically set the mode according to the user preference, use +# dark mode if no preference is set and TOGGLE allow to user to switch between +# light and dark mode via a button. +# The default value is: AUTO_LIGHT. +# This tag requires that the tag GENERATE_HTML is set to YES. + +HTML_COLORSTYLE = TOGGLE + # The HTML_COLORSTYLE_HUE tag controls the color of the HTML output. Doxygen # will adjust the colors in the style sheet and background images according to # this color. Hue is specified as an angle on a colorwheel, see @@ -1213,7 +1231,7 @@ HTML_EXTRA_FILES = # Minimum value: 0, maximum value: 359, default value: 220. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_COLORSTYLE_HUE = 220 +HTML_COLORSTYLE_HUE = 180 # The HTML_COLORSTYLE_SAT tag controls the purity (or saturation) of the colors # in the HTML output. For a value of 0 the output will use grayscales only. A @@ -1478,7 +1496,7 @@ DISABLE_INDEX = NO # The default value is: NO. # This tag requires that the tag GENERATE_HTML is set to YES. -GENERATE_TREEVIEW = NO +GENERATE_TREEVIEW = YES # The ENUM_VALUES_PER_LINE tag can be used to set the number of enum values that # doxygen will group on one line in the generated HTML documentation. @@ -1666,7 +1684,7 @@ EXTRA_SEARCH_MAPPINGS = # If the GENERATE_LATEX tag is set to YES, doxygen will generate LaTeX output. # The default value is: YES. -GENERATE_LATEX = YES +GENERATE_LATEX = NO # The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. If a # relative path is entered the value of OUTPUT_DIRECTORY will be put in front of diff --git a/docs/doxygen/DoxygenLayout.xml b/docs/doxygen/DoxygenLayout.xml new file mode 100644 index 0000000..4c64185 --- /dev/null +++ b/docs/doxygen/DoxygenLayout.xml @@ -0,0 +1,269 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/doxygen/generate.sh b/docs/doxygen/generate.sh new file mode 100755 index 0000000..45abd70 --- /dev/null +++ b/docs/doxygen/generate.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env bash +set -euo pipefail +cd $(dirname ${BASH_SOURCE[0]:-$0})/../.. +prefix=share/doc +mkdir -p $prefix +doxygen docs/doxygen/Doxyfile +echo """ +Generated documentation webpage: + file://$(pwd)/$(find $prefix -type f -name index.html) +""" diff --git a/docs/doxygen/mainpage.md b/docs/doxygen/mainpage.md new file mode 100644 index 0000000..32ffc18 --- /dev/null +++ b/docs/doxygen/mainpage.md @@ -0,0 +1,15 @@ +# Clas12root Code Documentation + +@todo Add more user documentation here! Populating CLAS12 Software documentation in general is an open service work project... + +## Commonly Used Classes + +- For reading HIPO files: + - `clas12::clas12reader`: for reading a single HIPO file + - `clas12root::HipoChain`: for reading multiple HIPO files +- For particles: + - `clas12::region_particle`: base class region particle + - `clas12::region_band`: for BAND + - `clas12::region_cdet`: for Central Detector + - `clas12::region_fdet`: for Forward Detector + - `clas12::region_ft`: for Forward Tagger (FT) diff --git a/iguana/Creators.h b/iguana/Creators.h deleted file mode 100644 index 3a66e6a..0000000 --- a/iguana/Creators.h +++ /dev/null @@ -1,55 +0,0 @@ -#pragma once - -#include "IguanaAlgo.h" -#include - -namespace clas12root{ - - - class Creators : public IguanaAlgo { - - public: - - iguana::physics::InclusiveKinematics* InclusiveKinematics(){return _incKin.get();} - - - void Use(const string& name){ - //Initialise the inclusive kinemtics algorithm - if(name=="physics::InclusiveKinematics"){ - _incKin.reset(new iguana::physics::InclusiveKinematics ); - AddAlgo(_incKin.get()); - } - - - } - - //Inclusive Kinematics - iguana::physics::InclusiveKinematicsVars doInclusiveKinematics(const clas12::region_particle *electron){ - return _incKin->ComputeFromLepton(electron->getPx(),electron->getPy(),electron->getPz(),_concurrent_key); //0 => concurrent key - } - - iguana::physics::InclusiveKinematicsVars doInclusiveKinematics(const FourVector& electron){ - return _incKin->ComputeFromLepton(electron.X(),electron.Y(),electron.Z(),_concurrent_key); //0 => concurrent key - } - iguana::physics::InclusiveKinematicsVars doInclusiveKinematics(const TLorentzVector& electron){ - return _incKin->ComputeFromLepton(electron.X(),electron.Y(),electron.Z(),_concurrent_key); //0 => concurrent key - } - //Reload - void Reload() const{ - auto run = C12()->runconfig()->getRun(); - if(_runNb==run) return; - _runNb=run; - - if(_incKin.get())_concurrent_key=_incKin->PrepareEvent(run); - } - - private: - - std::unique_ptr _incKin; - - mutable Long64_t _runNb =-1; - mutable int _concurrent_key=0; - - }; - -} diff --git a/iguana/Filters.h b/iguana/Filters.h deleted file mode 100644 index 6d9d5eb..0000000 --- a/iguana/Filters.h +++ /dev/null @@ -1,151 +0,0 @@ -#pragma once - -#include "IguanaAlgo.h" -#include -//#include -//#include -#include - -namespace clas12root{ - - - class Filters : public IguanaAlgo { - - public: - - iguana::clas12::ZVertexFilter* ZVertexFilter(){return _zVer.get();} - // iguana::clas12::FiducialFilter* FiducialFilter(){return _Fid.get();} - //iguana::clas12::PhotonGBTFilter* PhotonGBTFilter(){return _PhotGBT.get();} - - void Use(const string& name){ - - //Initialise the z-vertex filter algorithm - if(name=="clas12::ZVertexFilter"){ - _zVer.reset(new iguana::clas12::ZVertexFilter); - AddAlgo(_zVer.get()); - } - // if(name=="clas12::FiducialFilter"){ - // _Fid.reset(new iguana::clas12::FiducialFilter); - // AddAlgo(_Fid.get()); - // } - /* - if(name=="clas12::PhotonGBTFilter"){ - _PhotGBT.reset(new iguana::clas12::PhotonGBTFilter); - AddAlgo(_PhotGBT.get()); - } - */ - } - - //interface to vector transforms - void doFilters(std::function< bool ( const clas12::region_particle *p ) > fun) const{ - // bool filt = true; - // for(auto* p:parts) filt*=fun(p); - // return filt; - - //get the clas12root region particles vector for the event - auto& particles = C12()->getDetParticles(); - //save particles to delete - //(cannot to it actively when looping over same vector) - std::vector to_delete; - for(auto* pa:particles){ - if(fun(pa)==false){//failed filter must remove - to_delete.push_back(pa); - } - } - //remove failed particles from clas12root vector - if(to_delete.empty()==false){ - for(auto* pa:to_delete){ - //remove this particle pointer from detParticles - particles.erase(find(particles.begin(), particles.end(), pa)); - } - } - //all done - } - - ///////////////// - //specific filters - //zvertex filter - bool doZVertexFilter(const clas12::region_particle *p) const{ - return _zVer->Filter(p->par()->getVz(),p->par()->getPid(),p->par()->getStatus(),_concurrent_key);// 0 => concurrent key - } - void doZVertexFilter() const{ - doFilters([this](const clas12::region_particle *p){return doZVertexFilter(p);}); - } - /* - //////////////// - //fiducial filter - bool doFiducialFilter(const clas12::region_particle *p) const - { - iguana::clas12::FiducialFilter::traj_row_data traj; - auto c12=C12(); - traj.x1 = p->traj(clas12::DC,6)->getX(); - traj.x1 = p->traj(clas12::DC,6)->getY(); - traj.z1 = p->traj(clas12::DC,6)->getZ(); - traj.x2 = p->traj(clas12::DC,18)->getX(); - traj.x2 = p->traj(clas12::DC,18)->getY(); - traj.z2 = p->traj(clas12::DC,18)->getZ(); - traj.x3 = p->traj(clas12::DC,36)->getX(); - traj.x3 = p->traj(clas12::DC,36)->getY(); - traj.z3 = p->traj(clas12::DC,36)->getZ(); - traj.sector = p->getSector(); - - //Need FiducialFitler to be public! - // return _Fid->Filter(traj,c12->runconfig()->getTorus(),p->par()->getPid()); - return true; - } - void doFiducialFilter() const{ - doFilters([this](const clas12::region_particle *p){return doFiducialFilter(p);}); - } - */ - ////////////////// - //PhotonGBT filter - /* - bool doPhotonGBTFilter(const clas12::region_particle *p, const std::map& calo) const{ - - if(p->getPid()!=22) return true; - - //Need FiducialFitler to be public! - - //return _PhotGBT->Filter(*(p->par()),*(p->cal()),calo,p->getIndex(),C12()->getRunNumber() ); - - return true; - - } - void doPhotonGBTFilter() const{ - if(getNParticles()==0) return; - auto calo = _PhotGBT->GetCaloMap(*(C12()->getDetParticles(0)->cal(-1))); - doFilters([this,&calo](const clas12::region_particle *p){return doPhotonGBTFilter(p,calo);}); - } - */ - /////////////////////////////////////////////// - void doAllFilters(){ - - if(_zVer.get())doZVertexFilter(); - // if(_Fid.get())doFiducialFilter(); - // if(_PhotGBT.get())doPhotonGBTFilter(); - } - - - //Reload - void Reload() const{ - auto run = C12()->runconfig()->getRun(); - if(_runNb==run) return; - _runNb=run; - - if(_zVer.get())_concurrent_key=_zVer->PrepareEvent(run); - - } - - private: - - std::unique_ptr _zVer; - //std::unique_ptr _Fid; - //std::unique_ptr _PhotGBT; - - mutable Long64_t _runNb =-1; - mutable int _concurrent_key=0; - - }; - -} - diff --git a/iguana/Iguana.h b/iguana/Iguana.h deleted file mode 100644 index 50ba7ea..0000000 --- a/iguana/Iguana.h +++ /dev/null @@ -1,107 +0,0 @@ -#pragma once - -#include "clas12reader.h" - -#include "Transformers.h" -#include "Filters.h" -#include "Creators.h" - - -namespace clas12root{ - - - class Iguana { - - public: - - void SetClas12( const std::unique_ptr* c12 ){ - _transform.SetClas12(c12); - _filter.SetClas12(c12); - _creator.SetClas12(c12); - } - - template - OPTION_TYPE SetOptionAll(std::string const& key, const OPTION_TYPE val){ - - _transform.SetOptionAll(key,val); - _filter.SetOptionAll(key,val); - _creator.SetOptionAll(key,val); - - return val; - } - - void Start(){ - - _transform.Start(); - _filter.Start(); - _creator.Start(); - - } - - void Stop(){ - - _transform.Stop(); - _filter.Stop(); - _creator.Stop(); - - } - - - Transformers& GetTransformers(){return _transform;} - Filters& GetFilters(){return _filter;} - Creators& GetCreators(){return _creator;} - - - void GoodToBad(const ROOT::Math::PxPyPzMVector& good, TLorentzVector &bad){ - bad.SetXYZT(good.X(),good.Y(),good.Z(),good.T()); - } - - //multithreading - concurrent_key_t PrepareEvent(int const runnum) const{ - /* - if(o_runnum->NeedsHashing()) { - std::hash hash_ftn; - auto hash_key = hash_ftn(runnum); - - if(!o_runnum->HasKey(hash_key)){ - _concurrent_key=hash_key; - // Reload(runnum, hash_key); - Reload(); - } - } else { - if(o_runnum->IsEmpty() || o_runnum->Load(0) != runnum){ - //Reload(runnum, 0); - _concurrent_key=0; - Reload(); - } - } - */ - Reload(); - return 0; - } - int GetKey()const {return _concurrent_key;} - - void Reload() const{ - //_transforms->Reload(GetKey()); - //_filter->Reload(GetKey()); - _filter.Reload(); - _creator.Reload(); - } - private: - - clas12root::Transformers _transform; - clas12root::Filters _filter; - clas12root::Creators _creator; - - int _concurrent_key=0; - // mutable std::unique_ptr> o_runnum = {ConcurrentParamFactory::Create()}; - - - }; - - - - - - -} diff --git a/iguana/IguanaAlgo.h b/iguana/IguanaAlgo.h deleted file mode 100644 index b12c2d4..0000000 --- a/iguana/IguanaAlgo.h +++ /dev/null @@ -1,68 +0,0 @@ -#pragma once - -#include "clas12reader.h" -#include "region_particle.h" - -#include -#include - -#include -#include -#include - -namespace clas12root{ - - using FourVector = ROOT::Math::PxPyPzMVector; - using concurrent_key_t = std::size_t; - - class IguanaAlgo{ - - public: - - void SetClas12( const std::unique_ptr* c12){_c12=c12;} - - void Start(){ - - for(auto* algo: _algos){ - std::cout<<" start "<GetName()<Start(); - } - - } - void Stop(){ - - for(auto* algo: _algos){ - algo->Stop(); - } - - } - - - protected : - - friend class Iguana; - - void AddAlgo(iguana::Algorithm* algo){ - _algos.push_back(algo); - } - - template - OPTION_TYPE SetOptionAll(std::string const& key, const OPTION_TYPE val){ - - for(auto* algo: _algos){ - algo->SetOption(key,val); - } - return val; - } - - clas12::clas12reader* C12() const {return _c12->get();} - - private: - - std::vector _algos; - const std::unique_ptr* _c12=nullptr; - - - }; - -} diff --git a/iguana/Transformers.h b/iguana/Transformers.h deleted file mode 100644 index 97425f2..0000000 --- a/iguana/Transformers.h +++ /dev/null @@ -1,137 +0,0 @@ -#pragma once - -#include "IguanaAlgo.h" -#include -#include -#include - -namespace clas12root{ - - - class Transformers : public IguanaAlgo { - - public: - - iguana::clas12::MomentumCorrection* MomentumCorrection(){return _momCorr.get();} - - void Use(const string& name){ - - //Initialise the momentum corrections algorithm - if(name=="clas12::MomentumCorrection"){ - _momCorr.reset(new iguana::clas12::MomentumCorrection); - AddAlgo(_momCorr.get()); - } - else if(name=="clas12::FTEnergyCorrection"){ - _ftECorr.reset(new iguana::clas12::FTEnergyCorrection); - AddAlgo(_ftECorr.get()); - } - else{ - std::cerr<<"Iguana -> clas12root::Transformers routine "<Transform(p->getPx(),p->getPy(),p->getPz(),p->getSector(),p->getPid(),C12()->runconfig()->getTorus()); - - //set the FourVector X,Y,Z,M - p4.SetCoordinates(px,py,pz,p4.M()); - - } - - //FT momentum corrections - void doFTEnergyCorrection(const clas12::region_particle *p, FourVector& p4) const{ - - auto [px, py, pz, en] = _ftECorr->Transform(p->getPx(),p->getPy(),p->getPz(),p->getSector()); - //set the FourVector X,Y,Z,T - p4.SetXYZT(px,py,pz,en); - - } - - - - //Configure to automate the right correction for the particle - void doAllCorrections(std::vector parts, - std::vector p4s) const{ - - if(parts.size()!=p4s.size()){ - std::cerr<<"Iguana-> clas12root::Transformers::doCorrections : Different number of particles, "<getRegion() ){ - //Forward tagger corrections - case clas12::FT : - if(_ftECorr.get()){ - doFTEnergyCorrection(parts[i],*p4s[i]); - } - - break; - - case clas12::FD : - //Forward detector corrections - if(_momCorr.get()){ - doMomentumCorrection(parts[i],*p4s[i]); - } - - break; - - case clas12::CD : - - break; - - default: - break; - }//end region switch - }//particle loop - - } - - - - - //interface to vector transforms - /*void doCorrections(std::vector parts, - std::vector p4s, - std::function< void (const clas12::region_particle *p, TLorentzVector& p4) > fun) const{ - - for(uint i=0;iTransform(p->getPx(),p->getPy(),p->getPz(),p->getSector(),p->getPid(),C12()->runconfig()->getTorus()); - - p4.SetXYZM(px,py,pz,p4.M()); - - } - - void doMomentumCorrections(std::vector parts, - std::vector p4s){ - - return doCorrections(parts,p4s, - [this](const clas12::region_particle *p,TLorentzVector& p4){ - return doMomentumCorrection(p,p4); - }); - } - */ - - private: - - std::unique_ptr _momCorr; - std::unique_ptr _ftECorr; - - - }; - -} - diff --git a/tests/Ex11_Iguana_MomentumCorrection.C b/tests/Ex11_Iguana_MomentumCorrection.C deleted file mode 100644 index c507dac..0000000 --- a/tests/Ex11_Iguana_MomentumCorrection.C +++ /dev/null @@ -1,127 +0,0 @@ -#include "clas12reader.h" -#include "region_particle.h" -#include "HipoChain.h" -#include "Iguana.h" -#include -#include -#include -#include -#include -#include -#include -#include - -void Ex11_Iguana_MomentumCorrection(){ - - using FourVector = ROOT::Math::PxPyPzMVector; - // clas12databases::SetRCDBRootConnection("rcdb.root"); - gBenchmark->Start("iguana"); - //create hipo chain for looping over events - clas12root::HipoChain chain; - chain.Add("/cache/clas12/rg-a/production/recon/spring2019/torus-1/pass2/dst/recon/006616/rec_clas_006616.evio.00000-00004.hipo"); - chain.SetReaderTags({0}); //create clas12reader with just tag 0 events - //add some filters - auto config_c12=chain.GetC12Reader(); - config_c12->addExactPid(11,1); //exactly 1 electron - config_c12->addExactPid(211,1); //exactly 1 pi+ - config_c12->addExactPid(-211,1); //exactly 1 pi- - config_c12->ignoreBank("REC::CovMat"); - - - // create the chosen algorithms - clas12root::Iguana ig; - ig.SetClas12(chain.C12ptr());//connect to clas12reader - ig.GetTransformers().Use("clas12::MomentumCorrection"); - ig.GetTransformers().Use("clas12::FTEnergyCorrection"); - ig.GetFilters().Use("clas12::ZVertexFilter"); - ig.GetCreators().Use("physics::InclusiveKinematics"); - - ig.SetOptionAll("log", "debug"); - ig.Start(); - - - //decalre 4-vector objects - auto db=TDatabasePDG::Instance(); - FourVector p4beam(0,0,10.6,10.6); - FourVector p4target(0,0,0,db->GetParticle(2212)->Mass()); - FourVector p4el(0,0,0,db->GetParticle(11)->Mass()); - FourVector p4pip(0,0,0,db->GetParticle(211)->Mass()); - FourVector p4pim(0,0,0,db->GetParticle(-211)->Mass()); - TLorentzVector tp4el(0,0,0,db->GetParticle(11)->Mass()); - - TH1D hQ2{"Q2","Q^{2}",100,0,10}; - TH1D hx{"x","x",100,0,1}; - TH1D hy{"y","y",100,0,1}; - TH1D hW{"W","W",100,0,6}; - TH1D hcQ2{"cQ2","corrected Q^{2}",100,0,10}; - TH1D hcx{"cx","corrected x",100,0,1}; - TH1D hcy{"cy","corrected y",100,0,1}; - TH1D hcW{"cW","corrected W",100,0,6}; - - - //get run clas12reader - auto& c12=chain.C12ref(); - while ( chain.Next() ){ - - //auto ebeam = c12->rcdb()->current().beam_energy/1000; - //p4beam.SetXYZT(0,0,ebeam,ebeam); //approx. mass =0 - // get particles by type - // note we applied a filter to ensure size of all ==1 - auto electron=c12->getByID(11)[0]; - auto pip=c12->getByID(211)[0]; - auto pim=c12->getByID(-211)[0]; - //filter on z-vertices of the particles - // if( !(ig.GetFilters().doZVertexFilter({electron,pip,pim})) ) { - if( !(ig.GetFilters().doAllFilters({electron,pip,pim})) ) { - continue; - } - - //correct momentum and get 4-vectors - ig.GetTransformers().doAllCorrections({electron,pip,pim},{&p4el,&p4pip,&p4pim}); - - //calculate inclusive kinematics - auto kine = ig.GetCreators().doInclusiveKinematics(electron); - auto corrkine = ig.GetCreators().doInclusiveKinematics(p4el); - - hQ2.Fill(kine.Q2); - hx.Fill(kine.x); - hy.Fill(kine.y); - hW.Fill(kine.W); - - hcQ2.Fill(corrkine.Q2); - hcx.Fill(corrkine.x); - hcy.Fill(corrkine.y); - hcW.Fill(corrkine.W); - - //if you unfortunately prefer TLorentzVector, there is a transform function to help. - // ig.GoodToBad(p4el,tp4el); - // cout<<"check good to bad "<Stop("iguana"); - gBenchmark->Print("iguana"); - - auto* can = new TCanvas(); - can->Divide(2,2); - - can->cd(1); - hQ2.DrawCopy(); - hcQ2.SetLineColor(2); - hcQ2.DrawCopy("same"); - - can->cd(2); - hW.DrawCopy(); - hcW.SetLineColor(2); - hcW.DrawCopy("same"); - - can->cd(3); - hx.DrawCopy(); - hcx.SetLineColor(2); - hcx.DrawCopy("same"); - - can->cd(4); - hy.DrawCopy(); - hcy.SetLineColor(2); - hcy.DrawCopy("same"); - -}