Merge pull request #87 from dglazier/master

updates for iguana 0.8 and QADB 3.0

Clas12Banks/clas12databases.cpp

@@ -41,10 +41,12 @@ namespace clas12 {
     _myRcdbPath=other._myRcdbPath;
     _myCcdbPath=other._myCcdbPath;
     //_myQadbPath=other._myQadbPath;
-
+     _pass = other._pass;
+
     initDBs();
 
-   if(_qadb.get()){
+   if(other._qadb.get()){
+     initQA();
      qadb_requireOkForAsymmetry(other._qadbReqOKAsymmetry);
      qadb_requireGolden(other._qadbReqGolden);
      qadb_setQARequirements(other._qadbReqsQA);
@@ -57,11 +59,12 @@ namespace clas12 {
     _myRcdbPath=other._myRcdbPath;
     _myCcdbPath=other._myCcdbPath;
     //_myQadbPath=other._myQadbPath;
-
+    _pass = other._pass;
 
     initDBs();
 
-    if(_qadb.get()){
+    if(other._qadb.get()){
+      initQA();
       qadb_requireOkForAsymmetry(other._qadbReqOKAsymmetry);
       qadb_requireGolden(other._qadbReqGolden);
       qadb_setQARequirements(other._qadbReqsQA);
@@ -86,10 +89,10 @@ namespace clas12 {
     }
 
 
-    _qadb.reset( new qadb_reader{0} );
-    qadb_requireOkForAsymmetry(_qadbReqOKAsymmetry);
-    qadb_requireGolden(_qadbReqGolden);
-    qadb_setQARequirements(_qadbReqsQA);  
+    // _qadb.reset( new qadb_reader{0} );
+    // qadb_requireOkForAsymmetry(_qadbReqOKAsymmetry);
+    // qadb_requireGolden(_qadbReqGolden);
+    // qadb_setQARequirements(_qadbReqsQA);  
 
 }
   //update run number to each database
@@ -105,8 +108,11 @@ namespace clas12 {
      if(_ccdb.get()!=nullptr){
       _ccdb->setRun(runNb);
       _ccdb->updateTables();
-    }
-    if(_qadb.get()!=nullptr)_qadb->setRun(runNb);
+     }
+     if(_qadb.get()!=nullptr){
+       //     _qadb->setRun(runNb);
+       initQA();
+     }
     if(_rcdb.get()!=nullptr)_rcdb->readRun(runNb);
 
     close();

Clas12Banks/clas12databases.h

@@ -64,7 +64,12 @@ namespace clas12 {
     //qadb stuff
     //I would rather this was not needed here
     //but it is to make sure it gets passed to PROOF
-    void qadb_addQARequirement(string req){ _qadb->addQARequirement(req);_qadbReqsQA.push_back(req);};
+    void qadb_addQARequirement(string req){
+      if(_qadb.get()==nullptr){
+	throw std::runtime_error("clas12databases::qadb_addQARequirement - Need to call clas12reader::applyQA(cook) first");
+      }
+      _qadb->addQARequirement(req);_qadbReqsQA.push_back(req);
+    }
     void qadb_setQARequirements( std::vector<string> reqs){_qadb->setQARequirements(reqs); _qadbReqsQA = reqs;};
 
     void qadb_requireOkForAsymmetry(bool ok){_qadb->requireOkForAsymmetry(ok);_qadbReqOKAsymmetry=ok;};
@@ -75,13 +80,26 @@ namespace clas12 {
     void turnOffQADB(){_qadb.reset();}
     void turnOffRCDB(){_rcdb.reset();}
     void turnOffCCDB(){_ccdb.reset();}
+
+    void setPass(const string& pass){
+      _pass = pass;
+      initQA();
+    }
+    void initQA(){
+      std::cout<<"clas12databases::initQA pass = "<<_pass<<" run ="<<_runNb<<std::endl;
+     _qadb.reset( new qadb_reader{_pass,_runNb} );
+      qadb_requireOkForAsymmetry(_qadbReqOKAsymmetry);
+      qadb_requireGolden(_qadbReqGolden);
+      qadb_setQARequirements(_qadbReqsQA);  
+    }
   private:
 
    //names for copying to ROOT file for selector
     string _myRcdbPath;
     string  _myCcdbPath;
     //string  _myQadbPath;
 
+    string _pass; //!
     int _runNb={0};//!
 
     std::unique_ptr<rcdb_reader> _rcdb;//!

Clas12Banks/clas12reader.h

@@ -176,7 +176,7 @@ namespace clas12 {
     std::vector<region_part_ptr> getByID(int id);
     std::vector<region_part_ptr> getByRegion(int ir);
     std::vector<region_part_ptr> getByCharge(int ch);
-    
+
     const std::vector<short>& preCheckPids();
     const std::vector<short>& preCheckPidsOrCharge();
 
@@ -370,14 +370,9 @@ namespace clas12 {
     clas12databases* db(){return _db;};
 
     //clas12-qadb   
-    void applyQA() {
-      // if(_db)
-	if( _db->qa() )
-	  _applyQA=true;
-
-      if( _applyQA==false){
-	std::cout<<"Warning, clas12reader  applyQA() not valid"<<std::endl;
-      }
+    void applyQA(const string& pass) {
+      _applyQA=true;
+      _db->setPass(pass);
     }
 
     //double sumChargeFromQA();

Clas12Banks/qadb_reader.cpp

@@ -4,7 +4,7 @@ namespace clas12 {
 
 #ifdef CLAS_QADB
 
-  qadb_reader::qadb_reader(int runNb):_qa{}{_runNb=runNb;};
+  qadb_reader::qadb_reader(const string& pass,int runNb):_qa{pass}{_runNb=runNb;};
 
   ///////////////////////////////////////////////////////
   /// Add masks for all QA requirements
@@ -18,13 +18,13 @@ namespace clas12 {
   ///////////////////////////////////////////////////////
   ///Checks if an event passes all the QA requirements
   bool qadb_reader::passQAReqs(int evNb){
-    //  std::cout<<"DEBUG qadb_reader::passQAReqs "<<_runNb<<" "<<evNb<<std::endl;
+    // std::cout<<"DEBUG qadb_reader::passQAReqs "<<_runNb<<" "<<evNb<<std::endl;
     if(_runNb==0) return true;//e.g. simulation
 
     //isOkForAsymmetry already queries _QA, want to avoid doing it twice
     bool passAsymReq=true;
     bool queried=false;
-    // std::cout<<"DEBUG qadb_reader::passQAReqs "<<_runNb<<" "<<evNb<<"isgolden "<<isGolden(_runNb,evNb)<<" "<<isOkForAsymmetry(_runNb,evNb)<<std::endl;
+    // std::cout<<"DEBUG qadb_reader::passQAReqs "<<_runNb<<" "<<evNb<<"isgolden "<<_reqGolden<<" "<<isGolden(_runNb,evNb)<<" "<<isOkForAsymmetry(_runNb,evNb)<<std::endl;
     if(_reqOKAsymmetry){
       passAsymReq = isOkForAsymmetry(_runNb,evNb);
       queried=true;
@@ -38,7 +38,7 @@ namespace clas12 {
     if(passAsymReq && queried){
       //If an event is Golden it won't have other defects
       if(_reqGolden){
-	//std::cout<<"DEBUG qadb_reader::passQAReqs "<<_runNb<<" "<<evNb<<"isgolden "<<isGolden(_runNb,evNb)<<std::endl;
+	//	std::cout<<"DEBUG qadb_reader::passQAReqs "<<_runNb<<" "<<evNb<<"isgolden "<<isGolden(_runNb,evNb)<<std::endl;
 
 	//If the event passes the requirements, add charge
 	if(isGolden(_runNb,evNb)){

Clas12Banks/qadb_reader.h

@@ -16,7 +16,7 @@ namespace clas12 {
   class qadb_reader {
 
   public:
-    qadb_reader(int runNb=0);
+    qadb_reader(const string& pass,int runNb=0);
     //virtual ~qadb_reader()=default;  
 
 

RunRoot/Ex10_clas12DatabasesChain.C

@@ -42,9 +42,6 @@ void Ex10_clas12DatabasesChain(){
 
   //configure qadb requirements
    /*
-   * Several quality assurance requirements can be specified.
-   * requireOkForAsymmetry requires that an event was deemed
-   * suitable for asymmetry measurements.
    *
    * addRequirement requires that the event was not identified
    * as, for example, a marginal outlier. Several requirements
@@ -53,21 +50,14 @@ void Ex10_clas12DatabasesChain(){
    * See RGA analysis note and clas12-qadb github repository for
    * additional information.
    */
-  if(config_c12->qadb()!=nullptr){
-    config_c12->db()->qadb_requireOkForAsymmetry(true);
-    config_c12->db()->qadb_requireGolden(true);
-    config_c12->db()->qadb_addQARequirement("MarginalOutlier");
-    config_c12->db()->qadb_addQARequirement("TotalOutlier");
-   /*
-     * applyQA specifies to the clas12reader that quality assurance
-     * cuts will be applied, based on the .json file given as an 
-     * argument. This file should contain the Clas12 Quality Assurance
-     * database.
-     * i.e. clas12reader will only process events that pass quality assurance
-     * and ignore those which fail
-     */
-    config_c12->applyQA();
-  }
+  config_c12->applyQA(GETPASSSTRINGHERE);//GETPASSSTRINGHERE="latest", "pass1, "pass2",...
+  config_c12->db()->qadb_addQARequirement("MarginalOutlier");
+  config_c12->db()->qadb_addQARequirement("TotalOutlier");
+  config_c12->db()->qadb_addQARequirement("TerminalOutlier");
+  config_c12->db()->qadb_addQARequirement("MarginalOutlier");
+  config_c12->db()->qadb_addQARequirement("SectorLoss");
+  config_c12->db()->qadb_addQARequirement("LowLiveTime");
+
 
   gBenchmark->Start("db");
 

RunRoot/Ex11_Iguana_MomentumCorrection.C

@@ -22,9 +22,9 @@ void Ex11_Iguana_MomentumCorrection(){
   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->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");
 
 
@@ -62,20 +62,33 @@ void Ex11_Iguana_MomentumCorrection(){
   //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
-    // 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;
-    }
+    //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});
 
@@ -123,5 +136,6 @@ void Ex11_Iguana_MomentumCorrection(){
   hy.DrawCopy();
   hcy.SetLineColor(2);
   hcy.DrawCopy("same");
-
+
+  std::cout<<"analysed "<<hQ2.GetEntries()<<" events "<<std::endl;
 }

RunRoot/LoadIguana.C

@@ -7,6 +7,8 @@ void LoadIguana(){
   //external iguana linking
   TString IGUANA = gSystem->Getenv("IGUANA");
   gInterpreter->AddIncludePath(IGUANA+"/include");
+  //Some algorithms need extra paths
+  gInterpreter->AddIncludePath(IGUANA+"/../src/iguana/algorithms/clas12/PhotonGBTFilter/");
 
   gROOT->ProcessLine("#include <iguana/algorithms/AlgorithmBoilerplate.h>");
   gROOT->ProcessLine("#include <iguana/algorithms/AlgorithmSequence.h>");

iguana/Creators.h

@@ -25,20 +25,31 @@ namespace clas12root{
 
     //Inclusive Kinematics
     iguana::physics::InclusiveKinematicsVars doInclusiveKinematics(const clas12::region_particle *electron){
-      return _incKin->ComputeFromLepton(electron->getPx(),electron->getPy(),electron->getPz());
+      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());
+      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());
+      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<iguana::physics::InclusiveKinematics> _incKin;
-
+
+    mutable Long64_t _runNb =-1;
+    mutable int _concurrent_key=0;
+
   };
 
 }