From 191a3c00ee5fbc3512fc7b91a217eee3c66a506a Mon Sep 17 00:00:00 2001
From: wenxingfang <1473717798@qq.com>
Date: Fri, 6 Nov 2020 15:18:30 +0800
Subject: [PATCH] using NtupleSvc

---
 Digitisers/DCHDigi/src/DCHDigiAlg.cpp   | 116 ++++++++++++------------
 Digitisers/DCHDigi/src/DCHDigiAlg.h     |  39 ++++----
 Examples/options/tut_detsim_digi_SDT.py |  14 ++-
 3 files changed, 85 insertions(+), 84 deletions(-)

diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
index e20b6d10..d7c2949b 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
@@ -10,6 +10,9 @@
 #include <DD4hep/Objects.h>
 #include "DDRec/Vector3D.h"
 
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+
 
 #include <array>
 #include <math.h>
@@ -47,53 +50,44 @@ StatusCode DCHDigiAlg::initialize()
       return StatusCode::FAILURE;
   }
 
-  std::string s_output=m_Output; 
   if(m_WriteAna){
-      m_fout = new TFile(s_output.c_str(),"RECREATE"); 
-      m_tree = new TTree("evt","tree");
-      m_tree->Branch("chamber"  , &m_chamber  );
-      m_tree->Branch("layer"    , &m_layer    );
-      m_tree->Branch("cell"     , &m_cell     );
-      m_tree->Branch("cell_x"   , &m_cell_x   );
-      m_tree->Branch("cell_y"   , &m_cell_y   );
-      m_tree->Branch("simhit_x" , &m_simhit_x );
-      m_tree->Branch("simhit_y" , &m_simhit_y );
-      m_tree->Branch("simhit_z" , &m_simhit_z );
-      m_tree->Branch("hit_x"    , &m_hit_x    );
-      m_tree->Branch("hit_y"    , &m_hit_y    );
-      m_tree->Branch("hit_z"    , &m_hit_z    );
-      m_tree->Branch("dca"      , &m_dca      );
-      m_tree->Branch("hit_dE"   , &m_hit_dE   );
-      m_tree->Branch("hit_dE_dx", &m_hit_dE_dx);
+
+      NTuplePtr nt( ntupleSvc(), "MyTuples/DCH_digi_evt" );
+      if ( nt ) m_tuple = nt;
+      else {
+          m_tuple = ntupleSvc()->book( "MyTuples/DCH_digi_evt", CLID_ColumnWiseTuple, "DCH_digi_evt" );
+          if ( m_tuple ) {
+            m_tuple->addItem( "N_sim" , m_n_sim , 0, 100000 ).ignore();
+            m_tuple->addItem( "N_digi", m_n_digi, 0, 100000 ).ignore();
+            m_tuple->addItem( "simhit_x", m_n_sim, m_simhit_x).ignore();
+            m_tuple->addItem( "simhit_y", m_n_sim, m_simhit_y).ignore();
+            m_tuple->addItem( "simhit_z", m_n_sim, m_simhit_z).ignore();
+            m_tuple->addItem( "chamber" , m_n_digi, m_chamber  ).ignore();
+            m_tuple->addItem( "layer"   , m_n_digi, m_layer    ).ignore();
+            m_tuple->addItem( "cell"    , m_n_digi, m_cell     ).ignore();
+            m_tuple->addItem( "cell_x"  , m_n_digi, m_cell_x   ).ignore();
+            m_tuple->addItem( "cell_y"  , m_n_digi, m_cell_y   ).ignore();
+            m_tuple->addItem( "hit_x"    , m_n_digi,m_hit_x     ).ignore();
+            m_tuple->addItem( "hit_y"    , m_n_digi,m_hit_y     ).ignore();
+            m_tuple->addItem( "hit_z"    , m_n_digi,m_hit_z     ).ignore();
+            m_tuple->addItem( "dca"      , m_n_digi,m_dca       ).ignore();
+            m_tuple->addItem( "hit_dE"   , m_n_digi,m_hit_dE    ).ignore();
+            m_tuple->addItem( "hit_dE_dx", m_n_digi,m_hit_dE_dx ).ignore();
+          } 
+          else { // did not manage to book the N tuple....
+            error() << "    Cannot book N-tuple:" << long( m_tuple ) << endmsg;
+            return StatusCode::FAILURE;
+          }
+      }
   }
   std::cout<<"DCHDigiAlg::initialized"<< std::endl;
   return GaudiAlgorithm::initialize();
 }
 
