From 83adb5eadbf7b652a027cca7596ec97233fabc99 Mon Sep 17 00:00:00 2001
From: wenxingfang <1473717798@qq.com>
Date: Tue, 10 Nov 2020 16:19:49 +0800
Subject: [PATCH] replace ROOT Tree by NtupleSvc in GaudiPandora

---
 Examples/options/LCIO_read_pan.py             |  22 ++--
 Examples/options/tut_detsim_pandora.py        |  11 +-
 .../GaudiPandora/include/PandoraPFAlg.h       |  38 +++++-
 .../GaudiPandora/src/CaloHitCreator.cpp       |   2 +-
 .../Pandora/GaudiPandora/src/PandoraPFAlg.cpp | 108 ++++++++++++------
 5 files changed, 127 insertions(+), 54 deletions(-)

diff --git a/Examples/options/LCIO_read_pan.py b/Examples/options/LCIO_read_pan.py
index c85a3cae..ab9b93c8 100644
--- a/Examples/options/LCIO_read_pan.py
+++ b/Examples/options/LCIO_read_pan.py
@@ -8,9 +8,7 @@ dsvc = K4DataSvc("EventDataSvc")
 from Configurables import LCIOInput
 read = LCIOInput("read")
 read.inputs = [
-#"/cefs/data/FullSim/CEPC240/CEPC_v4/higgs/smart_final_states/E240.Pffh_invi.e0.p0.whizard195//ffh_inv.e0.p0.00001_1000_sim.slcio"
-#"/junofs/users/wxfang/CEPC/CEPCOFF/doReco/reco_output/nnh_aa.e0.p0.00010_000000_rec.slcio"
-"/junofs/users/wxfang/MyGit/tmp/fork_update_pandora/CEPCSW/Digi_sim_0.slcio"
+"/cefs/higgs/wxfang/cepc/Pandora/CaloDigi/gamma/Digi_sim_1.slcio"
 ]
 read.collections = [
         "MCParticle:MCParticle",
@@ -47,7 +45,11 @@ read.collections = [
 ##############################################################################
 from Configurables import GearSvc
 gearSvc  = GearSvc("GearSvc")
-gearSvc.GearXMLFile = "../Detector/DetCEPCv4/compact/FullDetGear.xml"
+gearSvc.GearXMLFile = "Detector/DetCEPCv4/compact/FullDetGear.xml"
+##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='LCIO_Pan_ana.root' OPT='NEW' TYP='ROOT'"]
 ##############################################################################
 from Configurables import PandoraPFAlg
 
@@ -55,6 +57,7 @@ pandoralg = PandoraPFAlg("PandoraPFAlg")
 pandoralg.use_dd4hep_geo     = False
 pandoralg.use_dd4hep_decoder = False
 pandoralg.use_preshower      = False
+pandoralg.WriteAna           = True
 pandoralg.collections = [
         "MCParticle:MCParticle",
         "CalorimeterHit:ECALBarrel",
@@ -77,9 +80,8 @@ pandoralg.collections = [
 pandoralg.WriteClusterCollection               = "PandoraClusters"              
 pandoralg.WriteReconstructedParticleCollection = "PandoraPFOs" 
 pandoralg.WriteVertexCollection                = "PandoraPFANewStartVertices"               
-pandoralg.AnaOutput = "Pandora_Ana.root"
 
-pandoralg.PandoraSettingsDefault_xml = "../Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml"
+pandoralg.PandoraSettingsDefault_xml = "Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml"
 #### Do not chage the collection name, only add or delete ###############
 pandoralg.TrackCollections      =  ["MarlinTrkTracks"]
 pandoralg.ECalCaloHitCollections=  ["ECALBarrel", "ECALEndcap", "ECALOther"]
@@ -119,16 +121,16 @@ pandoralg.AbsorberIntLengthOther= 0.006
 # write PODIO file
 from Configurables import PodioOutput
 write = PodioOutput("write")
-write.filename = "test.root"
+write.filename = "pan_test.root"
 write.outputCommands = ["keep *"]
 
 # ApplicationMgr
 from Configurables import ApplicationMgr
 ApplicationMgr(
-        #TopAlg = [read, pandoralg, write],
-        TopAlg = [read, pandoralg],
+        TopAlg = [read, pandoralg, write],
         EvtSel = 'NONE',
-        EvtMax = 1,
+        EvtMax = 10,
         ExtSvc = [dsvc, gearSvc],
+        HistogramPersistency = "ROOT",
         OutputLevel=INFO
 )
diff --git a/Examples/options/tut_detsim_pandora.py b/Examples/options/tut_detsim_pandora.py
index 98fd9617..82a314d8 100644
--- a/Examples/options/tut_detsim_pandora.py
+++ b/Examples/options/tut_detsim_pandora.py
@@ -108,7 +108,7 @@ example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
 from Configurables import CaloDigiAlg
 example_CaloDigiAlg = CaloDigiAlg("CaloDigiAlg")
 example_CaloDigiAlg.Scale = 1
-example_CaloDigiAlg.SimCaloHitCollection = "SimCalorimeterCol"
+example_CaloDigiAlg.SimCaloHitCollection = "EcalBarrelCollection"
 example_CaloDigiAlg.CaloHitCollection    = "ECALBarrel"
 example_CaloDigiAlg.CaloAssociationCollection    = "RecoCaloAssociation_ECALBarrel"
 ##############################################################################
@@ -116,12 +116,17 @@ from Configurables import GearSvc
 gearSvc  = GearSvc("GearSvc")
 gearSvc.GearXMLFile = "Detector/DetCEPCv4/compact/FullDetGear.xml"
 ##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='detsim_Pan_ana.root' OPT='NEW' TYP='ROOT'"]
+##############################################################################
 from Configurables import PandoraPFAlg
 
 pandoralg = PandoraPFAlg("PandoraPFAlg")
 pandoralg.use_dd4hep_geo     = True
 pandoralg.use_dd4hep_decoder = True
 pandoralg.use_preshower      = False
+pandoralg.WriteAna           = True
 pandoralg.collections = [
         "MCParticle:MCParticle",
         "CalorimeterHit:ECALBarrel",
@@ -144,7 +149,6 @@ pandoralg.collections = [
 pandoralg.WriteClusterCollection               = "PandoraClusters"              
 pandoralg.WriteReconstructedParticleCollection = "PandoraPFOs" 
 pandoralg.WriteVertexCollection                = "PandoraPFANewStartVertices"               
-pandoralg.AnaOutput = "Ana.root"
 
 pandoralg.PandoraSettingsDefault_xml = "Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml"
 #### Do not chage the collection name, only add or remove ###############
@@ -192,9 +196,10 @@ write.outputCommands = ["keep *"]
 # ApplicationMgr
 from Configurables import ApplicationMgr
 ApplicationMgr(
-        TopAlg = [genalg, detsimalg, example_CaloDigiAlg, pandoralg],
+        TopAlg = [genalg, detsimalg, example_CaloDigiAlg, pandoralg, write],
         EvtSel = 'NONE',
         EvtMax = 10,
         ExtSvc = [rndmengine, dsvc, geosvc, gearSvc,detsimsvc],
+        HistogramPersistency = "ROOT",
         OutputLevel=INFO
 )
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h
index 3f029fe3..67eb4835 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h
@@ -39,10 +39,13 @@
 #include "PfoCreator.h"
 #include "TrackCreator.h"
 
-#include "TROOT.h"
-#include "TTree.h"
-#include "TFile.h"
+//#include "TROOT.h"
+//#include "TTree.h"
+//#include "TFile.h"
 
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/NTuple.h"
+#include "GaudiKernel/MsgStream.h"
 
 /* PandoraPFAlg ========== <br>
  * 
@@ -254,6 +257,28 @@ protected:
   unsigned int                    m_nRun;                         ///< The run number
   unsigned int                    m_nEvent;                       ///< The event number
   //### For Ana #################
+  NTuple::Tuple* m_tuple = nullptr ;
+  NTuple::Item<long>   m_n_mc;
+  NTuple::Item<long>   m_n_rec;
+  NTuple::Item<int>   m_hasConversion;
+  NTuple::Array<int  > m_pReco_PID;    
+  NTuple::Array<float> m_pReco_mass;
+  NTuple::Array<float> m_pReco_energy;
+  NTuple::Array<float> m_pReco_px;
+  NTuple::Array<float> m_pReco_py;
+  NTuple::Array<float> m_pReco_pz;
+  NTuple::Array<float> m_pReco_charge;
+  NTuple::Array<int>   m_mc_p_size;
+  NTuple::Array<int>   m_mc_pid   ;
+  NTuple::Array<float> m_mc_mass  ;
+  NTuple::Array<float> m_mc_px    ;
+  NTuple::Array<float> m_mc_py    ;
+  NTuple::Array<float> m_mc_pz    ;
+  NTuple::Array<float> m_mc_charge;
+
+
+
+  /*
   TFile* m_fout;
   TTree* m_tree;
   std::vector<int  > m_pReco_PID;    
@@ -263,7 +288,6 @@ protected:
   std::vector<float> m_pReco_py;
   std::vector<float> m_pReco_pz;
   std::vector<float> m_pReco_charge;
-
   std::vector<int>   m_mc_p_size;
   std::vector<int>   m_mc_pid   ;
   std::vector<float> m_mc_mass  ;
@@ -272,11 +296,13 @@ protected:
   std::vector<float> m_mc_pz    ;
   std::vector<float> m_mc_charge;
   int m_hasConversion;
-
+  Gaudi::Property< std::string >              m_AnaOutput{ this, "AnaOutput", "Pan_Ana.root" };
+  */
+  
+  Gaudi::Property<bool> m_WriteAna {this, "WriteAna", false,"if do ana"};
   Gaudi::Property<bool> m_use_dd4hep_geo{this, "use_dd4hep_geo", false,"choose if use geo info from dd4hep"};
   Gaudi::Property<bool> m_use_dd4hep_decoder {this, "use_dd4hep_decoder", true,"if use decoder from dd4hep"};
   Gaudi::Property<bool> m_use_preshower {this, "use_preshower", false,"if use preshower layer for calorimeter"};
-  Gaudi::Property< std::string >              m_AnaOutput{ this, "AnaOutput", "Pan_Ana.root" };
   //######################
   std::map< std::string, std::string > m_collections;
   Gaudi::Property<std::vector<std::string>> m_readCols{this, "collections", {}, "Places of collections to read"};
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
index 818f78bf..a97fff97 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
@@ -789,7 +789,7 @@ void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *c
     const int physicalLayer(std::min(static_cast<int>(caloHitParameters.m_layer.Get()), static_cast<int>(layers.size()-1)));
     caloHitParameters.m_cellSize0 = layers[physicalLayer].cellSize0/dd4hep::mm;
     caloHitParameters.m_cellSize1 = layers[physicalLayer].cellSize1/dd4hep::mm;
-    std::cout<<"DD m_cellSize0="<<caloHitParameters.m_cellSize0.Get()<<",m_cellSize1="<<caloHitParameters.m_cellSize1.Get()<<std::endl;
+    //std::cout<<"DD m_cellSize0="<<caloHitParameters.m_cellSize0.Get()<<",m_cellSize1="<<caloHitParameters.m_cellSize1.Get()<<std::endl;
     double thickness = (layers[physicalLayer].inner_thickness+layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
     double nRadLengths = layers[physicalLayer].inner_nRadiationLengths;
     double nIntLengths = layers[physicalLayer].inner_nInteractionLengths;
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
index 6d458c60..bde13536 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
@@ -84,16 +84,47 @@ StatusCode PandoraPFAlg::initialize()
 {
 
   std::cout<<"init PandoraPFAlg"<<std::endl;
+  if(m_WriteAna){
 
+      NTuplePtr nt( ntupleSvc(), "MyTuples/Pan_reco_evt" );
+      if ( nt ) m_tuple = nt;
+      else {
+          m_tuple = ntupleSvc()->book( "MyTuples/Pan_reco_evt", CLID_ColumnWiseTuple, "Pan_reco_evt" );
+          if ( m_tuple ) {
+            m_tuple->addItem( "N_mc" , m_n_mc , 0, 1000 ).ignore();
+            m_tuple->addItem( "N_rec", m_n_rec, 0, 1000 ).ignore();
+            m_tuple->addItem( "m_pReco_PID"   , m_n_rec, m_pReco_PID    ).ignore();
+            m_tuple->addItem( "m_pReco_mass"  , m_n_rec, m_pReco_mass   ).ignore();
+            m_tuple->addItem( "m_pReco_energy", m_n_rec, m_pReco_energy ).ignore();
+            m_tuple->addItem( "m_pReco_px"    , m_n_rec, m_pReco_px     ).ignore();
+            m_tuple->addItem( "m_pReco_py"    , m_n_rec, m_pReco_py     ).ignore();
+            m_tuple->addItem( "m_pReco_pz"    , m_n_rec, m_pReco_pz     ).ignore();
+            m_tuple->addItem( "m_pReco_charge", m_n_rec, m_pReco_charge ).ignore();
+            m_tuple->addItem( "m_mc_p_size", m_n_mc  ,m_mc_p_size ).ignore();
+            m_tuple->addItem( "m_mc_pid"   , m_n_mc  ,m_mc_pid    ).ignore();
+            m_tuple->addItem( "m_mc_mass"  , m_n_mc  ,m_mc_mass   ).ignore();
+            m_tuple->addItem( "m_mc_px"    , m_n_mc  ,m_mc_px     ).ignore();
+            m_tuple->addItem( "m_mc_py"    , m_n_mc  ,m_mc_py     ).ignore();
+            m_tuple->addItem( "m_mc_pz"    , m_n_mc  ,m_mc_pz     ).ignore();
+            m_tuple->addItem( "m_mc_charge", m_n_mc  ,m_mc_charge ).ignore();
+            m_tuple->addItem( "m_hasConversion", m_hasConversion  ).ignore();
+          } 
+          else { // did not manage to book the N tuple....
+            error() << "    Cannot book N-tuple:" << long( m_tuple ) << endmsg;
+            return StatusCode::FAILURE;
+          }
+      }
+  }
+  /*
   std::string s_output =m_AnaOutput; 
   m_fout = new TFile(s_output.c_str(),"RECREATE"); 
   m_tree = new TTree("evt","tree");
-  m_tree->Branch("m_pReco_PID"   , &m_pReco_PID);
-  m_tree->Branch("m_pReco_mass"  , &m_pReco_mass);
+  m_tree->Branch("m_pReco_PID"   , &m_pReco_PID   );
+  m_tree->Branch("m_pReco_mass"  , &m_pReco_mass  );
   m_tree->Branch("m_pReco_energy", &m_pReco_energy);
-  m_tree->Branch("m_pReco_px"    , &m_pReco_px);
-  m_tree->Branch("m_pReco_py"    , &m_pReco_py);
-  m_tree->Branch("m_pReco_pz"    , &m_pReco_pz);
+  m_tree->Branch("m_pReco_px"    , &m_pReco_px    );
+  m_tree->Branch("m_pReco_py"    , &m_pReco_py    );
+  m_tree->Branch("m_pReco_pz"    , &m_pReco_pz    );
   m_tree->Branch("m_pReco_charge", &m_pReco_charge);
 
   m_tree->Branch("m_mc_p_size", &m_mc_p_size);
@@ -104,7 +135,7 @@ StatusCode PandoraPFAlg::initialize()
   m_tree->Branch("m_mc_pz"    , &m_mc_pz    );
   m_tree->Branch("m_mc_charge", &m_mc_charge);
   m_tree->Branch("m_hasConversion", &m_hasConversion);
-
+  */
   for ( const auto& col : m_readCols ) {
       auto seperater = col.find(':');
       std::string colType = col.substr(0, seperater);
@@ -334,8 +365,9 @@ StatusCode PandoraPFAlg::execute()
         PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pPfoCreator->CreateParticleFlowObjects(*m_CollectionMaps, m_ClusterCollection_w, m_ReconstructedParticleCollection_w, m_VertexCollection_w));
         
         StatusCode sc0 = CreateMCRecoParticleAssociation();
-        StatusCode sc = Ana();
-
+        if(m_WriteAna){
+            StatusCode sc = Ana();
+        }
         PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::Reset(*m_pPandora));
         this->Reset();
     }
@@ -358,10 +390,10 @@ StatusCode PandoraPFAlg::execute()
 
 StatusCode PandoraPFAlg::finalize()
 {
-  info() << "Finalized. Processed " << _nEvt << " events " <<",saved tree with entries="<<m_tree->GetEntries()<< endmsg;
-  m_fout->cd();
-  m_tree->Write();
-  m_fout->Close();
+  info() << "Finalized. Processed " << _nEvt << " events " << endmsg;
+  //m_fout->cd();
+  //m_tree->Write();
+  //m_fout->Close();
   delete m_pPandora;
   delete m_pGeometryCreator;
   delete m_pCaloHitCreator;
@@ -395,7 +427,7 @@ void PandoraPFAlg::Reset()
     m_pCaloHitCreator->Reset();
     m_pTrackCreator->Reset();
     m_pMCParticleCreator->Reset();
-
+    /*
     std::vector<int>()  .swap(m_pReco_PID   );
     std::vector<float>().swap(m_pReco_mass);
     std::vector<float>().swap(m_pReco_energy);
@@ -412,7 +444,7 @@ void PandoraPFAlg::Reset()
     std::vector<float>().swap(m_mc_pz    );
     std::vector<float>().swap(m_mc_charge);
     m_hasConversion = 0;
-
+    */
     m_CollectionMaps->clear();
 }
 
@@ -541,7 +573,9 @@ StatusCode PandoraPFAlg::updateMap()
 
 StatusCode PandoraPFAlg::Ana()
 {
-    int n_current = m_tree->GetEntries()+1;
+    m_n_mc = 0;
+    m_n_rec = 0;
+    m_hasConversion = 0;
     const edm4hep::ReconstructedParticleCollection* reco_col = m_ReconstructedParticleCollection_w.get();
     const edm4hep::MCRecoParticleAssociationCollection* reco_associa_col = m_MCRecoParticleAssociation_w.get();
     for(int i=0; i<reco_col->size();i++)
@@ -554,13 +588,14 @@ StatusCode PandoraPFAlg::Ana()
         const float mass = pReco.getMass();
         const float charge = pReco.getCharge();
         const int type = pReco.getType();
-        m_pReco_PID.push_back(type);
-        m_pReco_mass.push_back(mass);
-        m_pReco_charge.push_back(charge);
-        m_pReco_energy.push_back(energy);
-        m_pReco_px.push_back(px);
-        m_pReco_py.push_back(py);
-        m_pReco_pz.push_back(pz);
+        m_pReco_PID   [m_n_rec]=type;
+        m_pReco_mass  [m_n_rec]=mass;
+        m_pReco_charge[m_n_rec]=charge;
+        m_pReco_energy[m_n_rec]=energy;
+        m_pReco_px    [m_n_rec]=px;
+        m_pReco_py    [m_n_rec]=py;
+        m_pReco_pz    [m_n_rec]=pz;
+        m_n_rec ++ ;
         for(int j=0; j < reco_associa_col->size(); j++)
         {
             if(reco_associa_col->at(j).getRec().id() != pReco.id() ) continue;
@@ -568,20 +603,19 @@ StatusCode PandoraPFAlg::Ana()
         }
     }
     const edm4hep::MCParticleCollection*     MCParticle = nullptr;
-    StatusCode sc = StatusCode::SUCCESS;
-    sc =  getCol(m_mcParCol_r  , MCParticle );
-    if (NULL != MCParticle   )  
+    StatusCode sc =  getCol(m_mcParCol_r  , MCParticle );
+    if (NULL != MCParticle   )
     { 
         for(unsigned int i=0 ; i< MCParticle->size(); i++)
         {
-            m_mc_p_size.push_back(MCParticle->at(i).parents_size());
-            m_mc_pid   .push_back(MCParticle->at(i).getPDG());
-            m_mc_mass  .push_back(MCParticle->at(i).getMass());
-            m_mc_px    .push_back(MCParticle->at(i).getMomentum()[0]);
-            m_mc_py    .push_back(MCParticle->at(i).getMomentum()[1]);
-            m_mc_pz    .push_back(MCParticle->at(i).getMomentum()[2]);
-            m_mc_charge.push_back(MCParticle->at(i).getCharge());
-            //if(MCParticle->at(i).parents_size()==0) std::cout<<"MYDBUG evt="<<n_current<<", mc i="<<i<<",px="<<MCParticle->at(i).getMomentum()[0]<<",py="<<MCParticle->at(i).getMomentum()[1]<<",pz="<<MCParticle->at(i).getMomentum()[2]<<std::endl;
+            m_mc_p_size[m_n_mc] = MCParticle->at(i).parents_size();
+            m_mc_pid   [m_n_mc] = MCParticle->at(i).getPDG();
+            m_mc_mass  [m_n_mc] = MCParticle->at(i).getMass();
+            m_mc_px    [m_n_mc] = MCParticle->at(i).getMomentum()[0];
+            m_mc_py    [m_n_mc] = MCParticle->at(i).getMomentum()[1];
+            m_mc_pz    [m_n_mc] = MCParticle->at(i).getMomentum()[2];
+            m_mc_charge[m_n_mc] = MCParticle->at(i).getCharge();
+            m_n_mc ++ ;
             if (MCParticle->at(i).getPDG() != 22) continue;
             int hasEm = 0;
             int hasEp = 0;
@@ -593,7 +627,13 @@ StatusCode PandoraPFAlg::Ana()
             if(hasEm && hasEp) m_hasConversion=1;
         }
     }
-    m_tree->Fill();
+    //m_tree->Fill();
+    StatusCode status = m_tuple->write();
+    if ( status.isFailure() ) {
+        error() << "    Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
+        return StatusCode::FAILURE;
+    }
+
     return StatusCode::SUCCESS;
 }
 
-- 
GitLab