From b517fff4924287103679df26795d827f06b3f578 Mon Sep 17 00:00:00 2001
From: Guido Volpi <Guido.Volpi@cern.ch>
Date: Wed, 18 Jun 2014 10:11:52 +0200
Subject: [PATCH] Clustering is allowed.
Clsutering header name changed. (FastTrackSimWrap-04-00-02)
---
.../FastTrackSimWrap/DumpSp.h | 190 ++
.../FastTrackSimWrap/FTKRegionalWrapper.h | 70 +
.../FastTrackSimWrap/UniqueBarcode.h | 55 +
.../TrigFTK/FastTrackSimWrap/cmt/requirements | 78 +
.../share/FastTrackSimNewWrap_jobOptions.py | 32 +
...kSimRegionalWrap_64TowersIBL_jobOptions.py | 45 +
...rackSimRegionalWrap_64Towers_jobOptions.py | 40 +
...astTrackSimRegionalWrap_8Reg_jobOptions.py | 37 +
.../share/FastTrackSimWrap_jobOptions.py | 45 +
.../FastTrackSimWrap/share/submit_torque.pl | 236 ++
.../TrigFTK/FastTrackSimWrap/src/DumpSp.cxx | 1896 +++++++++++++++++
.../src/FTKRegionalWrapper.cxx | 219 ++
.../components/FastTrackSimWrap_entries.cxx | 14 +
.../src/components/FastTrackSimWrap_load.cxx | 4 +
14 files changed, 2961 insertions(+)
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/DumpSp.h
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/FTKRegionalWrapper.h
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/UniqueBarcode.h
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/cmt/requirements
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimNewWrap_jobOptions.py
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64TowersIBL_jobOptions.py
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64Towers_jobOptions.py
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_8Reg_jobOptions.py
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimWrap_jobOptions.py
create mode 100755 Trigger/TrigFTK/FastTrackSimWrap/share/submit_torque.pl
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/src/DumpSp.cxx
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/src/FTKRegionalWrapper.cxx
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_entries.cxx
create mode 100644 Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_load.cxx
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/DumpSp.h b/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/DumpSp.h
new file mode 100644
index 00000000000..3da6e1944e4
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/DumpSp.h
@@ -0,0 +1,190 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef FASTTRACKSIMWRAP_DUMPSP
+#define FASTTRACKSIMWRAP_DUMPSP
+
+// FASTTRACKSIMWRAP
+// ================================================================
+// This algorithm is used to extract all needed information from
+// Athena for use in the standalone FTK simulation. Eventually an
+// Athena-based FTK simulation will exist and this code will
+// disappear.
+// ================================================================
+// 2009-05-20 Antonio Boveia (boveia@hep.uchicago.edu)
+// access SDOs and dump GEANT truth for each silicon
+// channel. perform matching of offline reconstructed
+// tracks to truth using this info. substantial rewrite
+// and cleanup.
+// 2007-12-10 Giulio Usai <giulio.usai@cern.ch>
+// Monica Dunford <mdunford@hep.uchicago.edu>
+// original version.
+
+
+
+#include <string>
+#include <map>
+#include <bitset>
+#include <fstream>
+#include <boost/iostreams/filtering_stream.hpp>
+#include <boost/shared_ptr.hpp>
+
+//#include "GaudiKernel/AthAlgorithm.h"
+#include "AthenaBaseComps/AthAlgorithm.h"
+//#include "GaudiKernel/StatusCode.h"
+//#include "GaudiKernel/ToolHandle.h"
+//#include "GaudiKernel/MsgStream.h"
+//#include "GaudiKernel/ServiceHandle.h"
+//#include "GaudiKernel/IJobOptionsSvc.h"
+
+#include "GeneratorObjects/HepMcParticleLink.h"
+#include "TrkToolInterfaces/ITruthToTrack.h"
+#include "TrigDecisionTool/TrigDecisionTool.h"
+#include "InDetConditionsSummaryService/IInDetConditionsSvc.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrkTruthData/TrackTruthCollection.h"
+#include "TrkTrackSummaryTool/TrackSummaryTool.h"
+#include "TrkToolInterfaces/ITrackHoleSearchTool.h"
+#include "InDetPrepRawData/SiClusterContainer.h"
+
+class AtlasDetectorID;
+class StoreGateSvc;
+class ITruthParameters;
+class TruthSelector;
+class PixelID;
+class SCT_ID;
+class TRT_ID;
+class IBeamCondSvc;
+class EventID;
+
+namespace InDetDD {
+ class SiDetectorManager;
+}
+namespace HepPDT {
+ class ParticleDataTable;
+}
+namespace HepMC {
+ class GenParticle;
+}
+
+class
+DumpSp : public AthAlgorithm
+{
+
+public:
+
+ typedef std::map<unsigned int,HepMcParticleLink::ExtendedBarCode> RecoToTruthTrackMap;
+ typedef std::map<HepMcParticleLink::ExtendedBarCode,unsigned int> TruthToRecoTrackMap;
+ typedef std::multimap<HepMcParticleLink::ExtendedBarCode,float> TruthToRecoProbMap;
+ typedef std::map<Identifier,int> HitIndexMap;
+
+public:
+
+ DumpSp( const std::string& name , ISvcLocator* pSvcLocator );
+ virtual ~DumpSp() {}
+
+ StatusCode initialize(void);
+ StatusCode execute(void);
+ StatusCode finalize(void);
+
+private:
+
+ RecoToTruthTrackMap _rttTrackMap; // mapping from reconstructed track to mc barcode
+ TruthToRecoTrackMap _ttrTrackMap; // mapping from mc barcode to reconstructed track (index into m_tracksName TrackCollection)
+ TruthToRecoProbMap _ttrProbMap; // mapping from mc barcode to matching fraction for corresponding rttTrackMap entry
+
+ typedef enum { TAU_PARENT_BIT , B_PARENT_BIT , PION_PARENT_BIT , PION_IMMEDIATE_PARENT_BIT , NBITS } Bits;
+ typedef std::bitset<NBITS> ParentBitmask;
+
+ void build_matching_maps(); // build geant truth to reconstruction maps for this event
+ const ParentBitmask construct_truth_bitmap( const HepMC::GenParticle* p ) const;
+ void dump_spacepoints() const; // dump clustered spacepoints to text file
+ void dump_raw_silicon( HitIndexMap& hitIndexMap, HitIndexMap& pixelClusterIndexMap ) const; // dump raw silicon data to text file and populate hitIndexMap for rec. track processing
+ void dump_tracks( const HitIndexMap& hitIndexMap, const HitIndexMap& pixelClusterIndexMap ) const; // dump reconstructed tracks in new truth matching format
+ void dump_truth() const; // dump old track-truth matching format to text file
+ void dump_beamspot() const; // dump beamspot to wrapper output
+ void dump_bad_modules() const; // dump list of bad modules to the output
+ void dump_MBTS( ) const; // dump MBTS data for data/MC occupancy studies
+ void dump_vertex( ) const; // dump offline vertex info
+ void dump_lumibcid( const EventID* evid ) const; // dump lumiblock Number and bcid
+
+private:
+ // tools
+ // ================================================================
+
+
+ boost::shared_ptr<AtlasDetectorID> m_idHelper;
+ StoreGateSvc* m_storeGate;
+ StoreGateSvc* m_detStore;
+ StoreGateSvc* m_evtStore;
+
+ const PixelID* m_pixelId;
+ const SCT_ID* m_sctId;
+ const TRT_ID* m_trtId;
+
+ const InDetDD::SiDetectorManager* m_PIX_mgr;
+ const InDetDD::SiDetectorManager* m_SCT_mgr;
+ const InDetDD::SiDetectorManager* m_TRT_mgr;
+
+ const InDet::SiClusterContainer* m_pixelContainer;
+ const InDet::SiClusterContainer* m_sctContainer;
+
+ const HepPDT::ParticleDataTable* m_particleDataTable;
+
+ ToolHandle<Trk::ITruthToTrack> m_truthToTrack; //!< tool to create track parameters from a gen particle
+ ToolHandle<Trk::ITrackSummaryTool> m_trkSummaryTool;
+
+ //#ifdef HAVE_VERSION_15
+ ToolHandle<Trig::TrigDecisionTool> m_trigDecTool; // tool to retrieve trigger decision
+ ServiceHandle<IInDetConditionsSvc> m_pixelCondSummarySvc; // tool to retrieve pixel conditions db
+ ServiceHandle<IInDetConditionsSvc> m_sctCondSummarySvc; // tool to retrieve SCT conditions db
+ //#else
+ // ToolHandle<Trk::TruthToTrack> m_truthToTrack; //!< tool to create track parameters from a gen particle
+ //#endif
+ // ToolHandle<Trk::ITrackHoleSearchTool> m_holeSearchTool;
+
+ ServiceHandle<IBeamCondSvc> m_beamCondSvc;
+
+ // job configuration
+
+ std::string m_pixelClustersName;
+ std::string m_sctClustersName;
+ std::string m_pixelSpacePointsName;
+ std::string m_sctSpacePointsName;
+ std::string m_overlapSpacePointsName;
+ std::string m_tracksName;
+ std::string m_tracksTruthName;
+ double m_maxEta;
+ double m_minPt;
+
+ std::string m_spacePointsName;
+ std::string m_outFileName;
+ std::string m_outFileNameRawHits;
+
+ bool m_doTrigger;
+ bool m_doData;
+ bool m_doVertex;
+ bool m_doLumiBCID;
+ bool m_doBadMod;
+
+ bool m_dumpHitsOnTracks;
+ bool m_dumpSpacePoints;
+ bool m_dumpTruthIntersections;
+ bool m_dumpMBTSInfo;
+ bool m_useOfflineTrackSelectorTool;
+ bool m_outputBeamSpotToWrapper;
+ bool m_useSimpleCuts;
+
+ bool m_doRDODebug;
+
+
+ // output
+ // ================================================================
+
+ mutable boost::shared_ptr<boost::iostreams::filtering_ostream> ofl;
+ mutable boost::shared_ptr<boost::iostreams::filtering_ostream> oflraw;
+
+};
+
+#endif
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/FTKRegionalWrapper.h b/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/FTKRegionalWrapper.h
new file mode 100644
index 00000000000..0f2ec11cef1
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/FTKRegionalWrapper.h
@@ -0,0 +1,70 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef FTKRegionalWrapper_h
+#define FTKRegionalWrapper_h
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+
+#include "TrigFTKSim/FTK_SGHitInput.h"
+#include "TrigFTKSim/FTKRegionMap.h"
+#include "TrigFTKSim/FTKPMap.h"
+
+#include "TFile.h"
+#include "TTree.h"
+
+#include <string>
+
+class FTKRegionalWrapper : public AthAlgorithm {
+public:
+ FTKRegionalWrapper (const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~FTKRegionalWrapper ();
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+
+private:
+ ToolHandle<FTK_SGHitInputI> m_hitInputTool; // input handler
+
+ // variables to manage the distribution of the hits
+ bool m_IBLMode; // global FTK setup variable to handle IBL
+ std::string m_pmap_path; // path of the PMAP file
+ FTKPlaneMap *m_pmap; // pointer to the pmap object
+
+ // variables to manage the region maps
+ std::string m_rmap_path; // path of the region-map file
+ FTKRegionMap *m_rmap; // pointer to the RMAP object
+ int m_ntowers;
+
+ bool m_SaveRawHits; // flag to allow to store FTKRawHit collections, pmap non applied
+ bool m_SaveHits; // flag to allow to store FTKHit collections, pmap applied
+
+ /* the following flags were copied from TrigFTKSim/FTKDataInput.h on 11/02/14 check they are still up to date. Thanks, Guido*/
+ bool m_Clustering; // flag to enable the clustering
+ bool m_SaveClusterContent; // flag to allow the debug of the cluster content
+ bool m_DiagClustering; // true if diagonal clustering is used
+ bool m_SctClustering; // true if SCT clustering is performed
+ unsigned int m_PixelClusteringMode; // 1 means ToT information is used
+ // && 400/600um pixels are accounted for
+ // 0 is default used up to 2013
+ bool m_DuplicateGanged;
+ bool m_GangedPatternRecognition;
+
+ // variables related to the output fiels
+ std::string m_outpath; // output path
+ TFile *m_outfile; // ROOT file descriptor
+
+ TTree *m_hittree; // TTree for the hit storage
+ std::vector<FTKRawHit> *m_original_hits; // variables related to the FTKRawHit storage
+ std::vector<FTKHit> *m_logical_hits; // variables related to the FTKHit storage
+
+ TTree *m_evtinfo; // TTree with general event information
+ int m_run_number; // event's run number
+ int m_event_number; // event number
+
+ TTree *m_trackstree;
+ std::vector<FTKTruthTrack> m_truth_tracks;
+};
+
+#endif // FTKRegionalWrapper_h
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/UniqueBarcode.h b/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/UniqueBarcode.h
new file mode 100644
index 00000000000..d24c9189a73
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/FastTrackSimWrap/UniqueBarcode.h
@@ -0,0 +1,55 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef WRAP_UNIQUEBARCODE_H
+#define WRAP_UNIQUEBARCODE_H
+
+#include <cmath>
+#include <boost/tuple/tuple.hpp>
+
+class
+UniqueBarcode
+{
+
+private:
+
+ mutable bool _valid;
+ int _event_index;
+ int _barcode;
+
+public:
+
+ UniqueBarcode() : _valid( false ) {}
+
+ UniqueBarcode( const int& event_index , const int& barcode )
+ : _valid(true)
+ , _event_index( event_index )
+ , _barcode( barcode )
+ {}
+
+ UniqueBarcode( const signed long& ucode )
+ : _valid(true)
+ {
+ const boost::tuple<int,int> tmp( barcode_and_event_index(ucode) );
+ _barcode = tmp.get<0>();
+ _event_index = tmp.get<1>();
+ }
+
+ virtual ~UniqueBarcode() {}
+
+ static const signed long code( const int& event_index , const int& barcode ) { return( (100000l*event_index) + barcode ); }
+ static const boost::tuple<int,int> barcode_and_event_index( const signed long& ucode ) {
+ ldiv_t tmp = std::div( ucode , 100000l );
+ return boost::tuple<int,int>(tmp.rem,tmp.quot);
+ }
+
+ const int event_index() const { return _event_index; }
+ const int barcode() const { return _barcode; }
+
+ const signed long operator()() const { return code(_event_index,_barcode); }
+ const bool operator==(const UniqueBarcode& rhs) const { return( _valid && rhs._valid && _event_index==rhs._event_index && _barcode==rhs._barcode ); }
+
+};
+
+#endif // WRAP_UNIQUEBARCODE_H
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/cmt/requirements b/Trigger/TrigFTK/FastTrackSimWrap/cmt/requirements
new file mode 100644
index 00000000000..346e2d034ea
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/cmt/requirements
@@ -0,0 +1,78 @@
+package FastTrackSimWrap
+
+# This algorithm is used to extract all needed information from
+# Athena for use in the standalone FTK simulation. Eventually an
+# Athena-based FTK simulation will exist and this code will
+# disappear.
+
+author Monica Dunford <mdunford@hep.uchicago.edu>
+author Antonio Boveia <boveia@hep.uchicago.edu>
+author Guido Volpi <guido.volpi@cern.ch>
+
+use AthenaBaseComps AthenaBaseComps-* Control
+use AtlasROOT AtlasROOT-* External
+use TrigFTKSim TrigFTKSim-* Trigger/TrigFTK
+use AtlasReconstructionRunTime AtlasReconstructionRunTime-*
+use AtlasPolicy AtlasPolicy-*
+use GeneratorObjects GeneratorObjects-* Generators
+#use TruthHelper TruthHelper-* Generators/GenAnalysisTools
+#use HepMC HepMC-* Simulation
+use InDetPrepRawData InDetPrepRawData-* InnerDetector/InDetRecEvent
+use TrigCaloEvent TrigCaloEvent-* Trigger/TrigEvent
+use TrigDecisionTool TrigDecisionTool-* Trigger/TrigAnalysis
+use TrkTrack TrkTrack-* Tracking/TrkEvent
+use TrkTruthData TrkTruthData-* Tracking/TrkEvent
+#use TrkTruthToTrack TrkTruthToTrack-* Tracking/TrkTools
+use TrkTrackSummaryTool TrkTrackSummaryTool-* Tracking/TrkTools
+use TrkToolInterfaces TrkToolInterfaces-* Tracking/TrkTools
+#use TruthTools TruthTools-* Generators/GenAnalysisTools
+#use Particle Particle-* Reconstruction
+#use iPatRecEvent iPatRecEvent-* Reconstruction/iPat
+#use iPatTruthTrajectory iPatTruthTrajectory-* Reconstruction/iPat
+#use iPatTrackParameters iPatTrackParameters-* Reconstruction/iPat
+#use iPatInterfaces iPatInterfaces-* Reconstruction/iPat
+#use iPatTrajectory iPatTrajectory-* Reconstruction/iPat
+use InDetConditionsSummaryService InDetConditionsSummaryService-* InnerDetector/InDetConditions
+#use Boost v* LCG_Interfaces
+use AtlasBoost AtlasBoost-* External
+#use dcache_client v* LCG_Interfaces
+
+private
+use GaudiInterface GaudiInterface-* External
+use AtlasDetDescr AtlasDetDescr-* DetectorDescription
+use AtlasHepMC AtlasHepMC-* External
+use EventInfo EventInfo-* Event
+use HepPDT v* LCG_Interfaces
+use IdDict IdDict-* DetectorDescription
+use InDetBeamSpotService InDetBeamSpotService-* InnerDetector/InDetConditions
+use InDetIdentifier InDetIdentifier-* InnerDetector/InDetDetDescr
+use InDetRawData InDetRawData-* InnerDetector/InDetRawEvent
+use InDetReadoutGeometry InDetReadoutGeometry-* InnerDetector/InDetDetDescr
+use InDetSimData InDetSimData-* InnerDetector/InDetRawEvent
+use IdDictDetDescr IdDictDetDescr-* DetectorDescription
+use Identifier Identifier-* DetectorDescription
+use StoreGate StoreGate-* Control
+use TileIdentifier TileIdentifier-* TileCalorimeter
+use TrkSpacePoint TrkSpacePoint-* Tracking/TrkEvent
+use EventPrimitives EventPrimitives-* Event
+#use TrkParameters TrkParameters-* Tracking/TrkEvent
+use TrkRIO_OnTrack TrkRIO_OnTrack-* Tracking/TrkEvent
+use VxVertex VxVertex-* Tracking/TrkEvent
+end_private
+
+apply_pattern dual_use_library files=*.cxx
+apply_pattern declare_joboptions files="*jobOptions*.py"
+
+# fix missing libboost_iostreams interface in LCGCMT in 13.0.40
+#macro Boost_linkopts_iostreams " -lboost_iostreams-gcc-mt-$(Boost_file_version) " \
+# Darwin " -lboost_iostreams-d-$(Boost_file_version) " \
+# WIN32 " boost_iostreams-vc71-mt-$(Boost_file_version).lib "
+
+# fix missing libboost_iostreams in LCGCMT in 15.2.0
+#macro Boost_linkopts_iostreams " -lboost_iostreams-$(Boost_compiler_version)-mt " \
+# target-darwin " -lboost_iostreams-$(Boost_compiler_version)-$(Boost_file_version) " \
+# target-winxp " boost_iostreams-$(Boost_compiler_version)-$(Boost_file_version).lib "
+
+private
+macro_append Boost_linkopts " $(Boost_linkopts_iostreams) "
+
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimNewWrap_jobOptions.py b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimNewWrap_jobOptions.py
new file mode 100644
index 00000000000..d09f23c5876
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimNewWrap_jobOptions.py
@@ -0,0 +1,32 @@
+###############################################################
+#
+# DumpSp job options file
+#
+#==============================================================
+
+# Helper function from transforms
+from PyJobTransforms.trfUtils import findFile
+pmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_11L.pmap')
+print "Using PMAP:", pmap_path
+rmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_11L.tmap')
+print "Using RMAP:", rmap_path
+
+from AthenaCommon.AlgSequence import AlgSequence
+theJob = AlgSequence()
+
+
+from FastTrackSimWrap.FastTrackSimWrapConf import FTKRegionalWrapper
+if hasattr(runArgs,"outputNTUP_FTKIPFile") :
+ OutputNTUP_FTKIPFile = runArgs.outputNTUP_FTKIPFile
+else :
+ OutputNTUP_FTKIPFile = 'ftksim_newwrap_raw.root'
+
+print "Output file", OutputNTUP_FTKIPFile
+
+theJob += FTKRegionalWrapper(OutputLevel = DEBUG,
+ PMapPath = pmap_path,
+ RMapPath = rmap_path,
+ OutFileName = "ftksim_regionalwrap.root")
+
+print theJob
+###############################################################
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64TowersIBL_jobOptions.py b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64TowersIBL_jobOptions.py
new file mode 100644
index 00000000000..eb4bf1a9876
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64TowersIBL_jobOptions.py
@@ -0,0 +1,45 @@
+###############################################################
+#
+# DumpSp job options file
+#
+#==============================================================
+
+# Helper function from transforms
+from PyJobTransforms.trfUtils import findFile
+pmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_12Libl.pmap')
+print "Using PMAP:", pmap_path
+rmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_12Libl.tmap')
+print "Using RMAP:", rmap_path
+
+
+from AthenaCommon.AlgSequence import AlgSequence
+theJob = AlgSequence()
+
+
+from FastTrackSimWrap.FastTrackSimWrapConf import FTKRegionalWrapper
+if hasattr(runArgs,"outputNTUP_FTKIPFile") :
+ OutputNTUP_FTKIPFile = runArgs.outputNTUP_FTKIPFile
+else :
+ OutputNTUP_FTKIPFile = "ftksim_64Towers_wrap.root"
+
+from AthenaCommon.AppMgr import ToolSvc
+
+from TrigFTKSim.TrigFTKSimConf import FTK_SGHitInput
+FTKSGInput = FTK_SGHitInput( maxEta= 3.2, minPt= 0.8*GeV)
+FTKSGInput.OutputLevel = DEBUG
+FTKSGInput.ReadTruthTracks = True
+
+ToolSvc += FTKSGInput
+
+print "Output file", OutputNTUP_FTKIPFile
+wrapper = FTKRegionalWrapper(OutputLevel = DEBUG,
+ PMapPath = pmap_path,
+ RMapPath = rmap_path,
+ OutFileName = OutputNTUP_FTKIPFile)
+wrapper.IBLMode = True
+wrapper.HitInputTool = FTKSGInput
+#wrapper.Clustering = True
+theJob += wrapper
+
+print theJob
+###############################################################
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64Towers_jobOptions.py b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64Towers_jobOptions.py
new file mode 100644
index 00000000000..6fa13b6225f
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_64Towers_jobOptions.py
@@ -0,0 +1,40 @@
+###############################################################
+#
+# DumpSp job options file
+#
+#==============================================================
+
+# Helper function from transforms
+from PyJobTransforms.trfUtils import findFile
+pmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_11L.pmap')
+print "Using PMAP:", pmap_path
+rmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_11L.tmap')
+print "Using RMAP:", rmap_path
+
+
+from AthenaCommon.AlgSequence import AlgSequence
+theJob = AlgSequence()
+
+
+from FastTrackSimWrap.FastTrackSimWrapConf import FTKRegionalWrapper
+if hasattr(runArgs,"outputNTUP_FTKIPFile") :
+ OutputNTUP_FTKIPFile = runArgs.outputNTUP_FTKIPFile
+else :
+ OutputNTUP_FTKIPFile = "ftksim_64Towers_wrap.root"
+
+from AthenaCommon.AppMgr import ToolSvc
+
+from TrigFTKSim.TrigFTKSimConf import FTK_SGHitInput
+FTKSGInput = FTK_SGHitInput()
+FTKSGInput.ReadTruthTracks = True
+ToolSvc += FTKSGInput
+
+print "Output file", OutputNTUP_FTKIPFile
+
+theJob += FTKRegionalWrapper(OutputLevel = DEBUG,
+ PMapPath = pmap_path,
+ RMapPath = rmap_path,
+ OutFileName = OutputNTUP_FTKIPFile)
+
+print theJob
+###############################################################
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_8Reg_jobOptions.py b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_8Reg_jobOptions.py
new file mode 100644
index 00000000000..b073d616091
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimRegionalWrap_8Reg_jobOptions.py
@@ -0,0 +1,37 @@
+###############################################################
+#
+# DumpSp job options file
+#
+#==============================================================
+
+# Helper function from transforms
+from PyJobTransforms.trfUtils import findFile
+pmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_11L.pmap')
+print "Using PMAP:", pmap_path
+rmap_path = findFile(os.environ['DATAPATH'], 'ftk_configuration/map_files/raw_11L.rmap')
+print "Using RMAP:", rmap_path
+
+from AthenaCommon.AlgSequence import AlgSequence
+theJob = AlgSequence()
+
+
+from FastTrackSimWrap.FastTrackSimWrapConf import FTKRegionalWrapper
+if hasattr(runArgs,"outputNTUP_FTKIPFile") :
+ OutputNTUP_FTKIPFile = runArgs.outputNTUP_FTKIPFile
+else :
+ OutputNTUP_FTKIPFile = "ftksim_8Reg_wrap.root"
+
+from TrigFTKSim.TrigFTKSimConf import FTK_SGHitInput
+FTKSGInput = FTK_SGHitInput()
+FTKSGInput.ReadTruthTracks = True
+ToolSvc += FTKSGInput
+
+print "Output file", OutputNTUP_FTKIPFile
+
+theJob += FTKRegionalWrapper(OutputLevel = DEBUG,
+ PMapPath = pmap_path,
+ RMapPath = rmap_path,
+ OutFileName = OutputNTUP_FTKIPFile)
+
+print theJob
+###############################################################
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimWrap_jobOptions.py b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimWrap_jobOptions.py
new file mode 100644
index 00000000000..ec952dc9fe0
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/share/FastTrackSimWrap_jobOptions.py
@@ -0,0 +1,45 @@
+###############################################################
+#
+# DumpSp job options file
+#
+#==============================================================
+# top algorithm to be run and component library to be loaded. if
+# InDetTruthToTrack does not exist, you haven't configured the
+# rest of the job correctly. see the job options inside
+# submit_torque.pl
+#theApp.TopAlg += [ "DumpSp" ]
+#theApp.Dlls += [ "FastTrackSimWrap" ]
+#DumpSp = Algorithm( "DumpSp" )
+#DumpSp.maxEta = 3.3
+#DumpSp.minPt = 0.*GeV
+#DumpSp.OutFileName = "ftksim_wrap.dat.bz2"
+#DumpSp.OutFileNameRawHits = "ftksim_wrap_raw.dat.bz2"
+#DumpSp.OutputLevel = 2
+#DumpSp.tracksName = "Tracks"
+
+
+from AthenaCommon.AlgSequence import AlgSequence
+theJob = AlgSequence()
+
+
+from FastTrackSimWrap.FastTrackSimWrapConf import DumpSp
+
+if hasattr(runArgs,"outputTXT_FTKIPFile") :
+ OutputTXT_FTKIPFile = runArgs.outputTXT_FTKIPFile
+else :
+ OutputTXT_FTKIPFile = 'ftksim_wrap_raw.dat.bz2'
+
+
+theJob += DumpSp(OutputLevel = INFO,
+ maxEta = 3.3,
+ minPt = 0.8*GeV,
+ DoData = False,
+ OutFileName = "ftksim_wrap.dat.bz2",
+ OutFileNameRawHits = OutputTXT_FTKIPFile,
+ tracksName = "Tracks")
+
+print theJob
+
+
+
+###############################################################
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/share/submit_torque.pl b/Trigger/TrigFTK/FastTrackSimWrap/share/submit_torque.pl
new file mode 100755
index 00000000000..fe68e2f787d
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/share/submit_torque.pl
@@ -0,0 +1,236 @@
+#!/usr/bin/perl -w
+
+use strict;
+use Getopt::Long;
+use FileHandle;
+
+my $datestring = `date --iso-8601`;
+chomp( ${datestring} );
+$datestring =~ s/\-//g;
+
+my $sample = "";
+my ${run_dir} = `pwd`; chomp( ${run_dir} );
+my $maxEvents = -1;
+my $bzip_output = '';
+my $result_basedir = "";
+my $result_rundir = "";
+&GetOptions('s=s'=>\$sample,'d=s'=>\${run_dir},'n=i'=>\${maxEvents},'r=s'=>\${result_basedir},'bzip'=>\${bzip_output},'l=s'=>\${result_rundir});
+
+die "provide sample name." unless (${sample} ne "");
+die "number of events doesn't make sense" unless( ${maxEvents}!=0 && ${maxEvents}>=-1 );
+
+my $result_subdir = "${sample}/${datestring}";
+my $result_dir;
+if( !($result_basedir eq "") ) {
+ $result_dir = "${result_basedir}/${result_subdir}";
+} elsif( !($run_dir eq "") ) {
+ $result_dir = "${run_dir}/results/${result_subdir}";
+} else {
+ $result_dir = "./${result_subdir}";
+}
+if( $result_rundir eq "" ) {
+ $result_rundir = "${run_dir}/submission/${result_subdir}";
+}
+
+print "running sample ${sample}\nin ${run_dir}\n";
+if( $maxEvents != -1 ) {
+ print "on ${maxEvents}\n";
+}
+if( !($bzip_output eq '') ) { print "with bzip2 compression\n"; }
+print "with output to ${result_dir}\n";
+print "and batch/log output to ${result_rundir}\n";
+
+if( !(-e $result_dir) ) { system( "mkdir -p $result_dir" ); }
+if( !(-e $result_rundir) ) { system( "mkdir -p $result_rundir" ); }
+
+my @filenames;
+
+# define filenames
+if( $sample eq "test" ) {
+ @filenames = ( "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9069.pool.root" );
+}
+if( $sample eq "v1213indep" ) {
+ @filenames = ( "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9069.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9070.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9071.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9073.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9074.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9078.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9079.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9080.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9081.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9082.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9083.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9084.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9085.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9086.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9087.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9088.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9089.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9090.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9091.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9092.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9093.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9094.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9095.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9096.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9097.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9098.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9099.pool.root" ,
+ "dcap:///pnfs/uct3/data/users/brubaker/ftkprod1213_080712/digi/g4digi-9100.pool.root" );
+}
+if( $sample eq "cscwhbb1034" ) {
+ @filenames = ( "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00001.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00004.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00005.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00006.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00007.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00008.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00009.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00010.pool.root.3" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00011.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00013.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00015.pool.root.3" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00019.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00021.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00024.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00025.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00026.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00031.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00033.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00035.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00040.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00042.pool.root.3" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00043.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00047.pool.root.3" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00048.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00049.pool.root.3" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00050.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00051.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00052.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00054.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00055.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00059.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00061.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00062.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00063.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00066.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00068.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00073.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00076.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00080.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00081.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00082.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00085.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00086.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00087.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00089.pool.root.3" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00090.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00091.pool.root.2" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00094.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00100.pool.root.1" );
+}
+if( $sample eq "cscwhbb1034anton" ) {
+ @filenames = ( #"dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00001.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00013.pool.root.1" ,
+ "dcap:///pnfs/uct3/users/antonk/specialg_misal1_csc11.005850.WH120bb_pythia.digit.RDO.v13004003/RDO.022779._00015.pool.root.3" );
+}
+
+
+# create one batch job for each file
+my $jobid = 0;
+my $skipEvents = 0;
+my $nsleepsec=0;
+foreach my $filename (@filenames) {
+ ++$jobid;
+ $nsleepsec += 5;
+ my $batchname = "${result_rundir}/batch_${jobid}.sh";
+ my $logname = "${result_rundir}/batch_${jobid}.log";
+ my $joboptionsname = "${result_rundir}/jobOptions_${jobid}.py";
+ my $rawoutputname = "${result_dir}/ftksimwrap_hit_raw_${jobid}.dat";
+ my $outputname = "${result_dir}/ftksimwrap_hit_${jobid}.dat";
+ if( ${bzip_output} ne '' ) {
+ $rawoutputname = "${rawoutputname}.bz2";
+ $outputname = "${outputname}.bz2";
+ }
+ # remove output files if they exist
+ if( -e $rawoutputname ) { system( "rm -f ${rawoutputname}" ); }
+ if( -e $outputname ) { system( "rm -f ${outputname}" ); }
+ if( $jobid==1 ) {
+ $skipEvents = 147;
+ }
+ if( $jobid==224 ) {
+ $skipEvents = 147;
+ }
+ # write torque job script
+ my $fhbatch = FileHandle::new();
+ open $fhbatch, ">${batchname}";
+ print $fhbatch <<ENDBATCHSCRIPT;
+#!/bin/bash
+sleep ${nsleepsec}s
+cd ${run_dir}
+cd ..
+pwd
+hostname
+source setup_athena.sh
+source setup_athena.sh
+cd ${run_dir}
+nohup athena.py ${joboptionsname}
+ENDBATCHSCRIPT
+ system( "chmod u+x $batchname" );
+ # write job options
+ my $fhjoboptions = FileHandle::new();
+ open $fhjoboptions, ">${joboptionsname}";
+ print $fhjoboptions <<ENDJOBOPTIONS;
+# job options for 13.0.40 simulation
+OutputLevel = WARNING
+doTrkNtuple = False
+doPixelTrkNtuple = False
+doJiveXML = False
+doVP1 = False
+doWriteESD = False
+doWriteAOD = False
+doReadBS = False
+doAuditors = True
+
+import os
+if os.environ['CMTCONFIG'].endswith('-dbg'):
+ # --- do EDM monitor (debug mode only)
+ doEdmMonitor = True
+ # --- write out a short message upon entering or leaving each algorithm
+ doNameAuditor = True
+else:
+ doEdmMonitor = False
+ doNameAuditor = False
+include( "InDetRecExample/ConfiguredInDetFlags.py")
+InDetFlags = ConfiguredInDetFlags(xKalman = True,
+ iPatRec = True,
+ newTracking = True,
+ doTruth = True)
+InDetFlags.enableBackTracking()
+InDetFlags.enableV0Finder()
+InDetFlags.enableConversions()
+from AthenaCommon.DetFlags import DetFlags
+DetDescrVersion = "ATLAS-CSC-01-02-00"
+from TrkExTools.AtlasExtrapolator import AtlasExtrapolator
+AtlasExtrapolator = AtlasExtrapolator(name='Trk::Extrapolator')
+ToolSvc += AtlasExtrapolator
+include("InDetRecExample/InDetRec_all.py")
+# run fasttracksimwrap
+theApp.TopAlg += [ "DumpSp" ]
+theApp.Dlls += [ "FastTrackSimWrap" ]
+DumpSp = Algorithm( "DumpSp" )
+DumpSp.OutFileName = "${outputname}"
+DumpSp.OutFileNameRawHits = "${rawoutputname}"
+DumpSp.OutputLevel = 4
+DumpSp.minPt = 0.*GeV
+DumpSp.maxEta = 3.3
+DumpSp.TruthToTrackTool = InDetTruthToTrack
+# input file
+ServiceMgr.EventSelector.InputCollections = [ "${filename}" ]
+theApp.EvtMax = $maxEvents
+ServiceMgr.EventSelector.SkipEvents = $skipEvents
+ENDJOBOPTIONS
+ #
+ system( "qsub -l walltime=23:59:59 -N wrapper_${sample} -j oe -o ${logname} ${batchname}" );
+} # end for each input filename
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/src/DumpSp.cxx b/Trigger/TrigFTK/FastTrackSimWrap/src/DumpSp.cxx
new file mode 100644
index 00000000000..1c3f83e851b
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/src/DumpSp.cxx
@@ -0,0 +1,1896 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+// FASTTRACKSIMWRAP
+// ================================================================
+// This algorithm is used to extract all needed information from
+// Athena for use in the standalone FTK simulation. Eventually an
+// Athena-based FTK simulation will exist and this code will
+// disappear.
+// ================================================================
+//
+// 2011-03-22 Antonio Boveia (boveia@hep.uchicago.edu)
+// Ability to write text files directly to bzip2 using
+// boost::iostreams.
+// 2010-03-25 Antonio Boveia
+// Add SCT cluster data to wrapper output.
+// 2010-03-10 Joe Tuggle (jtuggle@hep.uchicago.edu)
+// Add pixel cluster data
+// 2009-05-20 Antonio Boveia
+// access SDOs and dump GEANT truth for each silicon
+// channel. perform matching of offline reconstructed
+// tracks to truth using this info. substantial rewrite
+// and cleanup.
+// 2007-12-10 Giulio Usai <giulio.usai@cern.ch>
+// Monica Dunford <mdunford@hep.uchicago.edu>
+// original version.
+
+// include before any HAVE_VERISON_15 checks below.
+#include "FastTrackSimWrap/DumpSp.h"
+
+
+#include <cmath>
+#include <cstdlib>
+#include <set>
+#include <map>
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <iostream>
+#include <iomanip>
+#include <limits>
+#include <utility>
+#include <fstream>
+#include <boost/bind.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/format.hpp>
+#include <boost/foreach.hpp>
+#include <boost/iostreams/filter/bzip2.hpp>
+#include <boost/iostreams/filtering_stream.hpp>
+#include <boost/iostreams/device/file.hpp>
+#include <boost/iostreams/concepts.hpp>
+#include <boost/iostreams/categories.hpp>
+#include <boost/iostreams/device/file.hpp>
+#include <boost/algorithm/string/predicate.hpp>
+
+#include "HepMC/GenEvent.h"
+#include "HepMC/GenVertex.h"
+#include "HepMC/GenParticle.h"
+
+#include "AtlasDetDescr/AtlasDetectorID.h"
+#include "IdDictDetDescr/IdDictManager.h"
+#include "InDetIdentifier/PixelID.h"
+#include "InDetIdentifier/SCT_ID.h"
+#include "InDetIdentifier/TRT_ID.h"
+#include "HepPDT/ParticleDataTable.hh"
+#include "HepPDT/ParticleData.hh"
+#include "InDetSimData/InDetSimDataCollection.h"
+#include "InDetSimData/SCT_SimHelper.h"
+#include "InDetSimData/PixelSimHelper.h"
+#include "InDetReadoutGeometry/SiCellId.h"
+#include "InDetReadoutGeometry/SCT_ModuleSideDesign.h"
+#include "TrkRIO_OnTrack/RIO_OnTrack.h"
+
+#include "EventInfo/EventInfo.h"
+#include "EventInfo/EventID.h"
+#include "StoreGate/DataHandle.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "TrkSpacePoint/SpacePoint.h"
+#include "TrkSpacePoint/SpacePointCLASS_DEF.h"
+#include "TrkSpacePoint/SpacePointCollection.h"
+#include "TrkSpacePoint/SpacePointContainer.h"
+#include "TrkSpacePoint/SpacePointOverlapCollection.h"
+#include "Identifier/IdentifierHash.h"
+#include "Identifier/Identifier.h"
+#include "IdDict/IdDictDefs.h"
+#include "InDetPrepRawData/SiClusterContainer.h"
+#include "InDetPrepRawData/SiClusterCollection.h"
+#include "GeneratorObjects/HepMcParticleLink.h"
+#include "GaudiKernel/IPartPropSvc.h"
+#include "TrkToolInterfaces/ITruthToTrack.h"
+#include "TrigCaloEvent/TrigT2MbtsBits.h"
+#include "TileIdentifier/TileTBID.h"
+#include "InDetConditionsSummaryService/IInDetConditionsSvc.h"
+#include "VxVertex/VxContainer.h"
+//#include "TrkTruthToTrack/TruthToTrack.h"
+#include "TrkToolInterfaces/ITruthToTrack.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrkTruthData/TrackTruthCollection.h"
+#include "TrkTrackSummaryTool/TrackSummaryTool.h"
+#include "TrkToolInterfaces/ITrackHoleSearchTool.h"
+
+#include "InDetRawData/InDetRawDataCollection.h"
+#include "InDetRawData/InDetRawDataContainer.h"
+#include "InDetRawData/InDetRawDataCLASS_DEF.h"
+#include "InDetReadoutGeometry/SiDetectorElementCollection.h"
+#include "InDetReadoutGeometry/SiDetectorElement.h"
+#include "GeneratorObjects/McEventCollection.h"
+
+#include "InDetReadoutGeometry/SiDetectorManager.h"
+#include "InDetReadoutGeometry/PixelDetectorManager.h"
+#include "InDetReadoutGeometry/SCT_DetectorManager.h"
+#include "InDetReadoutGeometry/TRT_DetectorManager.h"
+#include "InDetBeamSpotService/IBeamCondSvc.h"
+
+#include "EventPrimitives/EventPrimitivesHelpers.h"
+
+using namespace std;
+
+DumpSp::DumpSp(const string& name, ISvcLocator* pSvcLocator)
+ : AthAlgorithm( name , pSvcLocator )
+ , m_idHelper()
+ , m_storeGate( 0 )
+ , m_detStore( 0 )
+ , m_evtStore( 0 )
+ , m_pixelId( 0 )
+ , m_sctId( 0 )
+ , m_trtId( 0 )
+ , m_PIX_mgr( 0 )
+ , m_SCT_mgr( 0 )
+ , m_TRT_mgr( 0 )
+ , m_pixelContainer( 0 )
+ , m_sctContainer( 0 )
+ , m_particleDataTable( 0 )
+ , m_truthToTrack( "Trk::TruthToTrack/InDetTruthToTrack" )
+ , m_trkSummaryTool( "Trk::TrackSummaryTool/InDetTrackSummaryTool" )
+ , m_trigDecTool("Trig::TrigDecisionTool/TrigDecisionTool")
+ , m_pixelCondSummarySvc("PixelConditionsSummarySvc",name)
+ , m_sctCondSummarySvc("SCT_ConditionsSummarySvc",name)
+ // , m_holeSearchTool( "InDetHoleSearchTool" )
+ , m_beamCondSvc("BeamCondSvc",name)
+ , m_pixelClustersName( "PixelClusters" )
+ , m_sctClustersName( "SCT_Clusters" )
+ , m_pixelSpacePointsName( "PixelSpacePoints" )
+ , m_sctSpacePointsName( "SCT_SpacePoints" )
+ , m_overlapSpacePointsName( "OverlapSpacePoints" )
+ , m_tracksName( "ExtendedTracks" )
+ , m_tracksTruthName( "ExtendedTracksTruthCollection" )
+ , m_maxEta( 3.3 )
+ , m_minPt( 0. )
+ , m_spacePointsName( "" )
+ , m_outFileName( "ftksim_hits.dat.bz2" )
+ , m_outFileNameRawHits( "ftksim_raw_hits.dat.bz2" )
+ , m_doTrigger (false)
+ , m_doData (false)
+ , m_doVertex (false)
+ , m_doLumiBCID (false)
+ , m_doBadMod( true )
+ , m_dumpHitsOnTracks( false )
+ , m_dumpSpacePoints( false )
+ , m_dumpTruthIntersections( false )
+ , m_dumpMBTSInfo( false )
+ , m_useOfflineTrackSelectorTool( false )
+ , m_outputBeamSpotToWrapper( false )
+ , m_useSimpleCuts( true )
+ , m_doRDODebug( false )
+ , ofl()
+ , oflraw()
+{
+ declareProperty("maxEta", m_maxEta);
+ declareProperty("minPt", m_minPt);
+ declareProperty("pixelClustersName", m_pixelClustersName);
+ declareProperty("SCT_ClustersName", m_sctClustersName);
+ declareProperty("overlapSpacePointsName", m_overlapSpacePointsName);
+ declareProperty("pixelSpacePointsName", m_pixelSpacePointsName);
+ declareProperty("sctSpacePointsName", m_sctSpacePointsName);
+ declareProperty("OutFileName", m_outFileName);
+ declareProperty("OutFileNameRawHits", m_outFileNameRawHits);
+ declareProperty("dumpHitsOnTracks", m_dumpHitsOnTracks);
+ declareProperty("dumpSpacePoints", m_dumpSpacePoints);
+ declareProperty("dumpTruthIntersections", m_dumpTruthIntersections);
+ declareProperty("SummaryTool", m_trkSummaryTool);
+ declareProperty("tracksName" , m_tracksName);
+ declareProperty("tracksTruthName" , m_tracksTruthName);
+ declareProperty("TruthToTrackTool" , m_truthToTrack);
+ // declareProperty("HoleSearch" , m_holeSearchTool);
+ declareProperty("beamCondSvc" , m_beamCondSvc);
+ declareProperty("dumpMBTSInfo" , m_dumpMBTSInfo , "Dump X fields to output containing MBTS trigger info" );
+ declareProperty("useOfflineTrackSelectorTool" , m_useOfflineTrackSelectorTool);
+ declareProperty("outputBeamSpotToWrapper" , m_outputBeamSpotToWrapper);
+ declareProperty("useSimpleCuts" , m_useSimpleCuts);
+ //#ifdef HAVE_VERSION_15
+ declareProperty("TrigDecisionTool", m_trigDecTool );
+ declareProperty("PixelSummarySvc" , m_pixelCondSummarySvc);
+ declareProperty("SctSummarySvc" , m_sctCondSummarySvc);
+ declareProperty("DoTrigger" , m_doTrigger);
+ declareProperty("DoData" , m_doData);
+ declareProperty("DoVertex" , m_doVertex);
+ declareProperty("DoLumiBCID" , m_doLumiBCID);
+ declareProperty("DoBadMod" , m_doBadMod);
+ declareProperty("DoRDODebug" , m_doRDODebug);
+ //#endif
+}
+
+StatusCode
+DumpSp::initialize()
+{
+
+ ATH_MSG_INFO("DumpSp::initialize()");
+
+ if( service("StoreGateSvc", m_storeGate).isFailure() ) {
+ ATH_MSG_FATAL("StoreGate service not found");
+ return StatusCode::FAILURE;
+ }
+
+ if( m_truthToTrack.retrieve().isFailure() ) {
+ ATH_MSG_FATAL(m_truthToTrack << " truth to track tool not found");
+ return StatusCode::FAILURE;
+ } else {
+ ATH_MSG_INFO(m_truthToTrack << " retrieved");
+ }
+
+ // if ( m_holeSearchTool.retrieve().isFailure() ) {
+ // ATH_MSG_FATAL("Failed to retrieve tool " << m_holeSearchTool);
+ // return StatusCode::FAILURE;
+ // }
+
+ if( service("DetectorStore",m_detStore).isFailure() ) {
+ ATH_MSG_FATAL("DetectorStore service not found");
+ return StatusCode::FAILURE;
+ }
+
+ IPartPropSvc* partPropSvc = 0;
+ if( service("PartPropSvc", partPropSvc, true).isFailure() ) {
+ ATH_MSG_FATAL("particle properties service unavailable");
+ return StatusCode::FAILURE;
+ }
+ m_particleDataTable = partPropSvc->PDT();
+
+ // ID helpers
+ m_idHelper.reset( new AtlasDetectorID );
+ const IdDictManager* idDictMgr( 0 );
+ if( m_detStore->retrieve(idDictMgr, "IdDict").isFailure() || !idDictMgr ) {
+ ATH_MSG_ERROR( "Could not get IdDictManager !");
+ return StatusCode::FAILURE;
+ }
+ if( idDictMgr->initializeHelper( *m_idHelper.get() ) ) { // Returns 1 if there is a problem
+ ATH_MSG_ERROR( "Unable to initialize ID helper.");
+ return StatusCode::FAILURE;
+ }
+ if( m_detStore->retrieve(m_PIX_mgr, "Pixel").isFailure() ) {
+ ATH_MSG_ERROR( "Unable to retrieve Pixel manager from DetectorStore");
+ return StatusCode::FAILURE;
+ }
+ if( m_detStore->retrieve(m_pixelId, "PixelID").isFailure() ) {
+ ATH_MSG_ERROR( "Unable to retrieve Pixel helper from DetectorStore");
+ return StatusCode::FAILURE;
+ }
+ if( m_detStore->retrieve(m_SCT_mgr, "SCT").isFailure() ) {
+ ATH_MSG_ERROR( "Unable to retrieve SCT manager from DetectorStore");
+ return StatusCode::FAILURE;
+ }
+ if( m_detStore->retrieve(m_sctId, "SCT_ID").isFailure() ) {
+ ATH_MSG_ERROR( "Unable to retrieve SCT helper from DetectorStore");
+ return StatusCode::FAILURE;
+ }
+ // if( m_detStore->retrieve(m_TRT_mgr, "TRT").isFailure() ) {
+ // ATH_MSG_ERROR( "Unable to retrieve TRT manager from DetectorStore");
+ // return StatusCode::FAILURE;
+ // }
+
+ if( m_trkSummaryTool.retrieve().isFailure() ) {
+ ATH_MSG_FATAL("Failed to retrieve tool " << m_trkSummaryTool);
+ return StatusCode::FAILURE;
+ }
+
+ if(m_doTrigger || m_dumpMBTSInfo ){
+ if( m_trigDecTool.retrieve().isFailure() ) {
+ ATH_MSG_ERROR( "Unable to retrieve TrigDecsionTool");
+ } else {
+ ATH_MSG_INFO("Successfully retrieved TrigDecisionTool");
+ }
+ }
+ if( m_doBadMod) {
+ // Get PixelConditionsSummarySvc
+ if ( m_pixelCondSummarySvc.retrieve().isFailure() ) {
+ ATH_MSG_FATAL("Failed to retrieve tool " << m_pixelCondSummarySvc);
+ return StatusCode::FAILURE;
+ } else {
+ ATH_MSG_INFO("Retrieved tool " << m_pixelCondSummarySvc);
+ }
+ // Get SctConditionsSummarySvc
+ if ( m_sctCondSummarySvc.retrieve().isFailure() ) {
+ ATH_MSG_FATAL("Failed to retrieve tool " << m_sctCondSummarySvc);
+ return StatusCode::FAILURE;
+ } else {
+ ATH_MSG_INFO("Retrieved tool " << m_sctCondSummarySvc);
+ }
+ }
+
+
+ // open output streams
+ if( m_dumpSpacePoints ) {
+ ofl.reset( new boost::iostreams::filtering_ostream );
+ if( !ofl ) { return StatusCode::FAILURE; }
+ if( boost::algorithm::icontains( m_outFileName , ".bz2" ) ) {
+ boost::iostreams::bzip2_params params;
+ params.block_size = 9;
+ ofl->push( boost::iostreams::bzip2_compressor(params) );
+ }
+ ofl->push( boost::iostreams::file_sink(m_outFileName) ); // open the file
+ }
+ if( true ) {
+ oflraw.reset( new boost::iostreams::filtering_ostream );
+ if( !oflraw ) { return StatusCode::FAILURE; }
+ if( boost::algorithm::icontains( m_outFileNameRawHits , ".bz2" ) ) {
+ boost::iostreams::bzip2_params params;
+ params.block_size = 9;
+ oflraw->push( boost::iostreams::bzip2_compressor(params) );
+ }
+ oflraw->push( boost::iostreams::file_sink(m_outFileNameRawHits) ); // open the file
+ }
+
+ // jordan's code for the beamline
+ if ( m_beamCondSvc.retrieve().isFailure() && ( m_outputBeamSpotToWrapper || !m_useSimpleCuts ) ) {
+ ATH_MSG_FATAL("Failed to retrieve service " << m_beamCondSvc);
+ return StatusCode::FAILURE;
+ } else {
+ ATH_MSG_INFO("Retrieved service " << m_beamCondSvc);
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+DumpSp::execute()
+{
+ const EventInfo* eventInfo(0);
+ if( m_storeGate->retrieve(eventInfo).isFailure() ) {
+ ATH_MSG_ERROR( "Could not retrieve event info");
+ return StatusCode::FAILURE;
+ }
+ const EventID* eventID( eventInfo->event_ID() );
+
+ ATH_MSG_DEBUG( "entered execution for run " << eventID->run_number()
+ << " event " << eventID->event_number());
+
+ if(m_doTrigger || m_doData) {
+ //make trigger selection
+ const bool ok_passRandom = m_trigDecTool->isPassed("EF_rd0_filled_NoAlg");
+ if(!ok_passRandom){
+ ATH_MSG_DEBUG( " Event failed trigger requirement ");
+ return StatusCode::SUCCESS;
+ }else ATH_MSG_DEBUG( " Event passed trigger requirement! ");
+
+ }
+
+
+
+ // event header
+ if( m_dumpSpacePoints ) {
+ (*ofl) << "R\t" << eventID->run_number()<<'\n';
+ (*ofl) << "F\t" << eventID->event_number()<<'\t'<<eventInfo->averageInteractionsPerCrossing()<<'\t'<<eventInfo->actualInteractionsPerCrossing()<< '\n';
+ }
+ (*oflraw) << "R\t" << eventID->run_number()<<'\n';
+ (*oflraw) << "F\t" << eventID->event_number()<<'\t'<<eventInfo->averageInteractionsPerCrossing()<<'\t'<<eventInfo->actualInteractionsPerCrossing()<<'\n';
+ {
+ // dump bad modules once at the start of each lumi block.
+ static unsigned int last_lumi_block = static_cast<unsigned int>(~0);
+ if( last_lumi_block==static_cast<unsigned int>(~0) || eventID->lumi_block()!=last_lumi_block ) {
+ dump_bad_modules();
+ last_lumi_block = eventID->lumi_block();
+ }
+ }
+
+ // build truth-track geant matching
+
+ if( !m_doData)
+ build_matching_maps();
+
+ HitIndexMap hitIndexMap; // keep running index event-unique to each hit
+ HitIndexMap clusterIndexMap;
+ // get pixel and sct cluster containers
+ if( m_storeGate->retrieve(m_pixelContainer, m_pixelClustersName).isFailure() ) {
+ ATH_MSG_WARNING( "unable to retrieve the PixelCluster container " << m_pixelClustersName);
+ }
+ if( m_storeGate->retrieve(m_sctContainer, m_sctClustersName).isFailure() ) {
+ ATH_MSG_WARNING( "unable to retrieve the SCT_Cluster container " << m_sctClustersName);
+ }
+
+ if( m_dumpSpacePoints ) {
+ dump_spacepoints();
+ }
+
+ // dump raw silicon data
+ dump_raw_silicon( hitIndexMap , clusterIndexMap );
+
+ // dump ipat tracks using new format and interpretation code
+ dump_tracks( hitIndexMap , clusterIndexMap );
+
+ // dump truth tracks
+ if( !m_doData)
+ dump_truth();
+
+ // dump beamspot to wrapper output
+ dump_beamspot();
+
+
+ // dump data for occupancy studies
+ if( m_dumpMBTSInfo ) { dump_MBTS( ); }
+ // dump vertex data for occupancy studies
+ if( m_doLumiBCID || m_doData) { dump_lumibcid( eventID ); }
+ // dump vertex data for occupancy studies
+ if( m_doVertex || m_doData) { dump_vertex( ); }
+
+
+ // event footer
+ if( m_dumpSpacePoints ) { (*ofl) << "L\t" << eventID->event_number()<<'\n'; }
+ (*oflraw) << "L\t" << eventID->event_number()<<'\n';
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode
+DumpSp::finalize()
+{
+
+ ATH_MSG_INFO("finalized");
+
+ return StatusCode::SUCCESS;
+}
+
+// build ttr and rtt truth matching maps for this event
+
+void
+DumpSp::build_matching_maps()
+{
+ _ttrTrackMap.clear();
+ _ttrProbMap.clear();
+ _rttTrackMap.clear();
+
+ ATH_MSG_DEBUG( "building reconstruction-truth matching maps");
+
+ // retrieve necessary junk from Athena
+ const TrackCollection* RecCollection = 0;
+ if( m_storeGate->retrieve( RecCollection, m_tracksName ).isFailure() ) {
+ ATH_MSG_WARNING( "could not find TrackCollection " << m_tracksName);
+ RecCollection = 0;
+ }
+
+ const DataVector<Trk::Track>* trks = 0;
+ if( m_storeGate->retrieve(trks,m_tracksName).isFailure() ) {
+ ATH_MSG_WARNING( "could not find Trk::Track collection " << m_tracksName);
+ trks = 0;
+ }
+
+ const McEventCollection* SimTracks = 0;
+ if( m_storeGate->retrieve(SimTracks,"TruthEvent").isFailure() ) {
+ string key = "G4Truth";
+ if( m_storeGate->retrieve(SimTracks,key).isFailure() ) {
+ key = "";
+ if( m_storeGate->retrieve(SimTracks,key).isFailure() ) {
+ ATH_MSG_WARNING( "could not find the McEventCollection");
+ return;
+ }
+ }
+ }
+
+ if( !m_particleDataTable ) {
+ ATH_MSG_WARNING( " could not find the particle properties service");
+ }
+
+ const TrackTruthCollection* TruthMap = 0;
+ if( m_storeGate->retrieve(TruthMap,m_tracksTruthName).isFailure() ) {
+ ATH_MSG_WARNING( "could not find truth map " << m_tracksTruthName);
+ TruthMap = 0;
+ }
+
+ if( !RecCollection || !trks || !TruthMap ) {
+ ATH_MSG_WARNING( "missing necessary info for truth matching");
+ return;
+ }
+
+ // construct reconstruction to truth map.
+ for( DataVector<Trk::Track>::const_iterator trksItr=trks->begin(), trksItrF=trks->end(); trksItr!=trksItrF; ++trksItr ) {
+ // retrieve data
+ ElementLink<TrackCollection> tracklink;
+ tracklink.setElement( const_cast<Trk::Track*>(*trksItr) );
+ tracklink.setStorableObject( *RecCollection );
+ const ElementLink<TrackCollection> tracklink2( tracklink );
+ TrackTruthCollection::const_iterator found( TruthMap->find(tracklink2) );
+ if( found==TruthMap->end() ) { continue; }
+ TrackTruth trtruth( found->second );
+ HepMcParticleLink::ExtendedBarCode extBarcode(trtruth.particleLink().barcode(),trtruth.particleLink().eventIndex());
+ // update _ttrTrackMap with track index corresponding to the greatest figure of merit.
+ if( _ttrProbMap.find(extBarcode)==_ttrProbMap.end() ) {
+ // this is the only track matching this barcode so far
+ _ttrTrackMap.insert( pair<HepMcParticleLink::ExtendedBarCode,unsigned int>(extBarcode,distance(trks->begin(),trksItr)) );
+ } else {
+ // a reconstructed track match already exists. figure if this match is better, and keep it if so.
+ //
+ // retrieve probabilities for this barcode to match to each track
+ static vector<float> probs;
+ probs.clear();
+ transform( _ttrProbMap.lower_bound(extBarcode) , _ttrProbMap.upper_bound(extBarcode) , back_inserter(probs) ,
+ boost::bind(&TruthToRecoProbMap::value_type::second,_1) );
+ // determine track with highest figure-of-merit match and update
+ // the map of barcodes => reconstructed track indices (into
+ // track collection)
+ vector<float>::const_iterator imax = max_element(probs.begin(),probs.end());
+ if( imax!=probs.end() && (trtruth.probability() > (*imax)) ) {
+ _ttrTrackMap[extBarcode] = distance(trks->begin(),trksItr);
+ }
+ }
+ // update map from barcodes to highest reconstructed track matching probabilities
+ _ttrProbMap.insert( pair<HepMcParticleLink::ExtendedBarCode,float>(extBarcode,trtruth.probability()) );
+ // fill reco to truth map. allow multiple reconstructed tracks to point to the same barcode.
+ _rttTrackMap[ distance(trks->begin(),trksItr) ] = extBarcode;
+ } // end for each reconstructed track
+
+ ATH_MSG_DEBUG( "truth info from " << trks->size() << " tracks.");
+ ATH_MSG_DEBUG( "_rttTrackMap constructed with " << _rttTrackMap.size() << " entries.");
+ ATH_MSG_DEBUG( "_ttrTrackMap constructed with " << _ttrTrackMap.size() << " entries.");
+ ATH_MSG_DEBUG( "_ttrProbMap constructed with " << _ttrProbMap.size() << " entries.");
+
+}
+
+// dump truth, including geant matching data for offline reconstructed track collection
+
+void
+DumpSp::dump_truth() const
+{
+ ATH_MSG_DEBUG( "dumping truth information");
+
+ // retrieve truth tracks from athena
+ const McEventCollection* SimTracks = 0;
+ if( m_storeGate->retrieve(SimTracks,"TruthEvent").isFailure() ) {
+ string key = "G4Truth";
+ if( m_storeGate->retrieve(SimTracks,key).isFailure() ) {
+ key = "";
+ if( m_storeGate->retrieve(SimTracks,key).isFailure() ) {
+ ATH_MSG_WARNING( "could not find the McEventCollection");
+ return;
+ }
+ }
+ }
+
+ // is PDG info is available?
+ if( !m_particleDataTable ) {
+ ATH_MSG_ERROR( " could not find the particle properties service");
+ return;
+ }
+
+ // for debugging the calculation of d0_corr
+
+ bool showd0corrSuccess = true;
+
+ // dump each truth track
+
+ for( unsigned int ievt=0, fevt=(SimTracks ? SimTracks->size() : 0u); ievt!=fevt; ++ievt ) {
+ const HepMC::GenEvent* genEvent = SimTracks->at( ievt );
+ // retrieve the primary interaction vertex here. for now, use the dummy origin.
+ Amg::Vector3D primaryVtx(0.,0.,0.);
+ // the event should have signal process vertex unless it was generated as single particles.
+ // if it exists, use it for the primary vertex.
+ if( genEvent->signal_process_vertex() ) {
+ primaryVtx << genEvent->signal_process_vertex()->point3d().x() , genEvent->signal_process_vertex()->point3d().y() , genEvent->signal_process_vertex()->point3d().z();
+ ATH_MSG_INFO("using signal process vertex for eventIndex " << ievt << ":"
+ << primaryVtx.x() << "\t" << primaryVtx.y() << "\t" << primaryVtx.z()
+ );
+ }
+ for( HepMC::GenEvent::particle_const_iterator it=genEvent->particles_begin(), ft=genEvent->particles_end(); it!=ft; ++it ) {
+ const HepMC::GenParticle* const particle( *it );
+ const int pdgcode = particle->pdg_id();
+ // reject generated particles without a production vertex.
+ if( !particle->production_vertex() ) { continue; }
+ // reject neutral or unstable particles
+ const HepPDT::ParticleData* pd = m_particleDataTable->particle(abs(pdgcode));
+ if( !pd ) { continue; }
+ const float charge = pd->charge();
+ if( std::abs(charge)<0.5 ) { continue; }
+ if( particle->status()%1000!=1 ) { continue; }
+ // retrieve some particle kinematics
+ const float genPt = particle->momentum().perp()/1000.; // convert to pt in GeV.
+ const float genEta = particle->momentum().pseudoRapidity();
+ // reject non-fiducial particles
+ if( (genPt*1000.) < m_minPt ) { continue; }
+ if( std::abs(genEta) > m_maxEta ) { continue; }
+ // retrieve truth track parameters at perigee
+ boost::scoped_ptr<const Trk::TrackParameters> generatedTrackPerigee( m_truthToTrack->makePerigeeParameters(particle) );
+ const float m_track_truth_d0 = generatedTrackPerigee ? generatedTrackPerigee->parameters()[Trk::d0] : 999.;
+ const float m_track_truth_phi = generatedTrackPerigee ? generatedTrackPerigee->parameters()[Trk::phi0] : 999.;
+ const float m_track_truth_p = (generatedTrackPerigee && generatedTrackPerigee->parameters()[Trk::qOverP] != 0.) ?
+ generatedTrackPerigee->charge()/generatedTrackPerigee->parameters()[Trk::qOverP] : 10E7;
+ const float m_track_truth_x0 = generatedTrackPerigee ? generatedTrackPerigee->position().x() : 999.;
+ const float m_track_truth_y0 = generatedTrackPerigee ? generatedTrackPerigee->position().y() : 999.;
+ const float m_track_truth_z0 = generatedTrackPerigee ? generatedTrackPerigee->position().z() : 999.;
+ const float m_track_truth_q = generatedTrackPerigee ? generatedTrackPerigee->charge() : 0.;
+ const float m_track_truth_sinphi = generatedTrackPerigee ? std::sin(generatedTrackPerigee->parameters()[Trk::phi0]) : -1.;
+ const float m_track_truth_cosphi = generatedTrackPerigee ? std::cos(generatedTrackPerigee->parameters()[Trk::phi0]) : -1.;
+ const float m_track_truth_sintheta = generatedTrackPerigee ? std::sin(generatedTrackPerigee->parameters()[Trk::theta]) : -1.;
+ const float m_track_truth_costheta = generatedTrackPerigee ? std::cos(generatedTrackPerigee->parameters()[Trk::theta]) : -1.;
+ float truth_d0corr = m_track_truth_d0-( primaryVtx.y()*cos(m_track_truth_phi)-primaryVtx.x()*sin(m_track_truth_phi) );
+ float truth_zvertex = 0.;
+ if ( !m_useSimpleCuts ) { // determine d0_corr based on beam position from BeamCondSvc
+ truth_d0corr = m_track_truth_d0-( m_beamCondSvc->beamPos().y()*cos(m_track_truth_phi)-m_beamCondSvc->beamPos().x()*sin(m_track_truth_phi) );
+ truth_zvertex = m_beamCondSvc->beamPos().z();
+ if ( showd0corrSuccess ) {
+ ATH_MSG_DEBUG( "Beamspot from BeamCondSvc used to determine cuts in dump_truth()");
+ showd0corrSuccess = false;
+ }
+ }
+ const Amg::Vector3D startVertex(particle->production_vertex()->point3d().x(), particle->production_vertex()->point3d().y(), particle->production_vertex()->point3d().z());
+ // categorize particle (prompt, secondary, etc.) based on InDetPerformanceRTT/detector paper criteria.
+ bool isPrimary = true;
+ bool isDetPaperCut = true;
+ if( std::abs(truth_d0corr)>2. ) { isPrimary=false; }
+ if( particle->barcode()>100000 || particle->barcode()==0 ) { isPrimary=false; }
+ isDetPaperCut=isPrimary;
+ // debug line for production_vertex variable
+ if( false ) { ATH_MSG_DEBUG( "z correction -- bool particle->production_vertex() = " << particle->production_vertex()); }
+ //
+ if( isPrimary && particle->production_vertex() ) {
+ const Amg::Vector3D startVertex(particle->production_vertex()->point3d().x(), particle->production_vertex()->point3d().y(), particle->production_vertex()->point3d().z());
+ if( std::abs(startVertex.z() - truth_zvertex)>100. ) { isPrimary=false; }
+ // debug z vertex correction
+ if( false ) { ATH_MSG_DEBUG( "z correction -- startVertex.z() = " << startVertex.z() << " , truth_zvertex = "<< truth_zvertex); }
+ //double flight_length = -1.;
+ if( particle->end_vertex() ) {
+ Amg::Vector3D endVertex(particle->end_vertex()->point3d().x(), particle->end_vertex()->point3d().y(), particle->end_vertex()->point3d().z());
+ if( endVertex.perp()<400. && std::abs(endVertex.z())<2300. ) { isPrimary=false; }
+ //flight_length = startVertex.distance( endVertex );
+ }
+ } else {
+ isPrimary = false;
+ }
+ isDetPaperCut=isPrimary;
+ // print truth track info and geant matching for highest figure-of-merit track
+ int irecmatch = -1;
+ float precmatch = 0.;
+ HepMcParticleLink::ExtendedBarCode extBarcode2( particle->barcode(), ievt );
+ if( !_ttrProbMap.empty() ) {
+ TruthToRecoProbMap::const_iterator barcode=_ttrProbMap.find(extBarcode2);
+ if( barcode!=_ttrProbMap.end() ) {
+ vector<float> probs;
+ transform( _ttrProbMap.lower_bound(extBarcode2) , _ttrProbMap.upper_bound(extBarcode2) , back_inserter(probs) ,
+ boost::bind(&TruthToRecoProbMap::value_type::second,_1) );
+ vector<float>::const_iterator imax = max_element(probs.begin(),probs.end());
+ assert( imax!=probs.end() );
+ precmatch = *imax;
+ TruthToRecoTrackMap::const_iterator ibestrec = _ttrTrackMap.find( extBarcode2 );
+ assert( ibestrec!=_ttrTrackMap.end() );
+ irecmatch = ibestrec->second;
+ }
+ }
+ ParentBitmask parent_mask( construct_truth_bitmap( particle ) );
+ (*oflraw) << setiosflags(ios::scientific) << "T\t"
+ << setw(14) << setprecision(10) << m_track_truth_x0 << '\t'
+ << setw(14) << setprecision(10) << m_track_truth_y0 << '\t'
+ << setw(14) << setprecision(10) << m_track_truth_z0 << '\t'
+ << (int)m_track_truth_q << '\t'
+ << setw(14) << setprecision(10) << m_track_truth_p*(m_track_truth_cosphi*m_track_truth_sintheta) << '\t'
+ << setw(14) << setprecision(10) << m_track_truth_p*(m_track_truth_sinphi*m_track_truth_sintheta) << '\t'
+ << setw(14) << setprecision(10) << m_track_truth_p*(m_track_truth_costheta) << '\t'
+ << pdgcode << '\t'
+ << setw(14) << (int)irecmatch << '\t'
+ << setw(14) << setprecision(10) << precmatch << '\t'
+ << extBarcode2.eventIndex() << '\t'
+ << extBarcode2.barcode() << '\t'
+ << parent_mask.to_ulong() << '\t'
+ << isDetPaperCut << '\t'
+ << resetiosflags(ios::scientific) << '\n';
+ if( false && particle->production_vertex() && (particle->momentum().perp()/1000.)>10. &&
+ particle->status()%1000==1 && std::abs(charge)>0.5 ) {
+ using boost::format;
+ const HepMC::GenParticle* tau_parent = 0;
+ vector<const HepMC::GenParticle*> parents;
+ parents.push_back( particle );
+ while( !parents.empty() ) {
+ const HepMC::GenParticle* p = parents.back();
+ if( std::abs(p->pdg_id())==15 ) { tau_parent = p; break; }
+ parents.pop_back();
+ if( !(p->production_vertex()) ) { continue; }
+ copy( p->production_vertex()->particles_begin(HepMC::parents) , p->production_vertex()->particles_end(HepMC::parents) , back_inserter(parents) );
+ }
+ unsigned int nprongs = 0;
+ if( tau_parent && tau_parent->end_vertex() ) { nprongs = tau_parent->end_vertex()->particles_out_size(); }
+ cout << " Truth Track: "
+ << format("%|10d| %|10d| %|10d| %|2d| %|2d|") % ievt % particle->barcode() % particle->pdg_id() % (static_cast<bool>(tau_parent)) % nprongs
+ << endl;
+ if( std::abs(particle->pdg_id())>=11 && std::abs(particle->pdg_id())<=16 ) { cout << "a lepton!" << endl; }
+ }
+ } // end for each GenParticle in this GenEvent
+ } // end for each GenEvent
+}
+
+
+// dump spacepoints. does not provide geant truth matching data.
+
+void
+DumpSp::dump_spacepoints() const
+{
+ // dump spacepoints for each of the relevant detectors.
+ unsigned int nsp( 0u );
+ const SpacePointContainer* pixelSPContainer(0);
+ const SpacePointContainer* sctSPContainer(0);
+ const SpacePointOverlapCollection* overlapCollection(0);
+ if( m_storeGate->retrieve(pixelSPContainer,m_pixelSpacePointsName).isFailure() ) {
+ ATH_MSG_DEBUG( "Unable to retrieve PixelSpacePoint container");
+ }
+ if( m_storeGate->retrieve(sctSPContainer,m_sctSpacePointsName).isFailure() ) {
+ ATH_MSG_DEBUG( "Unable to retrieve SCT_SpacePoint container");
+ }
+ if( m_storeGate->retrieve(overlapCollection,m_overlapSpacePointsName).isFailure() ) {
+ ATH_MSG_DEBUG( "Unable to retrieve Overlap SpacePoint container");
+ }
+ if( !(pixelSPContainer || sctSPContainer || overlapCollection) ) {
+ ATH_MSG_WARNING( "Unable to retrieve any SpacePoint containers");
+ } else {
+ if( pixelSPContainer ) {
+ for( SpacePointContainer::const_iterator i=pixelSPContainer->begin(), f=pixelSPContainer->end(); i!=f; ++i ) {
+ for( SpacePointCollection::const_iterator j=(**i).begin(), jf=(**i).end(); j!=jf; ++j ) {
+ const Trk::SpacePoint* sp( *j );
+ ++nsp;
+ // dump pixel info
+ IdentifierHash idh = (*sp).elementIdList().first;
+ Identifier pid = m_pixelId->wafer_id( idh );
+ const Amg::Vector3D point = (*sp).globalPosition();
+ const double rad = point.perp();
+ const double phi = point.phi();
+ const double z = point.z();
+ const double x = rad*cos(phi);
+ const double y = rad*sin(phi);
+ (*ofl) << "H\t"
+ << setw(14) << setprecision(10) << x << '\t'
+ << setw(14) << setprecision(10) << y << '\t'
+ << setw(14) << setprecision(10) << z << '\t'
+ << 1 << '\t' // 1 pixel 0 sct
+ << m_pixelId->barrel_ec(pid) << '\t'
+ << m_pixelId->layer_disk(pid) << '\t'
+ << m_pixelId->phi_module(pid) << '\t'
+ << m_pixelId->eta_module(pid) << '\t'
+ << m_pixelId->phi_index(pid) << '\t'
+ << m_pixelId->eta_index(pid) << '\n';
+ } // end of pixel space point loop
+ } // end of pixel container loop
+ } // end of pixel info
+ if( sctSPContainer ) {
+ for( SpacePointContainer::const_iterator i=sctSPContainer->begin(), f=sctSPContainer->end(); i!=f; ++i ) {
+ for (SpacePointCollection::const_iterator j=(**i).begin(), jf=(**i).end(); j!=jf; ++j ) {
+ const Trk::SpacePoint* sp( *j );
+ ++nsp;
+ IdentifierHash idh = (*sp).elementIdList().first;
+ Identifier pid = m_sctId->wafer_id( idh );
+ const Amg::Vector3D point = (*sp).globalPosition();
+ const double rad = point.perp();
+ const double phi = point.phi();
+ const double z = point.z();
+ const double x = rad*cos(phi);
+ const double y = rad*sin(phi);
+ (*ofl) << "H\t"
+ << setw(14) << setprecision(10) << x << '\t'
+ << setw(14) << setprecision(10) << y << '\t'
+ << setw(14) << setprecision(10) << z << '\t'
+ << 0 << '\t' // 1 pixel 0 sct
+ << m_sctId->barrel_ec(pid) << '\t'
+ << m_sctId->layer_disk(pid) << '\t'
+ << m_sctId->phi_module(pid) << '\t'
+ << m_sctId->eta_module(pid) << '\t'
+ << m_sctId->side(pid) << '\t'
+ << m_sctId->strip(pid) << '\n';
+ } // end of sct space point loop
+ } // end of sct container loop
+ } // end of sct info
+ if( overlapCollection ) {
+ for( SpacePointCollection::const_iterator i=overlapCollection->begin(), f=overlapCollection->end(); i!=f; ++i ) {
+ const Trk::SpacePoint* sp( *i );
+ ++nsp;
+ IdentifierHash idh = (*sp).elementIdList().first;
+ Identifier pid = m_sctId->wafer_id( idh );
+ const Amg::Vector3D point = (*sp).globalPosition();
+ const double rad = point.perp();
+ const double phi = point.phi();
+ const double z = point.z();
+ const double x = rad*cos(phi);
+ const double y = rad*sin(phi);
+ (*ofl) << "h\t"
+ << setw(14) << setprecision(10) << x << '\t'
+ << setw(14) << setprecision(10) << y << '\t'
+ << setw(14) << setprecision(10) << z << '\t'
+ << 0 << '\t' // 1 pixel 0 sct
+ << m_sctId->barrel_ec(pid) << '\t'
+ << m_sctId->layer_disk(pid) << '\t'
+ << m_sctId->phi_module(pid) << '\t'
+ << m_sctId->eta_module(pid) << '\t'
+ << m_sctId->side(pid) << '\t'
+ << m_sctId->strip(pid) << '\n';
+ } // end of overlap space point loop
+ } // end of SCT overlap info
+ } // dump all spacepoints
+ ATH_MSG_INFO(nsp << " SpacePoints found");
+}
+
+// dump silicon channels with geant matching information.
+
+void
+DumpSp::dump_raw_silicon( HitIndexMap& hitIndexMap, HitIndexMap& clusterIndexMap ) const
+{
+ using namespace std;
+ unsigned int hitIndex = 0u;
+ unsigned int clusterIndex = 0u;
+ // compute some statistics
+ // unsigned long nchannels = 0ul;
+ // unsigned long nchannels_single_barcodes = 0ul;
+ // unsigned long nchannels_single_barcodes_primary = 0ul;
+ // unsigned long nchannels_single_barcodes_pt1gev = 0ul;
+ // unsigned long nchannels_single_barcodes_bjet = 0ul;
+ // unsigned long nchannels_multiple_barcodes = 0ul;
+ // unsigned long nchannels_multiple_barcodes_primary = 0ul;
+ // unsigned long nchannels_multiple_barcodes_pt1gev = 0ul;
+ // unsigned long nchannels_multiple_barcodes_bjet = 0ul;
+
+ const DataHandle<PixelRDO_Container> pixel_rdocontainer_iter;
+ const InDetSimDataCollection* pixelSimDataMap(0);
+ const bool have_pixel_sdo = m_storeGate->retrieve(pixelSimDataMap, "PixelSDO_Map").isSuccess();
+ if( m_storeGate->retrieve(pixel_rdocontainer_iter, "PixelRDOs").isSuccess() ) {
+ pixel_rdocontainer_iter->clID(); // anything to dereference the DataHandle
+ for( PixelRDO_Container::const_iterator iColl=pixel_rdocontainer_iter->begin(), fColl=pixel_rdocontainer_iter->end(); iColl!=fColl; ++iColl ) {
+ const InDetRawDataCollection<PixelRDORawData>* pixel_rdoCollection(*iColl);
+ if( !pixel_rdoCollection ) { continue; }
+ const int size = pixel_rdoCollection->size();
+ ATH_MSG_DEBUG( "Pixel InDetRawDataCollection found with " << size << " RDOs");
+ // loop on all RDOs
+ for( DataVector<PixelRDORawData>::const_iterator iRDO=pixel_rdoCollection->begin(), fRDO=pixel_rdoCollection->end(); iRDO!=fRDO; ++iRDO ) {
+ Identifier rdoId = (*iRDO)->identify();
+ // get the det element from the det element collection
+ const InDetDD::SiDetectorElement* sielement = m_PIX_mgr->getDetectorElement(rdoId); assert( sielement );
+ const InDetDD::SiLocalPosition localPos = sielement->localPositionOfCell(rdoId);
+ const InDetDD::SiLocalPosition rawPos = sielement->rawLocalPositionOfCell(rdoId);
+ const Amg::Vector3D gPos( sielement->globalPosition(localPos) );
+ // update map between pixel identifier and event-unique hit index.
+ // ganged pixels (nCells==2) get two entries.
+ hitIndexMap[rdoId] = hitIndex;
+ const int nCells = sielement->numberOfConnectedCells( sielement->cellIdOfPosition(rawPos) );
+ if( nCells==2 ) {
+ const InDetDD::SiCellId firstCell = sielement->cellIdOfPosition(rawPos);
+ const InDetDD::SiCellId tmpCell = sielement->connectedCell(firstCell,1);
+ const Identifier tmpId = sielement->identifierFromCellId(tmpCell);
+ hitIndexMap[tmpId] = hitIndex; // add second entry for ganged pixel ID
+ }
+ // if there is simulation truth available, try to retrieve the "most likely" barcode for this pixel.
+ const HepMC::GenParticle* best_parent = 0;
+ ParentBitmask parent_mask;
+ HepMcParticleLink::ExtendedBarCode best_extcode;
+ if( have_pixel_sdo && pixelSimDataMap ) {
+ InDetSimDataCollection::const_iterator iter( pixelSimDataMap->find(rdoId) );
+ // this might be the ganged pixel copy.
+ if( nCells>1 && iter==pixelSimDataMap->end() ) {
+ InDetDD::SiReadoutCellId SiRC( m_pixelId->phi_index(rdoId), m_pixelId->eta_index(rdoId) );
+ for( int ii=0; ii<nCells && iter==pixelSimDataMap->end(); ++ii ) {
+ iter = pixelSimDataMap->find(sielement->identifierFromCellId(sielement->design().connectedCell(SiRC,ii)));
+ }
+ } // end search for correct ganged pixel
+ // if SDO found for this pixel, associate the particle. otherwise leave unassociated.
+ if( iter!=pixelSimDataMap->end() ) {
+ const InDetSimData& sdo( iter->second );
+ const std::vector<InDetSimData::Deposit>& deposits( sdo.getdeposits() );
+ for( std::vector<InDetSimData::Deposit>::const_iterator iDep=deposits.begin(), fDep=deposits.end(); iDep!=fDep; ++iDep ) {
+ const HepMcParticleLink& particleLink( iDep->first );
+ //const InDetSimData::Deposit::second_type qdep( iDep->second ); // energy(charge) contributed by this particle
+ // RDO's without SDO's are delta rays or detector noise.
+ if( !particleLink.isValid() ) { continue; }
+ const HepMC::GenParticle* particle( particleLink );
+ const float genEta=particle->momentum().pseudoRapidity();
+ const float genPt=particle->momentum().perp(); // MeV
+ // reject unstable particles
+ if( particle->status()%1000!=1 ) { continue; }
+ // reject secondaries and low pT (<400 MeV) pileup
+ if( particle->barcode()>100000 || particle->barcode()==0 ) { continue; }
+ // reject far forward particles
+ if( fabs(genEta)>m_maxEta ) { continue; }
+ // "best_parent" is the highest pt particle
+ if( !best_parent || best_parent->momentum().perp()<genPt ) {
+ best_parent = particle;
+ best_extcode = HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() );
+ }
+ // bcs.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // if( particleLink.eventIndex()==0 ) {
+ // bcs_prim.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // }
+ // if( genPt>=1000. ) {
+ // bcs_pt1gev.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // }
+ // // check whether or not parent is a b
+ // {
+ // typedef pair<const HepMC::GenParticle*,unsigned int> Parent;
+ // vector<Parent> parents;
+ // parents.push_back( Parent(particle,0) );
+ // while( !parents.empty() ) {
+ // const HepMC::GenParticle* p = parents.back().first;
+ // const unsigned int level = parents.back().second;
+ // if( std::abs(p->pdg_id())==5 ) {
+ // bcs_bjet.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // break;
+ // }
+ // parents.pop_back();
+ // if( !(p->production_vertex()) ) { continue; }
+ // for( HepMC::GenVertex::particle_iterator i=p->production_vertex()->particles_begin(HepMC::parents), f=p->production_vertex()->particles_end(HepMC::parents); i!=f; ++i ) {
+ // parents.push_back( Parent(*i,level+1) );
+ // }
+ // }
+ // }
+ parent_mask |= construct_truth_bitmap( particle );
+ // check SDO
+ } // end for each contributing particle
+ } // end if truth found for this pixel
+ } // end if pixel truth available
+ ++hitIndex;
+ (*oflraw) << "S\t"
+ << setw(14) << setprecision(10)
+ << gPos.x() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.y() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.z() << '\t'
+ << 1 << '\t' // 1 pixel 0 sct
+ << m_pixelId->barrel_ec(rdoId) << '\t'
+ << m_pixelId->layer_disk(rdoId) << '\t'
+ << m_pixelId->phi_module(rdoId) << '\t'
+ << m_pixelId->eta_module(rdoId) << '\t'
+ << m_pixelId->phi_index(rdoId) << '\t'
+ << m_pixelId->eta_index(rdoId) << '\t'
+ << (*iRDO)->getToT() << '\t'
+ << (long)(best_parent ? best_extcode.eventIndex() : std::numeric_limits<long>::max()) << '\t'
+ << (long)(best_parent ? best_extcode.barcode() : std::numeric_limits<long>::max()) << '\t'
+ << setprecision(5) << static_cast<unsigned long>(std::ceil(best_parent ? best_parent->momentum().perp() : 0.)) << '\t' // particle pt in MeV
+ << parent_mask.to_ulong() << '\t'
+ << sielement->identifyHash() << '\t'
+ << '\n';
+ } // end for each RDO in the collection
+ } // for each pixel RDO collection
+ // dump all pixel RDO's and SDO's for debugging purposes
+ if(m_doRDODebug) {
+ for( PixelRDO_Container::const_iterator iColl=pixel_rdocontainer_iter->begin(), fColl=pixel_rdocontainer_iter->end(); iColl!=fColl; ++iColl ) {
+ const InDetRawDataCollection<PixelRDORawData>* pixel_rdoCollection(*iColl);
+ if( !pixel_rdoCollection ) { continue; }
+ for( DataVector<PixelRDORawData>::const_iterator iRDO=pixel_rdoCollection->begin(), fRDO=pixel_rdoCollection->end(); iRDO!=fRDO; ++iRDO ) {
+ Identifier rdoId = (*iRDO)->identify();
+ // get the det element from the det element collection
+ const InDetDD::SiDetectorElement* sielement = m_PIX_mgr->getDetectorElement(rdoId); assert( sielement);
+ const InDetDD::SiLocalPosition localPos = sielement->localPositionOfCell(rdoId);
+ const InDetDD::SiLocalPosition rawPos = sielement->rawLocalPositionOfCell(rdoId);
+ const Amg::Vector3D gPos( sielement->globalPosition(localPos) );
+ (*oflraw) << "# S\t"
+ << setw(14) << setprecision(10)
+ << gPos.x() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.y() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.z() << '\t'
+ << 1 << '\t' // 1 pixel 0 sct
+ << m_pixelId->barrel_ec(rdoId) << '\t'
+ << m_pixelId->layer_disk(rdoId) << '\t'
+ << m_pixelId->phi_module(rdoId) << '\t'
+ << m_pixelId->eta_module(rdoId) << '\t'
+ << m_pixelId->phi_index(rdoId) << '\t'
+ << m_pixelId->eta_index(rdoId) << '\t'
+ << (*iRDO)->getToT() << '\t'
+ << endl;
+ } // end for each pixel RDO
+ } // end for each pixel RDO collection
+ // dump SDO's
+ if( have_pixel_sdo && pixelSimDataMap ) {
+ for( InDetSimDataCollection::const_iterator i=pixelSimDataMap->begin(), f=pixelSimDataMap->end(); i!=f; ++i ) {
+ const Identifier sdoId( i->first );
+ const InDetSimData& sdo( i->second );
+ const vector<InDetSimData::Deposit>& deposits( sdo.getdeposits() );
+ (*oflraw) << "# s"
+ << " " << m_pixelId->barrel_ec(sdoId)
+ << " " << m_pixelId->layer_disk(sdoId)
+ << " " << m_pixelId->phi_module(sdoId)
+ << " " << m_pixelId->eta_module(sdoId)
+ << " " << m_pixelId->phi_index(sdoId)
+ << " " << m_pixelId->eta_index(sdoId)
+ << " " << PixelSimHelper::isNoise( sdo )
+ << " " << PixelSimHelper::isBelowThreshold( sdo )
+ << " " << PixelSimHelper::isDisabled( sdo )
+ << " " << PixelSimHelper::hasBadTOT( sdo )
+ << " " << deposits.size()
+ << endl;
+ for( vector<InDetSimData::Deposit>::const_iterator iDep=deposits.begin(), fDep=deposits.end(); iDep!=fDep; ++iDep ) {
+ const HepMcParticleLink& particleLink( iDep->first );
+ const InDetSimData::Deposit::second_type qdep( iDep->second ); // energy(charge) contributed by this particle
+ (*oflraw) << "# s q " << qdep << " " << particleLink.isValid() << " "
+ << (particleLink.isValid() ? particleLink.eventIndex() : std::numeric_limits<unsigned int>::max())
+ << (particleLink.isValid() ? particleLink.barcode() : std::numeric_limits<unsigned int>::max())
+ << endl;
+ } // end for each deposit in SDO
+ } // end for each pixel SDO
+ } // if have sdo info
+ } // end dump RDO's and SDO's for debugging purposes (m_doRDODebug)
+ } // dump raw pixel data
+
+ const InDetSimDataCollection* sctSimDataMap(0);
+ const bool have_sct_sdo = m_storeGate->retrieve(sctSimDataMap, "SCT_SDO_Map").isSuccess();
+ const DataHandle<SCT_RDO_Container> sct_rdocontainer_iter;
+ if( m_storeGate->retrieve(sct_rdocontainer_iter, "SCT_RDOs").isSuccess() ) {
+ sct_rdocontainer_iter->clID(); // anything to dereference the DataHandle
+ for( SCT_RDO_Container::const_iterator iColl=sct_rdocontainer_iter->begin(), fColl=sct_rdocontainer_iter->end(); iColl!=fColl; ++iColl ) {
+ const InDetRawDataCollection<SCT_RDORawData>* SCT_Collection(*iColl);
+ if( !SCT_Collection ) { continue; }
+ const int size = SCT_Collection->size();
+ ATH_MSG_DEBUG( "SCT InDetRawDataCollection found with " << size << " RDOs");
+ for( DataVector<SCT_RDORawData>::const_iterator iRDO=SCT_Collection->begin(), fRDO=SCT_Collection->end(); iRDO!=fRDO; ++iRDO ) {
+ const Identifier rdoId = (*iRDO)->identify();
+ // get the det element from the det element collection
+ const InDetDD::SiDetectorElement* sielement = m_SCT_mgr->getDetectorElement(rdoId);
+ const InDetDD::SCT_ModuleSideDesign& design = dynamic_cast<const InDetDD::SCT_ModuleSideDesign&>(sielement->design());
+ const InDetDD::SiLocalPosition localPos = design.positionFromStrip((*iRDO)->getStrip());
+ const Amg::Vector3D gPos = sielement->globalPosition(localPos);
+ hitIndexMap[rdoId] = hitIndex;
+ ++hitIndex;
+ // if there is simulation truth available, try to retrieve the
+ // "most likely" barcode for this strip.
+ const HepMC::GenParticle* best_parent = 0;
+ ParentBitmask parent_mask;
+ HepMcParticleLink::ExtendedBarCode best_extcode;
+ if( have_sct_sdo && sctSimDataMap ) {
+ InDetSimDataCollection::const_iterator iter( sctSimDataMap->find(rdoId) );
+ // if SDO found for this pixel, associate the particle
+ if( iter!=sctSimDataMap->end() ) {
+ const InDetSimData& sdo( iter->second );
+ const std::vector<InDetSimData::Deposit>& deposits( sdo.getdeposits() );
+ for( std::vector<InDetSimData::Deposit>::const_iterator iDep=deposits.begin(), fDep=deposits.end(); iDep!=fDep; ++iDep ) {
+ const HepMcParticleLink& particleLink( iDep->first );
+ // const InDetSimData::Deposit::second_type qdep( iDep->second ); // energy(charge) contributed by this particle
+ // RDO's without SDO's are delta rays or detector noise.
+ if( !particleLink.isValid() ) { continue; }
+ const HepMC::GenParticle* particle( particleLink );
+ const float genEta=particle->momentum().pseudoRapidity();
+ const float genPt=particle->momentum().perp(); // MeV
+ // reject unstable particles
+ if( particle->status()%1000!=1 ) { continue; }
+ // reject secondaries and low pt (<400 MeV) pileup truth
+ if( particle->barcode()>100000 || particle->barcode()==0 ) { continue; }
+ // reject far forward particles
+ if( fabs(genEta)>m_maxEta ) { continue; }
+ // "best_parent" is the highest pt particle
+ if( !best_parent || best_parent->momentum().perp()<genPt ) {
+ best_parent = particle;
+ best_extcode = HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() );
+ }
+ // bcs.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // if( particleLink.eventIndex()==0 ) {
+ // bcs_prim.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // }
+ // if( genPt>=1000. ) {
+ // bcs_pt1gev.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // }
+ // {
+ // typedef pair<const HepMC::GenParticle*,unsigned int> Parent;
+ // vector<Parent> parents;
+ // parents.push_back( Parent(particle,0) );
+ // while( !parents.empty() ) {
+ // const HepMC::GenParticle* p = parents.back().first;
+ // const unsigned int level = parents.back().second;
+ // if( std::abs(p->pdg_id())==5 ) {
+ // bcs_bjet.insert( HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() ) );
+ // break;
+ // }
+ // parents.pop_back();
+ // if( !(p->production_vertex()) ) { continue; }
+ // for( HepMC::GenVertex::particle_iterator i=p->production_vertex()->particles_begin(HepMC::parents), f=p->production_vertex()->particles_end(HepMC::parents); i!=f; ++i ) {
+ // parents.push_back( Parent(*i,level+1) );
+ // }
+ // }
+ // }
+ parent_mask |= construct_truth_bitmap( particle );
+ } // end for each contributing particle
+ } // end if truth found for this strip
+ } // end if sct truth available
+ (*oflraw) << "S\t"
+ << setw(14) << setprecision(10)
+ << gPos.x() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.y() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.z() << '\t'
+ << 0 << '\t' // 1 pixel 0 sct
+ << m_sctId->barrel_ec(rdoId) << '\t'
+ << m_sctId->layer_disk(rdoId) << '\t'
+ << m_sctId->phi_module(rdoId) << '\t'
+ << m_sctId->eta_module(rdoId) << '\t'
+ << m_sctId->side(rdoId) << '\t'
+ << m_sctId->strip(rdoId) << '\t'
+ << (*iRDO)->getGroupSize() << '\t'
+ << (long)(best_parent ? best_extcode.eventIndex() : std::numeric_limits<long>::max()) << '\t'
+ << (long)(best_parent ? best_extcode.barcode() : std::numeric_limits<long>::max()) << '\t'
+ << setprecision(5) << static_cast<unsigned long>(std::ceil(best_parent ? best_parent->momentum().perp() : 0.)) << '\t' // particle pt in MeV
+ << parent_mask.to_ulong() << '\t'
+ << sielement->identifyHash() << '\t'
+ << '\n';
+ } // end for each RDO in the strip collection
+ } // end for each strip RDO collection
+ // dump all RDO's and SDO's for a given event, for debugging purposes
+ if(m_doRDODebug) {
+ // dump SCT RDO's
+ for( SCT_RDO_Container::const_iterator iColl=sct_rdocontainer_iter->begin(), fColl=sct_rdocontainer_iter->end(); iColl!=fColl; ++iColl ) {
+ const InDetRawDataCollection<SCT_RDORawData>* SCT_Collection(*iColl);
+ if( !SCT_Collection ) { continue; }
+ const int size = SCT_Collection->size();
+ ATH_MSG_DEBUG( "SCT InDetRawDataCollection found with " << size << " RDOs");
+ for( DataVector<SCT_RDORawData>::const_iterator iRDO=SCT_Collection->begin(), fRDO=SCT_Collection->end(); iRDO!=fRDO; ++iRDO ) {
+ const Identifier rdoId = (*iRDO)->identify();
+ const InDetDD::SiDetectorElement* sielement = m_SCT_mgr->getDetectorElement(rdoId);
+ const InDetDD::SCT_ModuleSideDesign& design = dynamic_cast<const InDetDD::SCT_ModuleSideDesign&>(sielement->design());
+ const InDetDD::SiLocalPosition localPos = design.positionFromStrip((*iRDO)->getStrip());
+ const Amg::Vector3D gPos = sielement->globalPosition(localPos);
+ (*oflraw) << "# S\t"
+ << setw(14) << setprecision(10)
+ << gPos.x() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.y() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.z() << '\t'
+ << 0 << '\t' // 1 pixel 0 sct
+ << m_sctId->barrel_ec(rdoId) << '\t'
+ << m_sctId->layer_disk(rdoId) << '\t'
+ << m_sctId->phi_module(rdoId) << '\t'
+ << m_sctId->eta_module(rdoId) << '\t'
+ << m_sctId->side(rdoId) << '\t'
+ << m_sctId->strip(rdoId) << '\t'
+ << (*iRDO)->getGroupSize() << '\t'
+ << endl;
+ } // end for each SCT rdo
+ } // end for each SCT rdo collection
+ // dump SCT SDO's
+ if( have_sct_sdo && sctSimDataMap ) {
+ for( InDetSimDataCollection::const_iterator i=sctSimDataMap->begin(), f=sctSimDataMap->end(); i!=f; ++i ) {
+ const Identifier sdoId( i->first );
+ const InDetSimData& sdo( i->second );
+ const vector<InDetSimData::Deposit>& deposits( sdo.getdeposits() );
+ (*oflraw) << "# s"
+ << " " << m_sctId->barrel_ec(sdoId)
+ << " " << m_sctId->layer_disk(sdoId)
+ << " " << m_sctId->phi_module(sdoId)
+ << " " << m_sctId->eta_module(sdoId)
+ << " " << m_sctId->side(sdoId)
+ << " " << m_sctId->strip(sdoId)
+ << " " << SCT_SimHelper::isNoise( sdo )
+ << " " << SCT_SimHelper::isBelowThreshold( sdo )
+ << " " << SCT_SimHelper::isDisabled( sdo )
+ << " " << deposits.size()
+ << endl;
+ for( vector<InDetSimData::Deposit>::const_iterator iDep=deposits.begin(), fDep=deposits.end(); iDep!=fDep; ++iDep ) {
+ const HepMcParticleLink& particleLink( iDep->first );
+ const InDetSimData::Deposit::second_type qdep( iDep->second ); // energy(charge) contributed by this particle
+ (*oflraw) << "# s q " << qdep << " " << particleLink.isValid() << " "
+ << (particleLink.isValid() ? particleLink.eventIndex() : std::numeric_limits<unsigned int>::max())
+ << (particleLink.isValid() ? particleLink.barcode() : std::numeric_limits<unsigned int>::max())
+ << endl;
+ } // end for each deposit associated with this SDO
+ } // end for each SCT SDO
+ } // end if SDO's available, dump them
+ } // end dump all RDO's and SDO's for a given event (m_doRDODebug)
+ } // end dump raw SCT data
+
+ // FlagJT dump pixel clusters. They're in m_pixelContainer
+ m_pixelContainer->clID(); // anything to dereference the DataHandle
+ for( InDet::SiClusterContainer::const_iterator iColl=m_pixelContainer->begin(), fColl=m_pixelContainer->end(); iColl!=fColl; ++iColl ) {
+ const InDet::SiClusterCollection* pixelClusterCollection( *iColl );
+ if( !pixelClusterCollection ) {
+ ATH_MSG_DEBUG( "pixelClusterCollection not available!");
+ continue;
+ }
+ const unsigned int size = pixelClusterCollection->size();
+ ATH_MSG_DEBUG( "PixelClusterCollection found with " << size << " clusters");
+ for( DataVector<InDet::SiCluster>::const_iterator iCluster=pixelClusterCollection->begin(), fCluster=pixelClusterCollection->end(); iCluster!=fCluster; ++iCluster ) {
+ Identifier theId = (*iCluster)->identify();
+ const Amg::Vector3D gPos = (*iCluster)->globalPosition();
+ // if there is simulation truth available, try to retrieve the "most likely" barcode for this pixel cluster.
+ const HepMC::GenParticle* best_parent = 0;
+ ParentBitmask parent_mask;
+ HepMcParticleLink::ExtendedBarCode best_extcode;
+ if( have_pixel_sdo && pixelSimDataMap ) {
+ for( std::vector<Identifier>::const_iterator rdoIter = (*iCluster)->rdoList().begin();
+ rdoIter != (*iCluster)->rdoList().end(); rdoIter++ ) {
+ const InDetDD::SiDetectorElement* sielement = m_PIX_mgr->getDetectorElement(*rdoIter);
+ assert( sielement );
+ const InDetDD::SiLocalPosition rawPos = sielement->rawLocalPositionOfCell(*rdoIter);
+ const int nCells = sielement->numberOfConnectedCells( sielement->cellIdOfPosition(rawPos) );
+ InDetSimDataCollection::const_iterator iter( pixelSimDataMap->find(*rdoIter) );
+ // this might be the ganged pixel copy.
+ if( nCells>1 && iter==pixelSimDataMap->end() ) {
+ InDetDD::SiReadoutCellId SiRC( m_pixelId->phi_index(*rdoIter), m_pixelId->eta_index(*rdoIter) );
+ for( int ii=0; ii<nCells && iter==pixelSimDataMap->end(); ++ii ) {
+ iter = pixelSimDataMap->find(sielement->identifierFromCellId(sielement->design().connectedCell(SiRC,ii)));
+ }
+ } // end search for correct ganged pixel
+ // if SDO found for this pixel, associate the particle. otherwise leave unassociated.
+ if( iter!=pixelSimDataMap->end() ) {
+ const InDetSimData& sdo( iter->second );
+ const std::vector<InDetSimData::Deposit>& deposits( sdo.getdeposits() );
+ for( std::vector<InDetSimData::Deposit>::const_iterator iDep=deposits.begin(), fDep=deposits.end(); iDep!=fDep; ++iDep ) {
+ const HepMcParticleLink& particleLink( iDep->first );
+ //const InDetSimData::Deposit::second_type qdep( iDep->second ); // energy(charge) contributed by this particle
+ // RDO's without SDO's are delta rays or detector noise.
+ if( !particleLink.isValid() ) { continue; }
+ const HepMC::GenParticle* particle( particleLink );
+ const float genEta=particle->momentum().pseudoRapidity();
+ const float genPt=particle->momentum().perp(); // MeV
+ // reject unstable particles
+ if( particle->status()%1000!=1 ) { continue; }
+ // reject secondaries and low pT (<400 MeV) pileup
+ if( particle->barcode()>100000 || particle->barcode()==0 ) { continue; }
+ // reject far forward particles
+ if( fabs(genEta)>m_maxEta ) { continue; }
+ // "best_parent" is the highest pt particle
+ if( !best_parent || best_parent->momentum().perp()<genPt ) {
+ best_parent = particle;
+ best_extcode = HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() );
+ }
+ parent_mask |= construct_truth_bitmap( particle );
+ } // loop over deposits
+ } // if found SDO of pixel
+ } // loop over pixels in cluster
+ } // if we have pixel sdo's available
+
+ const InDetDD::SiDetectorElement* sielement( (*iCluster)->detectorElement() ); assert( sielement );
+
+ // compute hit position in units of channel number. local hit
+ // position is measured in a right-handed coordinate system
+ // centered on the center of a wafer
+ // (SiDetectorElement::center() in global coordinates) with axes
+ // SiDetectorElement::phiAxis() and SiDetectorElement::etaAxis()
+ // assume all rdos in cluster are on the same physical
+ // wafer. note that pixel pitch is not constant vs channel
+ // number, but the Trk::Surface models of disks and barrel
+ // (cylinder) treat it as such.
+
+ float localx=-1,localy=-1;
+ {
+ // const InDetDD::SiLocalPosition localCentroid( (*iCluster)->localPosition() );
+ // Amg::Vector2D corrLocalCentroid(localCentroid);
+ //
+ // corrLocalCentroid[Trk::distPhi] -= sielement->getLorentzCorrection();
+ //
+ // const InDetDD::SiCellId cellIdCentroid( sielement->cellIdOfPosition(corrLocalCentroid ));//localCentroid ) );
+ // //const Amg::Vector2D localCellCentroid( sielement->localPositionOfCell( cellIdCentroid ) );
+ //
+ //
+ // //using the non Lorentz shift-corrected value.
+ // const Amg::Vector2D localCellCentroid( sielement->localPositionOfCell( cellIdCentroid ) );
+ //
+ // float deltaphi = (*iCluster)->localPosition()[Trk::distPhi] - localCellCentroid[Trk::distPhi];
+ //
+ //
+ // //using the non Lorentz shift-corrected value.
+ // const Trk::LocalPosition localCellCentroid( sielement->localPositionOfCell( cellIdCentroid ) );
+ //
+ // float deltaphi = (*iCluster)->localPosition()[Trk::distPhi] - localCellCentroid[Trk::distPhi];
+ //
+ // float deltaeta = (*iCluster)->localPosition()[Trk::distEta] - localCellCentroid[Trk::distEta];
+ // float deltaxphi = deltaphi / sielement->phiPitch();
+ // float deltaxeta = deltaeta / sielement->etaPitch();
+ // // do not add 0.5 to center in middle of strip; match SCT leading spatial edge convention
+ // localx = cellIdCentroid.phiIndex() + deltaxphi; // + 0.5
+ // localy = cellIdCentroid.etaIndex() + deltaxeta; // + 0.5
+
+ // FlagAA: THIS CODE IS NOT BACKWARD COMPATIBLE! Commited on Sep 18th, 2013
+ // now stores the local position in millimiters without Lorentz Correction
+ // before it was storing a floating point coordinate in units of pixels
+ localx = (*iCluster)->localPosition()[Trk::distPhi] - sielement->getLorentzCorrection();;
+ localy = (*iCluster)->localPosition()[Trk::distEta];
+ }
+
+ (*oflraw) << "P\t"
+ << setw(14) << setprecision(10)
+ << gPos.x() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.y() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.z() << '\t'
+ << m_pixelId->barrel_ec(theId) << '\t'
+ << m_pixelId->layer_disk(theId) << '\t'
+ << m_pixelId->phi_module(theId) << '\t'
+ << m_pixelId->eta_module(theId) << '\t'
+ << m_pixelId->phi_index(theId) << '\t'
+ << m_pixelId->eta_index(theId) << '\t'
+ // Cluster local position
+ << localx << '\t'
+ << localy << '\t'
+ // Cluster width:
+ << (*iCluster)->width().colRow().x() << '\t' // width in phi?
+ << (*iCluster)->width().colRow().y() << '\t' // width in eta?
+ << (*iCluster)->rdoList().size() << '\t' // number of pixels in cluster
+ // Cluster truth
+ << (long)(best_parent ? best_extcode.eventIndex() : std::numeric_limits<long>::max()) << '\t'
+ << (long)(best_parent ? best_extcode.barcode() : std::numeric_limits<long>::max()) << '\t'
+ << setprecision(5) << static_cast<unsigned long>(std::ceil(best_parent ? best_parent->momentum().perp() : 0.)) << '\t' // particle pt in MeV
+ << parent_mask.to_ulong()
+ // Module identifier hash
+ << sielement->identifyHash()
+ //
+ << endl;
+ clusterIndexMap[theId] = clusterIndex++;
+ } // End loop over pixel clusters
+ } // End loop over pixel cluster collection
+
+ // dump SCT merged clusters.
+ m_sctContainer->clID(); // anything to dereference the DataHandle
+ for( InDet::SiClusterContainer::const_iterator iColl=m_sctContainer->begin(), fColl=m_sctContainer->end(); iColl!=fColl; ++iColl ) {
+ const InDet::SiClusterCollection* sctClusterCollection( *iColl );
+ if( !sctClusterCollection ) {
+ ATH_MSG_DEBUG( "sctClusterCollection not available!");
+ continue;
+ }
+ const unsigned int size = sctClusterCollection->size();
+ ATH_MSG_DEBUG( "SctClusterCollection found with " << size << " clusters");
+ for( DataVector<InDet::SiCluster>::const_iterator iCluster=sctClusterCollection->begin(), fCluster=sctClusterCollection->end(); iCluster!=fCluster; ++iCluster ) {
+ Identifier theId = (*iCluster)->identify();
+ const Amg::Vector3D gPos = (*iCluster)->globalPosition();
+ // if there is simulation truth available, try to retrieve the "most likely" barcode for this sct cluster.
+ const HepMC::GenParticle* best_parent = 0;
+ ParentBitmask parent_mask;
+ HepMcParticleLink::ExtendedBarCode best_extcode;
+ if( have_sct_sdo && sctSimDataMap ) {
+ for( std::vector<Identifier>::const_iterator rdoIter = (*iCluster)->rdoList().begin();
+ rdoIter != (*iCluster)->rdoList().end(); rdoIter++ ) {
+ const InDetDD::SiDetectorElement* sielement = m_SCT_mgr->getDetectorElement(*rdoIter);
+ assert( sielement );
+ const InDetDD::SiLocalPosition rawPos = sielement->rawLocalPositionOfCell(*rdoIter);
+ InDetSimDataCollection::const_iterator iter( sctSimDataMap->find(*rdoIter) );
+ // if SDO found for this cluster, associate the particle. otherwise leave unassociated.
+ if( iter!=sctSimDataMap->end() ) {
+ const InDetSimData& sdo( iter->second );
+ const std::vector<InDetSimData::Deposit>& deposits( sdo.getdeposits() );
+ for( std::vector<InDetSimData::Deposit>::const_iterator iDep=deposits.begin(), fDep=deposits.end(); iDep!=fDep; ++iDep ) {
+ const HepMcParticleLink& particleLink( iDep->first );
+ //const InDetSimData::Deposit::second_type qdep( iDep->second ); // energy(charge) contributed by this particle
+ // RDO's without SDO's are delta rays or detector noise.
+ if( !particleLink.isValid() ) { continue; }
+ const HepMC::GenParticle* particle( particleLink );
+ const float genEta=particle->momentum().pseudoRapidity();
+ const float genPt=particle->momentum().perp(); // MeV
+ // reject unstable particles
+ if( particle->status()%1000!=1 ) { continue; }
+ // reject secondaries and low pT (<400 MeV) pileup
+ if( particle->barcode()>100000 || particle->barcode()==0 ) { continue; }
+ // reject far forward particles
+ if( fabs(genEta)>m_maxEta ) { continue; }
+ // "best_parent" is the highest pt particle
+ if( !best_parent || best_parent->momentum().perp()<genPt ) {
+ best_parent = particle;
+ best_extcode = HepMcParticleLink::ExtendedBarCode( particleLink.barcode() , particleLink.eventIndex() );
+ }
+ parent_mask |= construct_truth_bitmap( particle );
+ } // loop over deposits
+ } // if found SDO of sct
+ } // loop over scts in cluster
+ } // if we have sct sdo's available
+
+ const InDetDD::SiDetectorElement* sielement( (*iCluster)->detectorElement() ); assert( sielement );
+
+ // compute hit position in units of strip number. local hit
+ // position is measured in a right-handed coordinate system
+ // centered on the center of a wafer
+ // (SiDetectorElement::center() in global coordinates) with axes
+ // SiDetectorElement::phiAxis() and SiDetectorElement::etaAxis()
+ // assume all rdos in cluster are on the same physical
+ // wafer. note that pitch is not necessarily a constant vs channel
+ // number, but the Trk::Surface models of disks and barrel
+ // (cylinder) treat it as such.
+
+ float localx=-1;
+ {
+ const InDetDD::SiLocalPosition localCentroid( (*iCluster)->localPosition() );
+ const InDetDD::SiCellId cellIdCentroid( sielement->cellIdOfPosition( localCentroid ) );
+ const Amg::Vector2D localCellCentroid( sielement->localPositionOfCell( cellIdCentroid ) );
+ float deltaphi = (*iCluster)->localPosition()[Trk::distPhi] - localCellCentroid[Trk::distPhi];
+ float deltaxphi = deltaphi / sielement->phiPitch();
+ // do not add 0.5 to center in middle of strip; match SCT leading spatial edge convention
+ localx = cellIdCentroid.phiIndex() + deltaxphi; // + 0.5
+ }
+
+ (*oflraw) << "C\t"
+ << setw(14) << setprecision(10)
+ << gPos.x() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.y() << '\t'
+ << setw(14) << setprecision(10)
+ << gPos.z() << '\t'
+ << 0 << '\t' // 1 pixel 0 sct
+ << m_sctId->barrel_ec(theId) << '\t'
+ << m_sctId->layer_disk(theId) << '\t'
+ << m_sctId->phi_module(theId) << '\t'
+ << m_sctId->eta_module(theId) << '\t'
+ << m_sctId->side(theId) << '\t'
+ // Cluster local position and width
+ << m_sctId->strip(theId) << '\t'
+ << localx << '\t'
+ << (*iCluster)->rdoList().size() << '\t' // number of strips in cluster
+ << (long)(best_parent ? best_extcode.eventIndex() : std::numeric_limits<long>::max()) << '\t'
+ << (long)(best_parent ? best_extcode.barcode() : std::numeric_limits<long>::max()) << '\t'
+ << setprecision(5) << static_cast<unsigned long>(std::ceil(best_parent ? best_parent->momentum().perp() : 0.)) << '\t' // particle pt in MeV
+ << parent_mask.to_ulong() << '\t'
+ << sielement->identifyHash() << '\t'
+ << endl;
+ clusterIndexMap[theId] = clusterIndex++;
+ } // End loop over sct clusters
+ } // End loop over sct cluster collection
+
+ // dump the statistics
+ // cout << boost::format("truth parent stats: nch = %|6d| %|6d| %|6d| %|6d| %|6d| fmt = %|5g| fmtp = %|5g| fmtp1 = %|5g| fmb = %|5g|")
+ // % nchannels
+ // % nchannels_single_barcodes
+ // % nchannels_single_barcodes_primary
+ // % nchannels_single_barcodes_pt1gev
+ // % nchannels_single_barcodes_bjet
+ // % (nchannels>0ul ? nchannels_multiple_barcodes/static_cast<double>(nchannels) : 0.)
+ // % (nchannels>0ul ? nchannels_multiple_barcodes_primary/static_cast<double>(nchannels) : 0.)
+ // % (nchannels>0ul ? nchannels_multiple_barcodes_pt1gev/static_cast<double>(nchannels) : 0.)
+ // % (nchannels>0ul ? nchannels_multiple_barcodes_bjet/static_cast<double>(nchannels) : 0.)
+ // << endl;
+
+} // end dump raw silicon data
+
+
+// dump list of bad pixel and SCT modules for this event
+
+void
+DumpSp::dump_bad_modules() const
+{
+ if( m_doBadMod) {
+ // dump list of bad pixel modules
+ for( InDetDD::SiDetectorElementCollection::const_iterator i=m_PIX_mgr->getDetectorElementBegin(), f=m_PIX_mgr->getDetectorElementEnd(); i!=f; ++i ) {
+ const InDetDD::SiDetectorElement* sielement( *i );
+ Identifier id = sielement->identify();
+ IdentifierHash idhash = sielement->identifyHash();
+ const bool is_bad = !(m_pixelCondSummarySvc->isGood( idhash ));
+ if( is_bad ) {
+ (*oflraw) << "B\t"
+ << 1 << '\t' // 1 pixel 0 sct
+ << m_pixelId->barrel_ec(id) << '\t'
+ << m_pixelId->layer_disk(id) << '\t'
+ << m_pixelId->phi_module(id) << '\t'
+ << m_pixelId->eta_module(id) << '\t'
+ // << m_pixelId->phi_index(id) << '\t'
+ // << m_pixelId->eta_index(id) << '\t'
+ // << is_bad
+ << std::endl;
+ }
+ } // end for each pixel module
+ // dump list of bad sct modules
+ for( InDetDD::SiDetectorElementCollection::const_iterator i=m_SCT_mgr->getDetectorElementBegin(), f=m_SCT_mgr->getDetectorElementEnd(); i!=f; ++i ) {
+ const InDetDD::SiDetectorElement* sielement( *i );
+ Identifier id = sielement->identify();
+ IdentifierHash idhash = sielement->identifyHash();
+ const bool is_bad = !(m_sctCondSummarySvc->isGood( idhash ));
+ if( is_bad ) {
+ (*oflraw) << "B\t"
+ << 0 << '\t' // 1 pixel 0 sct
+ << m_sctId->barrel_ec(id) << '\t'
+ << m_sctId->layer_disk(id) << '\t'
+ << m_sctId->phi_module(id) << '\t'
+ << m_sctId->eta_module(id) << '\t'
+ << m_sctId->side(id) << '\t'
+ // << is_bad
+ << std::endl;
+ }
+ } // end for each SCT module
+ } // end dump bad silicon modules
+} // end dump list of bad silicon modules
+
+
+
+// dump reconstructed tracks with geant matching information.
+
+void
+DumpSp::dump_tracks( const HitIndexMap& /*hitIndexMap*/, const HitIndexMap& clusterIndexMap ) const
+{
+ ATH_MSG_DEBUG( "getting the reconstructed track container.");
+ const DataVector<Trk::Track>* trks = 0;
+ if( m_storeGate->retrieve(trks,m_tracksName).isFailure() ) {
+ ATH_MSG_WARNING( "unable to retrieve reconstructed TrackCollection " << m_tracksName);
+ return;
+ }
+
+ if( !(trks->empty()) ) {
+ ATH_MSG_DEBUG( "found reconstructed track container with " << trks->size() << " tracks.");
+ ATH_MSG_DEBUG( "begin dumping reconstructed tracks");
+ }
+
+ unsigned int itrack = 0u;
+ for( DataVector<Trk::Track>::const_iterator iTrk=trks->begin(), iTrkF=trks->end(); iTrk!=iTrkF; ++iTrk, ++itrack ) {
+ const Trk::Track* track( *iTrk );
+ const Amg::Vector3D vertex(0.,0.,0.); // FIXME --- retrieve the fitted primary vertex position here
+ const double charge = track->perigeeParameters()->charge();
+ const Amg::Vector3D perigee_position( track->perigeeParameters()->position() );
+ const double z0 = perigee_position.z();
+ const double phi0 = track->perigeeParameters()->parameters()[Trk::phi];
+ const double qpt = charge*std::abs(((*track).perigeeParameters()->pT())/CLHEP::GeV);
+ // FIXME --- verify that this the correct d0 sign convention in
+ // the presence of beamline displacement. low priority until
+ // vertex!=0.
+ const Amg::Vector3D beam_to_perigee_vector( perigee_position - vertex );
+ const double d0_sign = std::sin( beam_to_perigee_vector.phi() - phi0 ) > 0. ? 1. : -1.;
+ const double d0 = d0_sign*( sqrt( pow(perigee_position.x()-vertex.x(),(int)2) + pow(perigee_position.y()-vertex.y(),(int)2) ) - vertex.perp() );
+ // apply track selection if requested
+ if( m_useOfflineTrackSelectorTool ) {
+ // v13 does not have an IDetailedTrackSelectorTool, so make the cuts manually here. use
+ // exactly the cuts configured by InDetRecExample, except 1) the track parameter extrapolation
+ // is the perigee extrapolation provided by the track, not the (potentially) different calculation
+ // employing a custom extrapolator, and 2) the cuts on the number of silicon hits are not corrected
+ // for dead SCT or pixel layers; this information is not available in v13 TrackSummary.
+ //
+ // get track parameters using the same extrapolation as DetailedTrackSelectorTool
+ // const Trk::Perigee* perigeeBeforeExtrapolation = dynamic_cast<const Trk::Perigee*>( track->perigeeParameters() );
+ const Trk::Perigee* extrapolatedPerigee = dynamic_cast<const Trk::Perigee*>( track->perigeeParameters() );
+ if( false ) {
+ // recompute the extrapolated perigee using a new extrapolator. the actual
+ // extrapolation is commented out, because this algorithm does not have
+ // an extrapolator configured.
+ Trk::PerigeeSurface perigeeSurface( vertex );
+ const Trk::TrackParameters* firstmeaspar = 0;
+ for( unsigned int i=0, f=track->trackParameters()->size(); i!=f; ++i ) {
+ if( dynamic_cast<const Trk::TrackParameters*>((*track->trackParameters())[i]) &&
+ !dynamic_cast<const Trk::Perigee*>((*track->trackParameters())[i]) ) {
+ firstmeaspar=(*track->trackParameters())[i];
+ break;
+ }
+ }
+ if( !firstmeaspar ) {
+ firstmeaspar = track->perigeeParameters();
+ if( !firstmeaspar ) { continue; } // no track selection if firstmeas + perigee does not exist !
+ }
+ //Trk::PropDirection dir( Trk::oppositeMomentum );
+ // for now, approximate the following statements:
+ // boost::scoped_ptr<const Trk::ParametersBase> extrapolatedParameters( firstmeaspar ?
+ // m_extrapolator->extrapolate(*firstmeaspar,perigeeSurface,dir,true,track->particleHypothesis() ) :
+ // 0 );
+ // extrapolatedPerigee = extrapolatedParameters ? dynamic_cast<const Trk::MeasuredPerigee*>(extrapolatedParameters.get()) : 0;
+ // with the default extrapolated parameters coming from the track.
+ }
+ if( !extrapolatedPerigee ) { continue; } // no track selection if firstmeas + perigee does not exist !
+ const AmgVector(5) &perigeeParms = extrapolatedPerigee->parameters();
+ if( std::abs(perigeeParms[Trk::d0]) > 1. ) { continue; }
+ if( std::abs(perigeeParms[Trk::z0]) > 1000. ) { continue; }
+ if( Amg::error(*(extrapolatedPerigee->covariance()),Trk::d0)>0.35 ) { continue; }
+ if( Amg::error(*(extrapolatedPerigee->covariance()),Trk::z0)>2.5 ) { continue; }
+ const Trk::FitQuality* trkQuality = track->fitQuality();
+ const double chi2 = trkQuality->chiSquared();
+ const int ndf = trkQuality->numberDoF();
+ if( chi2/double(ndf) > 3.5 ) { continue; }
+ boost::scoped_ptr<const Trk::TrackSummary> summary( m_trkSummaryTool->createSummary( *track ) );
+ const unsigned int nsct = summary->get(Trk::numberOfSCTHits);
+ const unsigned int npix = summary->get(Trk::numberOfPixelHits);
+ // const unsigned int nsct = (summary->get(Trk::numberOfSCTHits)+summary->get(Trk::numberOfSCTDeadSensors));
+ // const unsigned int npix = (summary->get(Trk::numberOfPixelHits) + summary->get(Trk::numberOfPixelDeadSensors));
+ if( nsct < 5 ) { continue; }
+ if( npix < 1 ) { continue; }
+ if( nsct+npix < 7 ) { continue; }
+ if( summary->get(Trk::numberOfSCTSharedHits)>1 ) { continue; }
+ // compute d0 and z0 significance.
+ // const double sinTheta = sin(perigeeParms[Trk::theta]);
+ //const double cosTheta = cos(perigeeParms[Trk::theta]);
+ const double d0wrtPriVtx = perigeeParms[Trk::d0];
+ //const double deltaZ = perigeeParms[Trk::z0];
+ //const double z0wrtPriVtx = deltaZ*sinTheta;
+ const double testtrackSigD0 = Amg::error(*(extrapolatedPerigee->covariance()),Trk::d0);
+ //const double testtrackSigZ0 = extrapolatedPerigee->localErrorMatrix().error(Trk::z0);
+ //const double testtrackSigTh = extrapolatedPerigee->localErrorMatrix().error(Trk::theta);
+ // error on IP:
+ //const double trackPhi = perigeeParms[Trk::phi];
+ //const double dIPdx = sin(trackPhi);
+ //const double dIPdy = -cos(trackPhi);
+ const double DD0 = testtrackSigD0*testtrackSigD0;
+ double newD0Err=0;
+ // const double DXX = dIPdx*dIPdx*vertex.errorPosition().covariance()[0][0];
+ // const double DYY = dIPdy*dIPdy*vertex.errorPosition().covariance()[1][1];
+ // const double DXY = 2.*dIPdx*dIPdy*vertex.errorPosition().covariance()[0][1];
+ // newD0Err = DD0 + DXX + DYY + DXY;
+ // }
+ newD0Err = DD0; // no vertex to correct
+ const double d0ErrwrtPriVtx = (newD0Err>0 ? sqrt(newD0Err) : -10e-9);
+ if( std::abs(d0wrtPriVtx/d0ErrwrtPriVtx)>4. ) { continue; }
+ }
+ // retrieve any truth info for this track.
+ RecoToTruthTrackMap::const_iterator irtt = _rttTrackMap.find( itrack );
+ HepMcParticleLink::ExtendedBarCode best_extcode;
+ float mc_frac = -0.001;
+ if( irtt!=_rttTrackMap.end() ) {
+ best_extcode = irtt->second;
+ TruthToRecoProbMap::const_iterator ittr = _ttrProbMap.find( irtt->second );
+ assert( ittr != _ttrProbMap.end() );
+ mc_frac = ittr->second;
+ }
+ // dump one line for each track
+ if( m_dumpSpacePoints ) {
+ (*ofl) << "E\t"
+ << setw(14) << setprecision(10) << d0 << '\t'
+ << setw(14) << setprecision(10) << z0 << '\t'
+ << setw(14) << setprecision(10) << phi0 << '\t'
+ << setw(14) << setprecision(10) << 1./tan(track->perigeeParameters()->parameters()[Trk::theta]) << '\t'
+ << setw(15) << setprecision(10) << qpt << '\t'
+ << setw(14) << (long)(mc_frac>=0. ? best_extcode.eventIndex() : std::numeric_limits<long>::max()) << '\t'
+ << setw(14) << (long)(mc_frac>=0. ? best_extcode.barcode() : std::numeric_limits<long>::max()) << '\t'
+ << setw(14) << setprecision(10) << mc_frac
+ << endl;
+ }
+ (*oflraw) << "E\t"
+ << setw(14) << setprecision(10) << d0 << '\t'
+ << setw(14) << setprecision(10) << z0 << '\t'
+ << setw(14) << setprecision(10) << phi0 << '\t'
+ << setw(14) << setprecision(10) << 1./tan(track->perigeeParameters()->parameters()[Trk::theta]) << '\t'
+ << setw(15) << setprecision(10) << qpt << '\t'
+ << setw(14) << (mc_frac>=0. ? best_extcode.eventIndex() : std::numeric_limits<long>::max()) << '\t'
+ << setw(14) << (mc_frac>=0. ? best_extcode.barcode() : std::numeric_limits<long>::max()) << '\t'
+ << setw(14) << setprecision(10) << mc_frac << '\t';
+ // dump the indices of the clusters used in this track
+ for( DataVector<const Trk::MeasurementBase>::const_iterator hitOnTrack = track->measurementsOnTrack()->begin();
+ hitOnTrack != track->measurementsOnTrack()->end(); hitOnTrack++ ) {
+ const Trk::RIO_OnTrack* hit = dynamic_cast<const Trk::RIO_OnTrack*>(*hitOnTrack);
+ if( !hit ) { continue; }
+ const Trk::PrepRawData* hitPRD( hit->prepRawData() );
+ if( hitPRD ) {
+ Identifier theId = hitPRD->identify();
+
+ if( !m_pixelId->is_pixel(theId) && ! m_sctId->is_sct(theId)){
+ ATH_MSG_VERBOSE( "an RIO on track which is not pix or sct");
+ continue;
+ }
+
+ // output the index of the hit
+ std::map<Identifier,int>::const_iterator theIndex = clusterIndexMap.find(theId);
+ if( theIndex != clusterIndexMap.end() ) {
+ (*oflraw) << theIndex->second << ' ';
+ ATH_MSG_VERBOSE( "a silicon cluster identifier was found in the index.");
+ } else {
+ ATH_MSG_VERBOSE( "a silicon cluster identifier was not found in the index.");
+
+ ATH_MSG_VERBOSE("Pixel "<< m_pixelId->barrel_ec(theId) << " "
+ << m_pixelId->layer_disk(theId) <<" "
+ << m_pixelId->phi_module(theId) <<" "
+ << m_pixelId->eta_module(theId) <<" "
+ << m_pixelId->phi_index(theId) << " "
+ << m_pixelId->eta_index(theId));
+ ATH_MSG_VERBOSE("SCT "<< m_sctId->barrel_ec(theId) << " "
+ << m_sctId->layer_disk(theId) <<" "
+ << m_sctId->phi_module(theId) <<" "
+ << m_sctId->eta_module(theId));
+
+ }
+ } else {
+ ATH_MSG_DEBUG( "a silicon cluster on track did not have any prepRawData.");
+ }
+ // if hit RIO on track
+ } // end loop over hits on tracks
+ (*oflraw) << endl;
+
+ // Alberto's code for dumping hits on tracks. this code needs to
+ // be updated to use the modern track classes.
+ if( m_dumpHitsOnTracks ) {
+ // // loop over list of HitOnTrack and dump them
+ // for( hit_const_iterator iHit=track->hit_list_begin(), fHit=track->hit_list_end(); iHit!=fHit; ++iHit ) {
+ // const HitOnTrack* pHit( *iHit );
+ // const Amg::Vector3D global_pos( pHit->position() );
+ // const Trk::TrkDetElementBase* detElem( pHit->detectorElement() );
+ // if( !detElem ) { continue; } // process only hits associated with a detector element
+ // const Identifier rdoId = detElem->identify();
+ // const InDetDD::SiDetectorElement* sielement = 0;
+ // // get the PrepRawData data( cluster for Pixel and SCT hits)
+ // const Trk::PrepRawData* pPrepClu = pHit->rio();
+ // // is this an SCT hit?
+ // sielement = m_SCT_mgr->getDetectorElement( rdoId );
+ // if( sielement && sielement->isSCT() ) { // hit is SCT
+ // // if there is no PrepRawData associated then the hit is
+ // // not a real hit seen by the detector, so skip it.
+ // if( !pPrepClu ) { continue; }
+ // // get local position w/ Lorentz correction
+ // const Amg::Vector2D* locPos = &(pPrepClu->localPosition());
+ // assert( locPos );
+ // // get list of associated hit identifiers
+ // const vector<Identifier>* pHitRdoList = &(pPrepClu->rdoList());
+ // assert( pHitRdoList );
+ // const Amg::Vector2D rawPos = *locPos - sielement->getLorentzCorrection();
+ // const Identifier tmpId = sielement->identifierOfPosition(rawPos);
+ // const Amg::Vector2D cellPos = sielement->localPositionOfCell(tmpId); // w/ Lorentz correction
+ // const Amg::Vector2D deltaPos = *locPos - cellPos; // delta w.r.t. center of cell
+ // (*oflraw) << "e\t"
+ // << setw(14) << setprecision(10)
+ // << global_pos[0] << '\t' << global_pos[1] << '\t' << global_pos[2] << '\t'
+ // << resetiosflags(ios::scientific)
+ // << 0 << '\t'
+ // << m_sctId->barrel_ec(rdoId) << '\t'
+ // << m_sctId->layer_disk(rdoId) << '\t'
+ // << m_sctId->phi_module(rdoId) << '\t'
+ // << m_sctId->eta_module(rdoId) << '\t'
+ // << m_sctId->side(rdoId) << '\t'
+ // << m_sctId->strip(tmpId) << '\t' // use tmpId here
+ // << setiosflags(ios::fixed) << setprecision(4)
+ // << deltaPos[0] << '\t' << deltaPos[1] << '\t' // delta w.r.t. center of cell (Lorentz correction cancels out)
+ // << resetiosflags(ios::scientific)
+ // << pHitRdoList->size() << '\n';
+ // // FlagAA: here the loop is over single strip hits.
+ // // It would be interesting to dump multiple hits
+ // // when they correspond to 2 separate raw hits.
+ // // E.g. when they cross the chip boundary and are thus clustered
+ // // separately at front-end level and into a single cluster offline.
+ // // This is left to be done later on.
+ // // int hitCounter=0;
+ // // for (pAssocHit = pHitRdoList->begin();
+ // // pAssocHit != pHitRdoList->end(); pAssocHit++) {
+ // // if (hitCounter<)
+ // // hitCounter++;
+ // // (*oflraw)<< "\t i=" << *pAssocHit;
+ // // (*oflraw)<< "\thId=" << hitIndexMap[*pAssocHit];
+ // // }
+ // } // end of SCT
+ // sielement = m_PIX_mgr->getDetectorElement( rdoId );
+ // if( sielement && sielement->isPixel() ) { // hit is pixel
+ // // if there is no PrepRawData associated then the hit is
+ // // not a real hit seen by the detector, so skip it.
+ // if( !pPrepClu ) { continue; }
+ // // get local position w/ Lorentz correction
+ // const Amg::Vector2D* locPos = &(pPrepClu->localPosition());
+ // assert( locPos );
+ // // get list of associated hit identifiers
+ // const vector<Identifier>* pHitRdoList = &(pPrepClu->rdoList());
+ // assert( pHitRdoList );
+ // const Amg::Vector2D rawPos = *locPos - sielement->getLorentzCorrection();
+ // const Identifier tmpId = sielement->identifierOfPosition(rawPos);
+ // const Amg::Vector2D cellPos = sielement->localPositionOfCell(tmpId); // w/ Lorentz correction
+ // const Amg::Vector2D deltaPos = *locPos - cellPos; // delta w.r.t. center of cell
+ // (*oflraw) << "e\t"
+ // << setw(14) << setprecision(10)
+ // << global_pos[0] << '\t' << global_pos[1] << '\t' << global_pos[2] << '\t'
+ // << resetiosflags(ios::scientific)
+ // << 1 << '\t'
+ // << m_pixelId->barrel_ec(rdoId) << '\t'
+ // << m_pixelId->layer_disk(rdoId) << '\t'
+ // << m_pixelId->phi_module(rdoId) << '\t'
+ // << m_pixelId->eta_module(rdoId) << '\t'
+ // << m_pixelId->phi_index(tmpId) << '\t' // use tmpId here
+ // << m_pixelId->eta_index(tmpId) << '\t'
+ // << setiosflags(ios::fixed) << setprecision(4)
+ // << deltaPos[0] << '\t' << deltaPos[1] << '\t' // delta w.r.t. center of cell (Lorentz correction cancels out)
+ // << resetiosflags(ios::scientific)
+ // << pHitRdoList->size();
+ // for( vector<Identifier>::const_iterator pAssocHit=pHitRdoList->begin(), fAssocHit=pHitRdoList->end(); pAssocHit!=fAssocHit; ++pAssocHit ) {
+ // (*oflraw) << '\t' << hitIndexMap.find(*pAssocHit)->second;
+ // }
+ // (*oflraw) << '\n';
+ // } // end if hit is pixel, dump
+ // } // end for each hit on this track
+ } // end dump hits on this track
+ } // end for each reconstructed track
+ if( !(trks->empty()) ) { ATH_MSG_DEBUG( "done dumping reconstructed tracks."); }
+}
+
+const DumpSp::ParentBitmask
+DumpSp::construct_truth_bitmap( const HepMC::GenParticle* particle ) const
+{
+ ParentBitmask result;
+ result.reset();
+ typedef pair<const HepMC::GenParticle*,unsigned int> Parent;
+ vector<Parent> parents;
+ parents.push_back( Parent(particle,0) );
+ while( !parents.empty() ) {
+ const HepMC::GenParticle* p = parents.back().first;
+ const unsigned int level = parents.back().second;
+ if( std::abs(p->pdg_id())==15 ) { result.set( TAU_PARENT_BIT , 1 ); }
+ if( std::abs(p->pdg_id())==5 ) { result.set( B_PARENT_BIT , 1 ); }
+ if( std::abs(p->pdg_id())==211 ) { result.set( PION_PARENT_BIT , 1 ); }
+ if( std::abs(p->pdg_id())==211 && level<=1 ) { result.set( PION_IMMEDIATE_PARENT_BIT , 1 ); }
+ if( result.count()==NBITS ) { break; }
+ parents.pop_back();
+ if( !(p->production_vertex()) ) { continue; }
+ for( HepMC::GenVertex::particle_iterator i=p->production_vertex()->particles_begin(HepMC::parents), f=p->production_vertex()->particles_end(HepMC::parents); i!=f; ++i ) {
+ parents.push_back( Parent(*i,level+1) );
+ }
+ }
+ return result;
+}
+
+
+void
+DumpSp::dump_beamspot() const
+{
+ if ( m_outputBeamSpotToWrapper ) { // output beam spot to wrapper
+ (*oflraw) << "V\t" << m_beamCondSvc->beamPos().x() << "\t"
+ << m_beamCondSvc->beamPos().y() << "\t"
+ << m_beamCondSvc->beamPos().z() << "\t"
+ << m_beamCondSvc->beamSigma(0) << "\t"
+ << m_beamCondSvc->beamSigma(1) << "\t"
+ << m_beamCondSvc->beamSigma(2) << "\t"
+ << m_beamCondSvc->beamSigmaXY() << "\t"
+ << m_beamCondSvc->beamTilt(0) << "\t"
+ << m_beamCondSvc->beamTilt(1) << "\n";
+ }
+}
+
+void
+DumpSp::dump_MBTS( ) const
+{
+ // variables needed
+ // number of vertices, MBTS trigger, Jinlong triggers (L1_MBTS_1 or L1_MBTS_1_1)
+ // BC ID, collision timing, etc.
+ string chains_regex(".*(mb|MB).*");
+ vector<string> chains = m_trigDecTool->getChainGroup(chains_regex)->getListOfTriggers();
+ int nChains = chains.size();
+ cout << nChains << " mb chains found." << endl;
+ BOOST_FOREACH( const string& chain , chains ) {
+ cout << " chain: " << chain << " ";
+ const bool ok_physics = m_trigDecTool->isPassed( chain , TrigDefs::Physics );
+ const bool ok_passthru = m_trigDecTool->isPassed( chain , TrigDefs::passedThrough );
+ cout << " " << ok_physics << " " << ok_passthru << endl;
+ }
+
+ const bool ok_physics_A = m_trigDecTool->isPassed( "L1_MBTS_1_1_Col" , TrigDefs::Physics );
+ const bool ok_passthru_A = m_trigDecTool->isPassed( "L1_MBTS_1_1_Col" , TrigDefs::passedThrough );
+ const bool ok_physics_B = m_trigDecTool->isPassed( "L1_MBTS_1_1_EMPTY" , TrigDefs::Physics );
+ const bool ok_passthru_B = m_trigDecTool->isPassed( "L1_MBTS_1_1_EMPTY" , TrigDefs::passedThrough );
+
+ const bool ok_physics_C = m_trigDecTool->isPassed( "L1_MBTS_1_1_UNPARIED" , TrigDefs::Physics );
+ const bool ok_passthru_C = m_trigDecTool->isPassed( "L1_MBTS_1_1_UNPAIRED" , TrigDefs::passedThrough );
+ const bool ok_physics_D = m_trigDecTool->isPassed( "L1_MBTS_1_1_UNPARIED1" , TrigDefs::Physics );
+ const bool ok_passthru_D = m_trigDecTool->isPassed( "L1_MBTS_1_1_UNPAIRED1" , TrigDefs::passedThrough );
+ const bool ok_physics_E = m_trigDecTool->isPassed( "L1_MBTS_1_1_UNPARIED2" , TrigDefs::Physics );
+ const bool ok_passthru_E = m_trigDecTool->isPassed( "L1_MBTS_1_1_UNPAIRED2" , TrigDefs::passedThrough );
+
+ const bool ok_physics_F = m_trigDecTool->isPassed( "L1_MBTS_1_1" , TrigDefs::Physics );
+ const bool ok_passthru_F = m_trigDecTool->isPassed( "L1_MBTS_1_1" , TrigDefs::passedThrough );
+ const bool ok_physics_G = m_trigDecTool->isPassed( "L1_MBTS_1" , TrigDefs::Physics );
+ const bool ok_passthru_G = m_trigDecTool->isPassed( "L1_MBTS_1" , TrigDefs::passedThrough );
+
+ // eXtra info
+ // mbts trigger bits
+ (*oflraw) << "X 1 "
+ << ok_physics_A << " " << ok_passthru_A << " "
+ << ok_physics_B << " " << ok_passthru_B << " "
+ << ok_physics_C << " " << ok_passthru_C << " "
+ << ok_physics_D << " " << ok_passthru_D << " "
+ << ok_physics_E << " " << ok_passthru_E << " "
+ << ok_physics_F << " " << ok_passthru_F << " "
+ << ok_physics_G << " " << ok_passthru_G << " "
+ << endl;
+
+ return;
+}
+
+
+
+void
+DumpSp::dump_lumibcid( const EventID* evid) const
+{
+
+ const int lb = evid->lumi_block();
+ const int bcid = evid->bunch_crossing_id();
+ (*oflraw) << "X 2 "
+ << lb << " " << bcid
+ << endl;
+
+ return;
+}
+
+void
+DumpSp::dump_vertex( ) const
+{
+
+ // number of z vertices
+ const VxContainer* vxes( 0 );
+ if( m_storeGate->retrieve( vxes, "VxPrimaryCandidate" ).isSuccess() ) {
+ if( !(vxes->empty()) ) {
+ for( VxContainer::const_iterator i=vxes->begin(), f=vxes->end(); i!=f; ++i ) {
+ const unsigned int vx_type = (*i)->vertexType();
+ const unsigned int vx_ntracks = (*i)->vxTrackAtVertex()->size();
+ const double vx_x = (*i)->recVertex().position().x();
+ const double vx_y = (*i)->recVertex().position().y();
+ const double vx_z = (*i)->recVertex().position().z();
+ (*oflraw) << "X 3 "
+ << vx_type << " "
+ << vx_ntracks << " "
+ << vx_x << " "
+ << vx_y << " "
+ << vx_z << " "
+ << endl;
+ }
+ }
+ }
+ return;
+
+}
+
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/src/FTKRegionalWrapper.cxx b/Trigger/TrigFTK/FastTrackSimWrap/src/FTKRegionalWrapper.cxx
new file mode 100644
index 00000000000..1d7aa4c443d
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/src/FTKRegionalWrapper.cxx
@@ -0,0 +1,219 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "FastTrackSimWrap/FTKRegionalWrapper.h"
+#include "TrigFTKSim/FTKDataInput.h"
+#include "TrigFTKSim/atlClustering.h"
+
+FTKRegionalWrapper::FTKRegionalWrapper (const std::string& name, ISvcLocator* pSvcLocator) :
+ AthAlgorithm(name, pSvcLocator),
+ m_hitInputTool("FTK_SGHitInput/FTK_SGHitInput"),
+ m_IBLMode(false),
+ m_pmap_path(""),
+ m_pmap(0x0),
+ m_rmap_path(""),
+ m_rmap(0x0),
+ m_ntowers(0),
+ m_SaveRawHits(true),
+ m_SaveHits(false),
+ m_Clustering(false),
+ m_SaveClusterContent(false),
+ m_DiagClustering(true),
+ m_SctClustering(false),
+ m_PixelClusteringMode(1),
+ m_DuplicateGanged(true),
+ m_GangedPatternRecognition(false),
+ m_outpath("ftksim_smartwrapper.root"),
+ m_outfile(0x0),
+ m_hittree(0x0)
+{
+ declareProperty("RMapPath",m_rmap_path);
+ declareProperty("PMapPath",m_pmap_path);
+ declareProperty("OutFileName",m_outpath);
+ declareProperty("HitInputTool",m_hitInputTool);
+ declareProperty("IBLMode",m_IBLMode);
+
+ // hit type options
+ declareProperty("SaveRawHits",m_SaveRawHits);
+ declareProperty("SaveHits",m_SaveHits);
+
+ // clustering options
+ declareProperty("Clustering",m_Clustering);
+ declareProperty("SaveClusterContent",m_SaveClusterContent);
+ declareProperty("DiagClustering",m_DiagClustering);
+ declareProperty("SctClustering",m_SctClustering);
+ declareProperty("PixelClusteringMode",m_PixelClusteringMode);
+ declareProperty("DuplicateGanged",m_DuplicateGanged);
+ declareProperty("GangedPatternRecognition",m_GangedPatternRecognition);
+}
+
+FTKRegionalWrapper::~FTKRegionalWrapper ()
+{
+ if (m_rmap) delete m_rmap;
+}
+
+StatusCode FTKRegionalWrapper::initialize()
+{
+ MsgStream log(msgSvc(), name());
+ log << MSG::INFO << "FTKRegionalWrapper::initialize()" << endreq;
+
+ // FTK library global setup variables
+ FTKSetup::getFTKSetup().setIBLMode(m_IBLMode);
+
+ log << MSG::INFO << "Read the logical layer definitions" << endreq;
+ // Look for the main plane-map
+ if (m_pmap_path.empty()) {
+ log << MSG::FATAL << "Main plane map definition missing" << endreq;
+ return StatusCode::FAILURE;
+ }
+ else {
+ m_pmap = new FTKPlaneMap(m_pmap_path.c_str());
+ if (!(*m_pmap)) {
+ log << MSG::FATAL << "Error using plane map: " << m_pmap_path << endreq;
+ return StatusCode::FAILURE;
+ }
+ }
+
+ // initialize the tower/region map
+ log << MSG::INFO << "Creating region map" << endreq;
+ m_rmap = new FTKRegionMap(m_pmap, m_rmap_path.c_str());
+ if (!(*m_rmap)) {
+ log << MSG::FATAL << "Error creating region map from: " << m_rmap_path.c_str() << endreq;
+ return StatusCode::FAILURE;
+ }
+
+ StatusCode schit = m_hitInputTool.retrieve();
+ if (schit.isFailure()) {
+ log << MSG::FATAL << "Could not retrieve FTK_SGHitInput tool" << endreq;
+ return StatusCode::FAILURE;
+ }
+ else {
+ log << MSG::INFO << "Setting FTK_SGHitInput tool" << endreq;
+ // set the pmap address to FTKDataInput to use in processEvent
+ m_hitInputTool->reference()->setPlaneMaps(m_pmap,0x0);
+ }
+
+ if (!m_SaveRawHits && !m_SaveHits) {
+ log << MSG::FATAL << "At least one hit format has to be saved: FTKRawHit or FTKHit" << endl;
+ return StatusCode::FAILURE;
+ }
+
+ /*
+ * prepare the output structure to store the hits and the other information
+ */
+
+ // create the output files
+ log << MSG::INFO << "Creating output file: " << m_outpath << endreq;
+ m_outfile = TFile::Open(m_outpath.c_str(),"recreate");
+
+ // create a TTree to store event information
+ m_evtinfo = new TTree("evtinfo","Events info");
+ m_evtinfo->Branch("RunNumber",&m_run_number,"RunNumber/I");
+ m_evtinfo->Branch("EventNumber",&m_event_number,"EventNumber/I");
+
+ // create and populate the TTree
+ m_hittree = new TTree("ftkhits","Raw hits for the FTK simulation");
+
+ // prepare a branch for each tower
+ m_ntowers = m_rmap->getNumRegions();
+
+ if (m_SaveRawHits) { // Save FTKRawHit data
+ m_original_hits = new vector<FTKRawHit>[m_ntowers];
+ for (int ireg=0;ireg!=m_ntowers;++ireg) { // towers loop
+ m_hittree->Branch(Form("RawHits%d.",ireg),&m_original_hits[ireg]);
+ } // end towers loop
+ }
+
+ if (m_SaveHits) {
+ m_logical_hits = new vector<FTKHit>[m_ntowers];
+ for (int ireg=0;ireg!=m_ntowers;++ireg) { // towers loop
+ m_hittree->Branch(Form("Hits%d.",ireg),&m_logical_hits[ireg]);
+ } // end towers loop
+ }
+
+ // create a TTree to store the truth tracks
+ m_trackstree = new TTree("truthtracks","Truth tracks");
+ // add the branch related to the truth tracks
+ m_trackstree->Branch("TruthTracks",&m_truth_tracks);
+
+ /* initialize the clustering global variables, decalred in TrigFTKSim/atlClusteringLNF.h */
+ SAVE_CLUSTER_CONTENT = m_SaveClusterContent;
+ DIAG_CLUSTERING = m_DiagClustering;
+ SCT_CLUSTERING = m_SctClustering;
+ PIXEL_CLUSTERING_MODE = m_PixelClusteringMode;
+ DUPLICATE_GANGED = m_DuplicateGanged;
+ GANGED_PATTERN_RECOGNITION = m_GangedPatternRecognition;
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode FTKRegionalWrapper::execute()
+{
+ // retrieve the pointer to the datainput object
+ FTKDataInput *datainput = m_hitInputTool->reference();
+
+ // reset the branches
+ for (int ireg=0;ireg!=m_ntowers;++ireg) {
+ if (m_SaveRawHits) m_original_hits[ireg].clear();
+ if (m_SaveHits) m_logical_hits[ireg].clear();
+ }
+
+ // TODO: test
+ // ask to read the data in the current event
+ datainput->readData();
+
+ // get the event information
+ m_run_number = datainput->runNumber(); // event's run number
+ m_event_number = datainput->eventNumber(); // event number
+ m_evtinfo->Fill();
+
+ // retrieve the original list of hits, the list is copied because the clustering will change it
+ vector<FTKRawHit> fulllist = datainput->getOriginalHits();
+
+ // if the clustering is requested it has to be done before the hits are distributed
+ if (m_Clustering ) {
+ atlClusteringLNF(fulllist);
+ }
+
+ // prepare to iterate on the input files
+ vector<FTKRawHit>::const_iterator ihit = fulllist.begin();
+ vector<FTKRawHit>::const_iterator ihitE = fulllist.end();
+
+ for (;ihit!=ihitE;++ihit) { // hit loop
+ const FTKRawHit &currawhit = *ihit;
+
+ //cout << "Hit " << currawhit.getHitType() << ": " << currawhit.getEventIndex() << " " << currawhit.getBarcode() << endl;
+ // calculate the equivalent hit
+ FTKHit hitref = currawhit.getFTKHit(m_pmap);
+
+ // check the region
+ for (int ireg=0;ireg!=m_ntowers;++ireg) {
+ if (m_rmap->isHitInRegion(hitref,ireg)) {
+ // if the equivalent hit is compatible with this tower the hit is saved
+ if (m_SaveRawHits) m_original_hits[ireg].push_back(currawhit);
+ if (m_SaveHits) m_logical_hits[ireg].push_back(hitref);
+ }
+ }
+ } // end hit loop
+
+ // fill the branches
+ m_hittree->Fill();
+
+ // get the list of the truth tracks
+ m_truth_tracks.clear();
+ const vector<FTKTruthTrack> &truthtracks = datainput->getTruthTrack();
+ m_truth_tracks.insert(m_truth_tracks.end(),truthtracks.begin(),truthtracks.end());
+ // Write the tracks
+ m_trackstree->Fill();
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode FTKRegionalWrapper::finalize()
+{
+ // close the output files
+ m_outfile->Write();
+ m_outfile->Close();
+ return StatusCode::SUCCESS;
+}
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_entries.cxx b/Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_entries.cxx
new file mode 100644
index 00000000000..ae488a04dbc
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_entries.cxx
@@ -0,0 +1,14 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+//#include "GaudiKernel/ServiceHandle.h"
+#include "InDetBeamSpotService/IBeamCondSvc.h"
+#include "FastTrackSimWrap/DumpSp.h"
+#include "FastTrackSimWrap/FTKRegionalWrapper.h"
+
+DECLARE_ALGORITHM_FACTORY( DumpSp )
+DECLARE_ALGORITHM_FACTORY( FTKRegionalWrapper )
+DECLARE_FACTORY_ENTRIES( FastTrackSimWrap )
+{
+ DECLARE_ALGORITHM( DumpSp );
+ DECLARE_ALGORITHM( FTKRegionalWrapper );
+}
+
diff --git a/Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_load.cxx b/Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_load.cxx
new file mode 100644
index 00000000000..de50274baed
--- /dev/null
+++ b/Trigger/TrigFTK/FastTrackSimWrap/src/components/FastTrackSimWrap_load.cxx
@@ -0,0 +1,4 @@
+//
+//
+#include "GaudiKernel/LoadFactoryEntries.h"
+LOAD_FACTORY_ENTRIES( FastTrackSimWrap )
--
GitLab