scaler reader and fix HipoChain to run sequntially through files with…

… clas12reader

Clas12Banks/clas12reader.cpp

@@ -5,14 +5,38 @@
  */
 
 #include "clas12reader.h"
-#include "scaler_reader.h"
 
 namespace clas12 {
 
-  clas12reader::clas12reader(std::string filename,std::vector<long> tags){
+  clas12reader::clas12reader(std::string filename,std::vector<long> tags):_filename(filename){
     cout<<" clas12reader::clas12reader reading "<<filename.data()<<endl;
     _reader.setTags(tags);
-    _reader.open(filename.data()); //keep a pointer to the reader
+
+    if(_filename.empty()==false)initReader();
+  }
+  ///////////////////////////////////////////////////
+  ///copy constructor
+  ///opens a new reader
+  ///Can give alternative filename
+  clas12reader::clas12reader(const clas12reader &other,std::string filename,std::vector<long> tags):_filename(filename){
+    cout<<" clas12reader::clas12reader reading "<<filename.data()<<endl;
+
+    //if default filename take same file as original
+    if(_filename.empty())_filename=other._filename;
+    if(other._scalReader.get()) scalerReader();
+
+    _reader.setTags(tags);
+    initReader();
+
+     _givenPids=other._givenPids;
+    _pidSelect=other._pidSelect;
+    _pidSelectExact=other._pidSelectExact;
+    _zeroOfRestPid=other._zeroOfRestPid;
+    _useFTBased=other._useFTBased;;
+
+  }
+  void clas12reader::initReader(){
+   _reader.open(_filename.data()); //keep a pointer to the reader
 
     // hipo::dictionary  factory;
     _reader.readDictionary(_factory);
@@ -44,13 +68,11 @@ namespace clas12 {
       _bcher.reset(new cherenkov{_factory.getSchema("REC::Cherenkov")});
     if(_factory.hasSchema("REC::ForwardTagger"))
       _bft.reset(new forwardtagger{_factory.getSchema("REC::ForwardTagger")});
- 
+
     makeListBanks();
 
-    scaler_reader sr(filename);
-    _runBeamCharge = sr.getBeamCharge();
-
   }
+
   ///////////////////////////////////////////////////////
   ///read the data
   void clas12reader::clearEvent(){

Clas12Banks/clas12reader.h

@@ -36,6 +36,8 @@
 #include "region_fdet.h"
 #include "region_cdet.h"
 #include "region_ft.h"
+#include "scaler_reader.h"
+
 #include "dictionary.h"
 
 #include <algorithm>
@@ -56,6 +58,8 @@ namespace clas12 {
 
     clas12reader()=default;
     clas12reader(std::string filename,std::vector<long> tags=std::vector<long>());
+    clas12reader(const clas12reader &other,std::string filename="",std::vector<long> tags=std::vector<long>());
+
     virtual ~clas12reader()=default;
 
     hipo::reader& getReader(){return _reader;}
@@ -146,15 +150,28 @@ namespace clas12 {
       return ( pattern & (1<<k)) != 0;
     }
 
+    scalerreader_ptr scalerReader(){
+      if(_scalReader.get()) return _scalReader.get();
+      _scalReader.reset(new scaler_reader(_filename));
+      _runBeamCharge = _scalReader->getBeamCharge();
+      return _scalReader.get();
+    }
     double getRunBeamCharge() const noexcept{ return _runBeamCharge;}
     double getCurrApproxCharge(){return _runBeamCharge*_nevent/_reader.getEntries();}
+
+    protected:
+
+    void initReader();
 
-  private:
 
+    private:
+
     void hipoRead(){
       _reader.read(_event);
       _isRead=true;
     }
+
+    std::string _filename;
 
     //reader
     hipo::reader     _reader;
@@ -180,6 +197,8 @@ namespace clas12 {
     ft_uptr _bft;
     vtp_uptr _bvtp;
 
+
+
     std::vector<std::unique_ptr<hipo::bank> > _addBanks; //owns additional banks
     std::vector<hipo::bank* > _allBanks; 
 
@@ -193,21 +212,25 @@ namespace clas12 {
     //this vector links to raw ptrs, does not own
     std::vector<region_part_ptr> _detParticles;
 
-    std::vector<short> _pids;
-    std::vector<short> _givenPids;
-    std::map<short,short> _pidSelect;
-    std::map<short,short> _pidSelectExact;
 
     double _runBeamCharge{0};
     long _nevent{0};
     ushort _nparts{0};
     ushort _n_rfdets{0};
     ushort _n_rcdets{0};
     ushort _n_rfts{0};
+
+    std::vector<short> _pids;
+    bool _isRead{false};
+
+    //members that need copied in constructor
+    scalerreader_uptr _scalReader;
+    std::vector<short> _givenPids;
+    std::map<short,short> _pidSelect;
+    std::map<short,short> _pidSelectExact;
     bool _zeroOfRestPid{false};
     bool _useFTBased{false};
-    bool _isRead{false};
-  };
+   };
   //helper functions
 
   //filter vectors via a lambda function

Clas12Banks/scaler_reader.cpp

@@ -1,24 +1,81 @@
 #include "scaler_reader.h"
+#include <algorithm>
 
 namespace clas12{
 
   scaler_reader::scaler_reader(std::string filename){
-    _reader.setTags(1);
-    _reader.open(filename.data()); //keep a pointer to the reader
-    _reader.readDictionary(_factory);
-    _factory.show();
-    _bscal.reset(new scaler{_factory.getSchema("RUN::scaler")});
-    _bscal->show();
-    while(_reader.next()){
-      _reader.read(_event);
-      _event.getStructure(*_bscal.get());
+    if(filename.empty())return;
+    hipo::reader     reader_;
+    hipo::event      event_;
+    hipo::dictionary  factory_;
+
+    reader_.setTags(1);
+    reader_.open(filename.data()); //keep a pointer to the reader
+    reader_.readDictionary(factory_);
+
+    _bscal.reset(new scaler{factory_.getSchema("RUN::scaler")});
+    _brun.reset(new runconfig{factory_.getSchema("RUN::config")});
+
+    while(reader_.next()){
+      reader_.read(event_);
+      event_.getStructure(*_bscal.get());
       if(_bscal->getRows()==0)continue; //only RUN::scaler
+      event_.getStructure(*_brun.get());
 
+      _nScalReads++;
       auto cup=_bscal->getFCupGated();
       if(cup>_maxCup) _maxCup=cup;
       if(cup<_minCup) _minCup=cup;
+        _beamCharges.push_back(std::make_pair(_brun->getEvent(),cup));
+   }
+    std::sort(_beamCharges.begin(),_beamCharges.end());
+
+    // int ccc=0;
+    _deltaCharges.resize(_nScalReads);
+    for(int i=0;i<_nScalReads;i++){
+      if(i)_deltaCharges[i]=_beamCharges[i].second-_beamCharges[i-1].second;
+      else _deltaCharges[i]=0;
     }
+    std::cout<<"Number of scaler reads "<<_nScalReads<<" "<<getBeamCharge()<<std::endl;
+  }
+
+  size_t  scaler_reader::addLongCounter(){
+    auto index=_longCounter.size();
+    std::vector<int64_t> vec(_nScalReads);
+    _longCounter.push_back(std::move(vec));
+    return index;
+  }
+  size_t  scaler_reader::addDoubleCounter(){
+    auto index=_doubleCounter.size();
+    std::vector<double> vec(_nScalReads);
+   _doubleCounter.push_back(std::move(vec));
+    return index;
   }
 
+  size_t scaler_reader::findPosition(size_t event){
 
+    if(event==_lastEvent) return _lastPosition;
+    _lastEvent=event;//new event
+
+    if(_lastPosition>0){
+      //beamCharges is sorted in event number
+      if(event> _beamCharges[_lastPosition-1].first
+	 && event<_beamCharges[_lastPosition].first)
+	return _lastPosition; //still in the same scaler read
+    }
+
+
+    auto itr=std::upper_bound(_beamCharges.begin(),_beamCharges.end(),
+			      evch_pair(event,0),
+			      [](const evch_pair& entry,const evch_pair& rhs){
+				return entry.first<rhs.first;
+			      }
+			      );
+    _lastPosition = std::distance(_beamCharges.begin(),itr);
+
+    //For standard events this should not happen, may happen with tag!=0
+    if(_lastPosition>=_nScalReads) _lastPosition=_nScalReads-1;
+    return  _lastPosition;
+  }
+
 }

Clas12Banks/scaler_reader.h

@@ -4,28 +4,51 @@
 #include "dictionary.h"
 #include "event.h"
 #include "scaler.h"
+#include "runconfig.h"
 #include <string>
 #include <float.h>
 
 namespace clas12 {
 
+  using evch_pair = std::pair<int,float>;
+
   class scaler_reader{
 
   public:
     scaler_reader(std::string filename);
 
     double getBeamCharge() const noexcept{return _maxCup-_minCup;}
 
+    size_t  addLongCounter();
+    size_t  addDoubleCounter();
+    size_t findPosition(size_t event);
+
+    void incrementLong(size_t index,size_t event,int64_t val){_longCounter[index][findPosition(event)]+=val;}
+    void incrementDouble(size_t index,size_t event,double val){_doubleCounter[index][findPosition(event)]+=val;}
+
+    std::vector<int64_t>&  getLongCounter(size_t index){return _longCounter[index];}
+    std::vector<double>&  getDoubleCounter(size_t index){return _doubleCounter[index];}
+    std::vector<float>&  getDeltaCharges(){return _deltaCharges;}
+    bool validCharge(size_t ev){ return validChargePos(findPosition(ev));}
+    //requres each accumated scaler>0 and delta charge>0
+    bool validChargePos(size_t pos){return pos>0 ?_beamCharges[pos].second>0&&_beamCharges[pos-1].second>0&&_deltaCharges[pos]>0 : false;}
+
   private:
-    //reader
-    hipo::reader     _reader;
-    hipo::event      _event;
-    hipo::dictionary  _factory;
-
+
     scaler_uptr _bscal;
+    runconfig_uptr _brun;
+
     double _maxCup{0};
     double _minCup{DBL_MAX};
-
+    std::vector<std::vector<double>> _doubleCounter{0};
+    std::vector<evch_pair> _beamCharges{0};
+    std::vector<std::vector<int64_t>> _longCounter{0};
+    std::vector<float>_deltaCharges{0};
+    long _nScalReads{0};
+    size_t _lastPosition{0};
+    size_t _lastEvent{0};
   };
+  using scalerreader_ptr=clas12::scaler_reader*;
+  using scalerreader_uptr=std::unique_ptr<clas12::scaler_reader>;
 
 }

Clas12Root/HipoChain.cpp

@@ -50,4 +50,31 @@ namespace clas12root {
     }
     return n-1;
   }
+
+  Bool_t HipoChain::Next(){
+    auto more = _c12->next();
+    if(more) return kTRUE;
+    // if(_c12.get()==nullptr)
+    if(_idxFile>0)_totBeamCharge+=_c12->getRunBeamCharge();
+    if(NextFile()==false) return kFALSE;
+    more = _c12->next();
+    if(more) return kTRUE;
+
+    //no more delete reader
+    _c12.reset();
+    //check if another file
+    return Next();
+
+  }
+
+  Bool_t HipoChain::NextFile(){
+    std::cout<<"HipoChain::NextFile() "<<_idxFile<<" out of "<<GetNFiles()<<std::endl;
+   if(_idxFile>=GetNFiles())
+      return kFALSE;//no more files
+    std::cout<<"HipoChain::NextFile() "<<GetFileName(_idxFile)<<std::endl;
+    //open next file
+   _c12.reset(new clas12::clas12reader{*_c12.get(),GetFileName(_idxFile++).Data(),_readerTags});
+    _c12ptr = _c12.get();	
+    return kTRUE;
+  }
 }

Clas12Root/HipoChain.h

@@ -1,6 +1,8 @@
 #pragma once
 
 
+#include "clas12reader.h"
+
 #include <TNamed.h>
 #include <TChain.h>
 #include <TObjArray.h>
@@ -18,7 +20,7 @@ namespace clas12root {
 
 
     void Add(TString name);
-    Int_t GetNFiles(){return _ListOfFiles->GetEntries();}
+    Int_t GetNFiles() const {return _ListOfFiles->GetEntries();}
     Long64_t GetNRecords();
 
     TString GetFileName(Int_t i){
@@ -32,13 +34,38 @@ namespace clas12root {
     };
     Long64_t GetRecordsToHere(Int_t ifile);
     Int_t GetFileFromRecord(Long64_t irec);
+
+    //Functions for sequential processing of chains of files
+    Bool_t Next();
+    Bool_t NextFile();
+    void SetReaderTags(std::vector<long> tags){_readerTags=tags;}
+
+    clas12::clas12reader* GetC12Reader() {
+      if( (_c12ptr=_c12.get()) )
+	return _c12ptr;
+      _c12.reset(new clas12::clas12reader{""});
+      _c12ptr = _c12.get();
+      return  _c12ptr;
+    }
+
+    void AddBeamCharge(Double_t bc){_totBeamCharge+=bc;}
+    Double_t TotalBeamCharge() const noexcept{return _totBeamCharge;}
+
   private :
     TChain _tchain;
 
-    TObjArray* _ListOfFiles=nullptr; //owned by _tchain
+    TObjArray* _ListOfFiles{nullptr}; //owned by _tchain
 
-    Long64_t _Nrecords=-1;
+    std::unique_ptr<clas12::clas12reader> _c12;
+    clas12::clas12reader* _c12ptr{nullptr};
+
+    Long64_t _Nrecords{-1};
     std::vector<Int_t> _fileRecords;
+
+    std::vector<long> _readerTags;
+    Int_t _idxFile{0};
+
+    Double_t _totBeamCharge{0};
 
     ClassDef(clas12root::HipoChain,1);
   };