-void DCHDigiAlg::Reset()
-{
-
-  std::vector<int  >().swap( m_chamber   );
-  std::vector<int  >().swap( m_layer     );
-  std::vector<int  >().swap( m_cell      );
-  std::vector<float>().swap( m_cell_x    );
-  std::vector<float>().swap( m_cell_y    );
-  std::vector<float>().swap( m_simhit_x  );
-  std::vector<float>().swap( m_simhit_y  );
-  std::vector<float>().swap( m_simhit_z  );
-  std::vector<float>().swap( m_hit_x     );
-  std::vector<float>().swap( m_hit_y     );
-  std::vector<float>().swap( m_hit_z     );
-  std::vector<float>().swap( m_dca       );
-  std::vector<float>().swap( m_hit_dE    );
-  std::vector<float>().swap( m_hit_dE_dx );
-
-}
 
 StatusCode DCHDigiAlg::execute()
 {
   info() << "Processing " << _nEvt << " events " << endmsg;
-  Reset();
   std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> > id_hits_map;
   edm4hep::TrackerHitCollection* Vec   = w_DigiDCHCol.createAndPut();
   edm4hep::MCRecoTrackerAssociationCollection* AssoVec   = w_AssociationCol.createAndPut();
@@ -116,6 +110,8 @@ StatusCode DCHDigiAlg::execute()
       }
   }
 
+  m_n_sim = 0;
+  m_n_digi = 0 ;
   for(std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> >::iterator iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++)
   {
     unsigned long long wcellid = iter->first;
@@ -165,9 +161,10 @@ StatusCode DCHDigiAlg::execute()
         asso.setWeight(iter->second.at(i).getEDep()/tot_edep);
 
         if(m_WriteAna){
-            m_simhit_x.push_back(pos.x());
-            m_simhit_y.push_back(pos.y());
-            m_simhit_z.push_back(pos.z());
+            m_simhit_x[m_n_sim] = pos.x();
+            m_simhit_y[m_n_sim] = pos.y();
+            m_simhit_z[m_n_sim] = pos.z();
+            m_n_sim ++ ;
         }
     }
     
@@ -178,34 +175,37 @@ StatusCode DCHDigiAlg::execute()
     trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm
 
     if(m_WriteAna){
-        m_chamber.push_back(chamber);
-        m_layer  .push_back(layer  );
-        m_cell      .push_back(cellID);
-        m_cell_x    .push_back(Wstart.x());
-        m_cell_y    .push_back(Wstart.y());
-        m_hit_x     .push_back(pos_x);
-        m_hit_y     .push_back(pos_y);
-        m_hit_z     .push_back(pos_z);
-        m_dca       .push_back(min_distance);
-        m_hit_dE    .push_back(trkHit.getEDep());
-        m_hit_dE_dx .push_back(trkHit.getEdx() );
+        m_chamber  [m_n_digi] = chamber;
+        m_layer    [m_n_digi] = layer  ;
+        m_cell     [m_n_digi] = cellID;
+        m_cell_x   [m_n_digi] = Wstart.x();
+        m_cell_y   [m_n_digi] = Wstart.y();
+        m_hit_x    [m_n_digi] = pos_x;
+        m_hit_y    [m_n_digi] = pos_y;
+        m_hit_z    [m_n_digi] = pos_z;
+        m_dca      [m_n_digi] = min_distance;
+        m_hit_dE   [m_n_digi] = trkHit.getEDep();
+        m_hit_dE_dx[m_n_digi] = trkHit.getEdx() ;
+        m_n_digi ++ ;
     }
   }
+
+
   std::cout<<"output digi DCHhit size="<< Vec->size() <<std::endl;
   _nEvt ++ ;
 
-  if(m_WriteAna) m_tree->Fill();
-
+  if(m_WriteAna){
+      StatusCode status = m_tuple->write();
+      if ( status.isFailure() ) {
+        error() << "    Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
+        return StatusCode::FAILURE;
+      }
+  }
   return StatusCode::SUCCESS;
 }
 
 StatusCode DCHDigiAlg::finalize()
 {
   info() << "Processed " << _nEvt << " events " << endmsg;
-  if(m_WriteAna){
-      m_fout->cd();
-      m_tree->Write();
-      m_fout->Close();
-  }
   return GaudiAlgorithm::finalize();
 }
diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h
index 7e38bf25..339f85ba 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.h
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h
@@ -2,6 +2,7 @@
 #define DCH_DIGI_ALG_H
 
 #include "FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
 #include "GaudiAlg/GaudiAlgorithm.h"
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/TrackerHitCollection.h"
@@ -14,9 +15,6 @@
 #include "DetSegmentation/GridDriftChamber.h"
 
 #include "TVector3.h"
-#include "TROOT.h"
-#include "TTree.h"
-#include "TFile.h"
 
 
 
@@ -39,7 +37,6 @@ public:
   /** Called after data processing for clean up.
    */
   virtual StatusCode finalize() ;
