From d00abf169b04988b5228386b3ed26b0c98cb3116 Mon Sep 17 00:00:00 2001
From: lintao <lintao51@gmail.com>
Date: Sun, 29 Nov 2020 11:40:16 +0800
Subject: [PATCH] WIP: add a python script and avoid the out of index.

---
 Analysis/TotalInvMass/src/TotalInvMass.cc     | 46 +++++++++++++------
 Examples/options/tut_analysis_TotalInvMass.py | 43 +++++++++++++++++
 2 files changed, 74 insertions(+), 15 deletions(-)
 create mode 100644 Examples/options/tut_analysis_TotalInvMass.py

diff --git a/Analysis/TotalInvMass/src/TotalInvMass.cc b/Analysis/TotalInvMass/src/TotalInvMass.cc
index 46891f1f..536739b2 100644
--- a/Analysis/TotalInvMass/src/TotalInvMass.cc
+++ b/Analysis/TotalInvMass/src/TotalInvMass.cc
@@ -59,7 +59,7 @@ TotalInvMass::TotalInvMass(const std::string& name, ISvcLocator* svcLoc)
 }
 
 StatusCode TotalInvMass::initialize() {
-
+    info() << "TotalInvMass::initializing..." << endmsg;
     // printParameters();
 
     TFile *tree_file=new TFile(_treeFileName.value().c_str(),(_overwrite ? "RECREATE" : "UPDATE"));
@@ -166,12 +166,15 @@ StatusCode TotalInvMass::initialize() {
 
     _Num = 0;
 
+    info() << "TotalInvMass::initializd" << endmsg;
+
      return GaudiAlgorithm::initialize();
 
 }
 
 StatusCode TotalInvMass::execute() 
 {		
+    info() << "TotalInvMass::executing..." << endmsg;
 
     EVENT::LCEvent* evtP = nullptr;
 
@@ -262,7 +265,7 @@ StatusCode TotalInvMass::execute()
     HcalHitColl.push_back("HCALBarrel");
     HcalHitColl.push_back("HCALEndcap");
     HcalHitColl.push_back("HCALOther");
-	
+
     try{
 
         for(int t = 0; t< int(hdl_EcalHitColl.size()); t++) {
@@ -501,31 +504,40 @@ StatusCode TotalInvMass::execute()
             P[2] = a_RecoP.getMomentum()[2];
 
             if(a_RecoP.clusters_size() > 0) {
+
                 auto currClu = a_RecoP.getClusters(0);
                 CluEn = currClu.getEnergy();
 
-                if(!Charge) {
+                if((!Charge) and (currClu.subdetectorEnergies_size() > 2)) {
                     CluEnCom[0] = currClu.getSubdetectorEnergies(0);
                     CluEnCom[1] = currClu.getSubdetectorEnergies(1);	
                     NeCaloE_a[0] += currClu.getSubdetectorEnergies(0);
                     NeCaloE_a[1] += currClu.getSubdetectorEnergies(1);
                 }
 
-                _EcalCluE += currClu.getSubdetectorEnergies(0);
-                _HcalCluE += currClu.getSubdetectorEnergies(1);
+                if (currClu.subdetectorEnergies_size() > 2) {
+                    _EcalCluE += currClu.getSubdetectorEnergies(0);
+                    _HcalCluE += currClu.getSubdetectorEnergies(1);
+                }
+
 
                 if(Charge) {
                     TrkSumEn += Energy; 
-                                
-                    auto a_Trk = a_RecoP.getTracks(0);
-                    TrackHit = a_Trk.trackerHits_size();
-                    StartPos[0] = (a_Trk.getTrackerHits(0)).getPosition()[0];
-                    StartPos[1] = (a_Trk.getTrackerHits(0)).getPosition()[1];
-                    StartPos[2] = (a_Trk.getTrackerHits(0)).getPosition()[2];
 
-                    EndPos[0] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[0];
-                    EndPos[1] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[1];
-                    EndPos[2] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[2];	
+                    if (a_RecoP.tracks_size()>0) {
+                        auto a_Trk = a_RecoP.getTracks(0);
+                        TrackHit = a_Trk.trackerHits_size();
+
+                        if (TrackHit>2) {
+                            StartPos[0] = (a_Trk.getTrackerHits(0)).getPosition()[0];
+                            StartPos[1] = (a_Trk.getTrackerHits(0)).getPosition()[1];
+                            StartPos[2] = (a_Trk.getTrackerHits(0)).getPosition()[2];
+
+                            EndPos[0] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[0];
+                            EndPos[1] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[1];
+                            EndPos[2] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[2];	
+                        }
+                    }
 
                     if( Energy > CluEn + sqrt(Energy)) {
                         ElargeP[0] += Energy; 
@@ -537,6 +549,7 @@ StatusCode TotalInvMass::execute()
                         EsmallP[0] += Energy;
                         EsmallP[1] += CluEn; 
                     }
+
                 }
             }
 
@@ -545,7 +558,6 @@ StatusCode TotalInvMass::execute()
         }
     }catch (lcio::DataNotAvailableException err) { }
 
-
     try{
         //LCCollection* col_RecoPandora = evtP->getCollection( "PandoraPFOs" );			
         //for(int i2 = 0; i2 < col_RecoPandora->getNumberOfElements(); i2++)
@@ -681,6 +693,10 @@ StatusCode TotalInvMass::execute()
     _outputTree->Fill();
     _Num++;
     // }  	  
+
+
+    info() << "TotalInvMass::execute done" << endmsg;
+
     return StatusCode::SUCCESS;
 
 }	
diff --git a/Examples/options/tut_analysis_TotalInvMass.py b/Examples/options/tut_analysis_TotalInvMass.py
new file mode 100644
index 00000000..8bbdc13c
--- /dev/null
+++ b/Examples/options/tut_analysis_TotalInvMass.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+
+from Gaudi.Configuration import *
+
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+# read LCIO files
+from Configurables import k4LCIOInput
+lcioinput = k4LCIOInput("k4LCIOInput")
+lcioinput.inputs = [
+"/cefs/higgs/yudan/CEPC240/Reco_tpc_1800/qqh/Reco_qqh__00001.slcio"
+]
+lcioinput.collections = [
+    "MCParticle:MCParticle",
+    "CalorimeterHit:ECALBarrel",
+    "CalorimeterHit:ECALEndcap",
+    "CalorimeterHit:HCALBarrel",
+    "CalorimeterHit:HCALEndcap",
+    "CalorimeterHit:HCALOther",
+    "ReconstructedParticle:AncientPFOs",
+    "ReconstructedParticle:ArborLICHPFOs"
+]
+
+from Configurables import TotalInvMass
+total_inv_mass = TotalInvMass("TotalInvMass")
+
+# # write PODIO file
+# from Configurables import PodioOutput
+# write = PodioOutput("write")
+# write.filename = "test.root"
+# write.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr(
+        TopAlg = [lcioinput, total_inv_mass],
+        EvtSel = 'NONE',
+        EvtMax = 10,
+        ExtSvc = [dsvc],
+        OutputLevel=INFO #DEBUG
+)
+
-- 
GitLab