-  void Reset();
  
 protected:
 
@@ -47,22 +44,23 @@ protected:
   typedef std::vector<float> FloatVec;
   int _nEvt ;
 
-  TFile* m_fout;
-  TTree* m_tree;
-  std::vector<int  > m_chamber   ;
-  std::vector<int  > m_layer     ;
-  std::vector<int  > m_cell      ;
-  std::vector<float> m_cell_x    ;
-  std::vector<float> m_cell_y    ;
-  std::vector<float> m_simhit_x  ;
-  std::vector<float> m_simhit_y  ;
-  std::vector<float> m_simhit_z  ;
-  std::vector<float> m_hit_x     ;
-  std::vector<float> m_hit_y     ;
-  std::vector<float> m_hit_z     ;
-  std::vector<float> m_dca       ;
-  std::vector<float> m_hit_dE    ;
-  std::vector<float> m_hit_dE_dx ;
+  NTuple::Tuple* m_tuple = nullptr ;
+  NTuple::Item<long>   m_n_sim;
+  NTuple::Item<long>   m_n_digi;
+  NTuple::Array<int  > m_chamber   ;
+  NTuple::Array<int  > m_layer     ;
+  NTuple::Array<int  > m_cell      ;
+  NTuple::Array<float> m_cell_x    ;
+  NTuple::Array<float> m_cell_y    ;
+  NTuple::Array<float> m_simhit_x  ;
+  NTuple::Array<float> m_simhit_y  ;
+  NTuple::Array<float> m_simhit_z  ;
+  NTuple::Array<float> m_hit_x     ;
+  NTuple::Array<float> m_hit_y     ;
+  NTuple::Array<float> m_hit_z     ;
+  NTuple::Array<float> m_dca       ;
+  NTuple::Array<float> m_hit_dE    ;
+  NTuple::Array<float> m_hit_dE_dx ;
 
 
 
@@ -79,7 +77,6 @@ protected:
   Gaudi::Property<float> m_mom_threshold { this, "mom_threshold", 0};// GeV
   Gaudi::Property<bool>  m_WriteAna { this, "WriteAna", false};
 
-  Gaudi::Property<std::string> m_Output { this, "output", "ana_DCH_digi.root"};
 
   // Input collections
   DataHandle<edm4hep::SimTrackerHitCollection> r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};
diff --git a/Examples/options/tut_detsim_digi_SDT.py b/Examples/options/tut_detsim_digi_SDT.py
index a234b179..6511cf25 100644
--- a/Examples/options/tut_detsim_digi_SDT.py
+++ b/Examples/options/tut_detsim_digi_SDT.py
@@ -79,7 +79,7 @@ gun.ThetaMaxs = [90] # rad; 45deg
 gun.PhiMins = [90] # rad; 0deg
 gun.PhiMaxs = [90] # rad; 360deg
 #gun.PhiMins = [0] # rad; 0deg
-#gun.PhiMaxs = [0] # rad; 360deg
+#gun.PhiMaxs = [360] # rad; 360deg
 
 # stdheprdr = StdHepRdr("StdHepRdr")
 # stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
@@ -152,6 +152,10 @@ elif dedxoption == "BetheBlochEquationDedxSimTool":
     dedx_simtool.scale = 10
     dedx_simtool.resolution = 0.0001
 
+##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='DCH_digi_ana.root' OPT='NEW' TYP='ROOT'"]
 ##############################################################################
 from Configurables import DCHDigiAlg
 dCHDigiAlg = DCHDigiAlg("DCHDigiAlg")
@@ -162,14 +166,13 @@ dCHDigiAlg.SimDCHitCollection = "DriftChamberHitsCollection"
 dCHDigiAlg.DigiDCHitCollection = "DigiDCHitsCollection"
 dCHDigiAlg.AssociationCollection = "DCHAssociationCollectio"
 dCHDigiAlg.WriteAna  = True
-dCHDigiAlg.output  = "ana_DCH_digi_v1.root"
 
 ##############################################################################
 # POD I/O
 ##############################################################################
 from Configurables import PodioOutput
 out = PodioOutput("outputalg")
-out.filename = "detsim_digi_DCH.root"
+out.filename = "digi_DCH.root"
 out.outputCommands = ["keep *"]
 
 ##############################################################################
@@ -177,9 +180,10 @@ out.outputCommands = ["keep *"]
 ##############################################################################
 
 from Configurables import ApplicationMgr
-#ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, out],
-ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg],
+ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, out],
                 EvtSel = 'NONE',
                 EvtMax = 10,
                 ExtSvc = [rndmengine, dsvc, geosvc],
+                HistogramPersistency = "ROOT",
+                OutputLevel=INFO
 )
-- 
GitLab