diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
index bdae4e06d5b5d51123a704ddd77d88f4ce44ab3c..49a186f47a7204009e4e08efd2402d4805ccac6e 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
@@ -60,6 +60,7 @@ atlas_add_root_dictionary( ISF_FastCaloSimParametrizationLib
ISF_FastCaloSimParametrization/TFCSNNLateralShapeParametrization.h
ISF_FastCaloSimParametrization/TFCSSimpleLateralShapeParametrization.h
ISF_FastCaloSimParametrization/TreeReader.h
+ ISF_FastCaloSimParametrization/FCS_Cell.h
Root/LinkDef.h
EXTERNAL_PACKAGES ROOT HepPDT XercesC CLHEP HepMC Geant4 )
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/FCS_Cell.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/FCS_Cell.h
similarity index 80%
rename from Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/FCS_Cell.h
rename to Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/FCS_Cell.h
index 560d45d6c361b5d02301737a7ed04d0b39f5f88b..0999fc7ccaf6feedae7f49d43f6541416a47349b 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/FCS_Cell.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/FCS_Cell.h
@@ -85,42 +85,6 @@ struct FCS_matchedcellvector //this is the matched structure for the whole event
inline float scalingfactor(){float cellsum=0.; float hitsum=0.; for (unsigned int i=0; i<m_vector.size(); i++){cellsum+=m_vector[i].cell.energy;for (unsigned int j=0; j<m_vector[i].hit.size(); j++){hitsum+=m_vector[i].hit[j].hit_energy;};}; return cellsum/hitsum;}; //doesn't check for 0!
};
-struct FCS_truth : public TLorentzVector
-{
- std::vector<double> TTC_entrance_eta;
- std::vector<double> TTC_entrance_phi;
- std::vector<double> TTC_entrance_r;
- std::vector<double> TTC_entrance_z;
- std::vector<double> TTC_back_eta;
- std::vector<double> TTC_back_phi;
- std::vector<double> TTC_back_r;
- std::vector<double> TTC_back_z;
- double TTC_IDCaloBoundary_eta;
- double TTC_IDCaloBoundary_phi;
- double TTC_IDCaloBoundary_r;
- double TTC_IDCaloBoundary_z;
- double TTC_Angle3D;
- double TTC_AngleEta;
- int barcode;
- int vtxbarcode;
- int pdgid;
- ClassDef (FCS_truth, 1)
-};
-
-
-#ifdef __CINT__
-#pragma link C++ struct FCS_cell+;
-#pragma link C++ struct FCS_hit+;
-#pragma link C++ struct FCS_g4hit+;
-#pragma link C++ struct std::vector<FCS_hit>+;
-#pragma link C++ struct std::vector<FCS_g4hit>+;
-#pragma link C++ struct FCS_matchedcell+;
-#pragma link C++ struct std::vector<FCS_matchedcell>+;
-#pragma link C++ struct FCS_matchedcellvector+;
-#pragma link C++ class std::vector<Float_t>+;
-#pragma link C++ struct FCS_truth+;
-#pragma link C++ struct std::vector<FCS_truth>+;
-#endif
#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ISF_HitAnalysis.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ISF_HitAnalysis.h
index f2972e3f380b5e014e661053878ba2ba742ca231..f07595ec60cd6469e10d0dfac5cb6ee742ea612a 100755
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ISF_HitAnalysis.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ISF_HitAnalysis.h
@@ -28,6 +28,8 @@
#include "ISF_FastCaloSimParametrization/IFastCaloSimCaloExtrapolation.h"
#include "ISF_FastCaloSimParametrization/IFastCaloSimGeometryHelper.h"
#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
+#include "ISF_FastCaloSimParametrization/FCS_Cell.h"
+
namespace Trk
{
@@ -40,6 +42,7 @@ class ICaloSurfaceHelper;
#include <string>
#include <Rtypes.h>
+#include <TLorentzVector.h>
//#include "TH1.h"
/* *************************************************************
@@ -79,13 +82,10 @@ class ISF_HitAnalysis : public AthAlgorithm {
virtual StatusCode geoInit(IOVSVC_CALLBACK_ARGS);
virtual StatusCode updateMetaData(IOVSVC_CALLBACK_ARGS);
- //bool get_calo_etaphi(std::vector<Trk::HitInfo>* hitVector,CaloCell_ID_FCS::CaloSample sample);
- bool get_calo_etaphi(std::vector<Trk::HitInfo>* hitVector,int sample,int subpos=SUBPOS_MID);
- bool get_calo_surface(std::vector<Trk::HitInfo>* hitVector);
- bool rz_cylinder_get_calo_etaphi(std::vector<Trk::HitInfo>* hitVector, double cylR, double cylZ, Amg::Vector3D& pos, Amg::Vector3D& mom);
-
IFastCaloSimGeometryHelper* GetCaloGeometry() const {return &(*m_CaloGeometryHelper);};
+ const static int MAX_LAYER = 25;
+
private:
/** a handle on Store Gate for access to the Event Store */
@@ -135,12 +135,22 @@ class ISF_HitAnalysis : public AthAlgorithm {
std::vector<Long64_t>* m_g4hit_cellidentifier;
std::vector<float>* m_g4hit_samplingfraction;
std::vector<int>* m_g4hit_sampling;
- //Ok, this won't work, ROOT won't let me save a custom object which it doesn't know about
- //std::vector<zh_matchedcell>* m_matched_cells;
+
+ //CaloHitAna variables
+ FCS_matchedcellvector* m_oneeventcells; //these are all matched cells in a single event
+ FCS_matchedcellvector* m_layercells[MAX_LAYER]; //these are all matched cells in a given layer in a given event
+
+ Float_t m_total_cell_e = 0;
+ Float_t m_total_hit_e = 0;
+ Float_t m_total_g4hit_e = 0;
+
+ std::vector<Float_t>* m_final_cell_energy;
+ std::vector<Float_t>* m_final_hit_energy;
+ std::vector<Float_t>* m_final_g4hit_energy;
+
TTree * m_tree;
std::string m_ntupleFileName;
- std::string m_ntupleDirName;
std::string m_ntupleTreeName;
std::string m_metadataTreeName;
std::string m_geoFileName;
@@ -246,6 +256,13 @@ class ISF_HitAnalysis : public AthAlgorithm {
double m_CaloBoundaryR;
double m_CaloBoundaryZ;
double m_calomargin;
+ bool m_saveAllBranches;
+ bool m_doAllCells;
+ bool m_doLayers;
+ bool m_doLayerSums;
+ bool m_doG4Hits;
+ Int_t m_TimingCut;
+
std::string m_MC_DIGI_PARAM;
std::string m_MC_SIM_PARAM;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h
index 5df2b861b1c9c3ba0b4aa25dc8416540c799ae48..2dde8bd432bf47b1d1ecf1a91c4e35974204f65e 100755
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h
@@ -9,6 +9,18 @@
#pragma link C++ class MeanAndRMS;
#pragma link C++ class TreeReader;
#pragma link C++ class EnergyParametrizationValidation;
+#pragma link C++ struct FCS_cell+;
+#pragma link C++ struct FCS_hit+;
+#pragma link C++ struct FCS_g4hit+;
+#pragma link C++ struct std::vector<FCS_hit>+;
+#pragma link C++ struct std::vector<FCS_g4hit>+;
+#pragma link C++ struct FCS_matchedcell+;
+#pragma link C++ struct std::vector<FCS_matchedcell>+;
+#pragma link C++ struct FCS_matchedcellvector+;
+#pragma link C++ class std::vector<Float_t>+;
+#pragma link C++ struct FCS_truth+;
+#pragma link C++ struct std::vector<FCS_truth>+;
+#pragma link C++ class std::vector<std::vector<float> >+;
#ifndef CaloGeometryFromFile_h
#pragma link C++ class CaloGeometryLookup;
#pragma link C++ class CaloGeometry;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple.py
index e3fcc92a05aaa96197aa4446c69531debce89341..d36533ce642ca38994748bddc626b4ed68a7c537 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple.py
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_ntuple.py
@@ -16,7 +16,8 @@ from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
#athenaCommonFlags.FilesInput = glob( "ESD_*root" )
#athenaCommonFlags.FilesInput = ["/afs/cern.ch/user/c/cmills/public/pions20GeV_fulldet.ESD.pool.root"]
#athenaCommonFlags.FilesInput = ["/afs/cern.ch/user/c/cmills/public/pions20GeV_z0150_fulldet.ESD.pool.root"]
-athenaCommonFlags.FilesInput = ["root://eosatlas//eos/atlas/user/z/zhubacek/FastCaloSim/ForMichael/ESD_evgen_calo__211_E50000_50000_eta20_25_Evts0-5500_vz_0_origin_calo.pool.root"]
+#athenaCommonFlags.FilesInput = ["root://eosatlas//eos/atlas/user/z/zhubacek/FastCaloSim/ForMichael/ESD_evgen_calo__211_E50000_50000_eta20_25_Evts0-5500_vz_0_origin_calo.pool.root"]
+athenaCommonFlags.FilesInput = ["/afs/cern.ch/work/a/ahasib/public/photon.50GeV.ESD.pool.root"]
ServiceMgr.EventSelector.InputCollections = athenaCommonFlags.FilesInput() # This is stupid and redundant, but necessary
@@ -56,7 +57,13 @@ ISF_HitAnalysis.CaloBoundaryZ = 3549.5 #before: 3475.0
ISF_HitAnalysis.CaloMargin=100 #=10cm
ISF_HitAnalysis.NTruthParticles = 1 # Copy only one truth particle to the ntuples for now
#ISF_HitAnalysis.OutputLevel = WARNING
-ISF_HitAnalysis.OutputLevel = ERROR
+ISF_HitAnalysis.SaveAllBranches = False
+ISF_HitAnalysis.DoAllCells = False
+ISF_HitAnalysis.DoLayers = True
+ISF_HitAnalysis.DoLayerSums = True
+ISF_HitAnalysis.DoG4Hits = False
+ISF_HitAnalysis.TimingCut = 999999
+ISF_HitAnalysis.OutputLevel = DEBUG
#############################
##### NEW TRACKING SETUP ####
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/ISF_HitAnalysis.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/ISF_HitAnalysis.cxx
index 5636cccef78f289a1838b5757e1c8c20fcd206e8..f5de1ae83d100bf39ae83a01ce1509545169c133 100755
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/ISF_HitAnalysis.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/ISF_HitAnalysis.cxx
@@ -45,7 +45,6 @@
#include "TString.h"
#include "TVector3.h"
#include <sstream>
-//#include "TLorentzVector.h"
// For MC Truth information:
#include "GeneratorObjects/McEventCollection.h"
@@ -54,7 +53,6 @@
//####################
#include "GaudiKernel/ListItem.h"
#include "CaloDetDescr/CaloDepthTool.h"
-//#include "RecoToolInterfaces/IExtrapolateToCaloTool.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkSurfaces/CylinderSurface.h"
#include "TrkSurfaces/DiscSurface.h"
@@ -75,10 +73,6 @@
#include <functional>
#include <iostream>
-
-//const std::string MC_DIGI_PARAM = "/Digitization/Parameters";
-//const std::string MC_SIM_PARAM = "/Simulation/Parameters";
-
ISF_HitAnalysis::ISF_HitAnalysis(const std::string& name, ISvcLocator* pSvcLocator)
: AthAlgorithm(name, pSvcLocator)
//, m_storeGate(0)
@@ -119,10 +113,14 @@ ISF_HitAnalysis::ISF_HitAnalysis(const std::string& name, ISvcLocator* pSvcLocat
, m_g4hit_cellidentifier(0)
, m_g4hit_samplingfraction(0)
, m_g4hit_sampling(0)
- // , m_matched_cells(0)
+ , m_total_cell_e(0)
+ , m_total_hit_e(0)
+ , m_total_g4hit_e(0)
+ , m_final_cell_energy(0)
+ , m_final_hit_energy(0)
+ , m_final_g4hit_energy(0)
, m_tree(0)
- , m_ntupleFileName("/ntuples/file1")
- , m_ntupleDirName("ISF_HitAnalysis")
+ , m_ntupleFileName("ISF_HitAnalysis")
, m_ntupleTreeName("CaloHitAna")
, m_metadataTreeName("MetaData")
, m_geoFileName("ISF_Geometry")
@@ -139,25 +137,6 @@ ISF_HitAnalysis::ISF_HitAnalysis(const std::string& name, ISvcLocator* pSvcLocat
, m_ptruth_pt(0)
, m_ptruth_p(0)
, m_pdgid(0)
-
- //######################
-
-/*
- , m_TTC_entrance_eta(0)
- , m_TTC_entrance_phi(0)
- , m_TTC_entrance_r(0)
- , m_TTC_entrance_z(0)
- , m_TTC_back_eta(0)
- , m_TTC_back_phi(0)
- , m_TTC_back_r(0)
- , m_TTC_back_z(0)
- , m_TTC_IDCaloBoundary_eta(0)
- , m_TTC_IDCaloBoundary_phi(0)
- , m_TTC_IDCaloBoundary_r(0)
- , m_TTC_IDCaloBoundary_z(0)
- , m_TTC_Angle3D(0)
- , m_TTC_AngleEta(0)
-*/
, m_newTTC_entrance_eta(0)
, m_newTTC_entrance_phi(0)
, m_newTTC_entrance_r(0)
@@ -191,7 +170,6 @@ ISF_HitAnalysis::ISF_HitAnalysis(const std::string& name, ISvcLocator* pSvcLocat
//Note that m_xxx are pointers to vectors set to 0, not set to empty vector! see note around TBranch
{
declareProperty("NtupleFileName", m_ntupleFileName);
- declareProperty("NtupleDirectoryName", m_ntupleDirName);
declareProperty("NtupleTreeName", m_ntupleTreeName);
declareProperty("GeoFileName", m_geoFileName);
declareProperty("MetadataTreeName", m_metadataTreeName);
@@ -219,6 +197,16 @@ ISF_HitAnalysis::ISF_HitAnalysis(const std::string& name, ISvcLocator* pSvcLocat
declareProperty("MetaDataSim", m_MC_SIM_PARAM );
declareProperty("MetaDataDigi", m_MC_DIGI_PARAM );
+ declareProperty("SaveAllBranches", m_saveAllBranches = false);
+ declareProperty("DoAllCells", m_doAllCells = false);
+ declareProperty("DoLayers", m_doLayers = true);
+ declareProperty("DoLayerSums", m_doLayerSums = true);
+ declareProperty("DoG4Hits", m_doG4Hits = false);
+ declareProperty("TimingCut", m_TimingCut = 999999);
+
+
+
+
m_surfacelist.resize(0);
m_surfacelist.push_back(CaloCell_ID_FCS::PreSamplerB);
m_surfacelist.push_back(CaloCell_ID_FCS::PreSamplerE);
@@ -371,22 +359,6 @@ StatusCode ISF_HitAnalysis::updateMetaData( IOVSVC_CALLBACK_ARGS_P( I, keys ) )
StatusCode ISF_HitAnalysis::initialize()
{
ATH_MSG_INFO( "Initializing ISF_HitAnalysis" );
-
- /*
- // get a handle of StoreGate for access to the Detector Store
- sc = service("DetectorStore", m_detStore);
- if (sc.isFailure()) {
- ATH_MSG_ERROR( "ZH Unable to retrieve pointer to Detector Store");
- return sc;
- }
- // get a handle of StoreGate for access to the Event Store
- sc = service("StoreGateSvc", m_storeGate);
- if (sc.isFailure()) {
- ATH_MSG_ERROR( "Unable to retrieve pointer to StoreGateSvc" );
- return sc;
- }
- */
-
//
// Register the callback(s):
//
@@ -483,10 +455,8 @@ StatusCode ISF_HitAnalysis::initialize()
return StatusCode::FAILURE;
}
//#########################
-
- //
- m_tree = new TTree( TString(m_ntupleTreeName), "CaloHitAna" );
- std::string fullNtupleName = "/"+m_ntupleFileName+"/"+m_ntupleDirName+"/"+m_ntupleTreeName;
+ m_tree = new TTree("FCS_ParametrizationInput", "FCS_ParametrizationInput");
+ std::string fullNtupleName = "/"+m_ntupleFileName+"/"+m_ntupleTreeName;
sc = m_thistSvc->regTree(fullNtupleName, m_tree);
if (sc.isFailure() || !m_tree )
{
@@ -498,22 +468,98 @@ StatusCode ISF_HitAnalysis::initialize()
if (m_tree)
{
ATH_MSG_INFO("Successfull registered TTree: " << fullNtupleName);
- //this actually creates the vector itself! And only if it succeeds! Note that the result is not checked! And the code is probably leaking memory in the end
- m_tree->Branch("HitX", &m_hit_x);
- m_tree->Branch("HitY", &m_hit_y);
- m_tree->Branch("HitZ", &m_hit_z);
- m_tree->Branch("HitE", &m_hit_energy);
- m_tree->Branch("HitT", &m_hit_time);
- m_tree->Branch("HitIdentifier", &m_hit_identifier);
- m_tree->Branch("HitCellIdentifier", &m_hit_cellidentifier);
- m_tree->Branch("HitIsLArBarrel", &m_islarbarrel);
- m_tree->Branch("HitIsLArEndCap", &m_islarendcap);
- m_tree->Branch("HitIsHEC", &m_islarhec);
- m_tree->Branch("HitIsFCAL", &m_islarfcal);
- m_tree->Branch("HitIsTile", &m_istile);
- m_tree->Branch("HitSampling", &m_hit_sampling);
- m_tree->Branch("HitSamplingFraction", &m_hit_samplingfraction);
+ //initialize the variables before creating the branches
+ m_hit_x = new std::vector<float>;
+ m_hit_y = new std::vector<float>;
+ m_hit_z = new std::vector<float>;
+ m_hit_energy = new std::vector<float>;
+ m_hit_time = new std::vector<float>;
+ m_hit_identifier = new std::vector<Long64_t>;
+ m_hit_cellidentifier = new std::vector<Long64_t>;
+ m_islarbarrel = new std::vector<bool>;
+ m_islarendcap = new std::vector<bool>;
+ m_islarhec = new std::vector<bool>;
+ m_islarfcal = new std::vector<bool>;
+ m_istile = new std::vector<bool>;
+ m_hit_sampling = new std::vector<int>;
+ m_hit_samplingfraction = new std::vector<float>;
+
+ m_truth_energy = new std::vector<float>;
+ m_truth_px = new std::vector<float>;
+ m_truth_py = new std::vector<float>;
+ m_truth_pz = new std::vector<float>;
+ m_truth_pdg = new std::vector<int>;
+ m_truth_barcode = new std::vector<int>;
+ m_truth_vtxbarcode = new std::vector<int>;
+
+ m_cell_identifier = new std::vector<Long64_t>;
+ m_cell_energy = new std::vector<float>;
+ m_cell_sampling = new std::vector<int>;
+
+ m_g4hit_energy = new std::vector<float>;
+ m_g4hit_time = new std::vector<float>;
+ m_g4hit_identifier = new std::vector<Long64_t>;
+ m_g4hit_cellidentifier = new std::vector<Long64_t>;
+ m_g4hit_samplingfraction = new std::vector<float>;
+ m_g4hit_sampling = new std::vector<int>;
+
+ m_total_cell_e = 0;
+ m_total_hit_e = 0;
+ m_total_g4hit_e = 0;
+
+ m_final_cell_energy = new std::vector<Float_t>;
+ m_final_hit_energy = new std::vector<Float_t>;
+ m_final_g4hit_energy = new std::vector<Float_t>;
+
+ m_newTTC_entrance_eta = new std::vector<std::vector<float> >;
+ m_newTTC_entrance_phi = new std::vector<std::vector<float> >;
+ m_newTTC_entrance_r = new std::vector<std::vector<float> >;
+ m_newTTC_entrance_z = new std::vector<std::vector<float> >;
+ m_newTTC_back_eta = new std::vector<std::vector<float> >;
+ m_newTTC_back_phi = new std::vector<std::vector<float> >;
+ m_newTTC_back_r = new std::vector<std::vector<float> >;
+ m_newTTC_back_z = new std::vector<std::vector<float> >;
+ m_newTTC_mid_eta = new std::vector<std::vector<float> >;
+ m_newTTC_mid_phi = new std::vector<std::vector<float> >;
+ m_newTTC_mid_r = new std::vector<std::vector<float> >;
+ m_newTTC_mid_z = new std::vector<std::vector<float> >;
+ m_newTTC_IDCaloBoundary_eta = new std::vector<float>;
+ m_newTTC_IDCaloBoundary_phi = new std::vector<float>;
+ m_newTTC_IDCaloBoundary_r = new std::vector<float>;
+ m_newTTC_IDCaloBoundary_z = new std::vector<float>;
+ m_newTTC_Angle3D = new std::vector<float>;
+ m_newTTC_AngleEta = new std::vector<float>;
+
+ // Optional branches
+ if(m_saveAllBranches){
+ m_tree->Branch("HitX", &m_hit_x);
+ m_tree->Branch("HitY", &m_hit_y);
+ m_tree->Branch("HitZ", &m_hit_z);
+ m_tree->Branch("HitE", &m_hit_energy);
+ m_tree->Branch("HitT", &m_hit_time);
+ m_tree->Branch("HitIdentifier", &m_hit_identifier);
+ m_tree->Branch("HitCellIdentifier", &m_hit_cellidentifier);
+ m_tree->Branch("HitIsLArBarrel", &m_islarbarrel);
+ m_tree->Branch("HitIsLArEndCap", &m_islarendcap);
+ m_tree->Branch("HitIsHEC", &m_islarhec);
+ m_tree->Branch("HitIsFCAL", &m_islarfcal);
+ m_tree->Branch("HitIsTile", &m_istile);
+ m_tree->Branch("HitSampling", &m_hit_sampling);
+ m_tree->Branch("HitSamplingFraction", &m_hit_samplingfraction);
+
+ m_tree->Branch("CellIdentifier", &m_cell_identifier);
+ m_tree->Branch("CellE", &m_cell_energy);
+ m_tree->Branch("CellSampling", &m_cell_sampling);
+
+ m_tree->Branch("G4HitE", &m_g4hit_energy);
+ m_tree->Branch("G4HitT", &m_g4hit_time);
+ m_tree->Branch("G4HitIdentifier", &m_g4hit_identifier);
+ m_tree->Branch("G4HitCellIdentifier", &m_g4hit_cellidentifier);
+ m_tree->Branch("G4HitSamplingFraction",&m_g4hit_samplingfraction);
+ m_tree->Branch("G4HitSampling", &m_g4hit_sampling);
+ }
+ //CaloHitAna output variables
m_tree->Branch("TruthE", &m_truth_energy);
m_tree->Branch("TruthPx", &m_truth_px);
m_tree->Branch("TruthPy", &m_truth_py);
@@ -522,36 +568,34 @@ StatusCode ISF_HitAnalysis::initialize()
m_tree->Branch("TruthBarcode", &m_truth_barcode);
m_tree->Branch("TruthVtxBarcode", &m_truth_vtxbarcode);
- m_tree->Branch("CellIdentifier", &m_cell_identifier);
- m_tree->Branch("CellE", &m_cell_energy);
- m_tree->Branch("CellSampling", &m_cell_sampling);
-
- m_tree->Branch("G4HitE", &m_g4hit_energy);
- m_tree->Branch("G4HitT", &m_g4hit_time);
- m_tree->Branch("G4HitIdentifier", &m_g4hit_identifier);
- m_tree->Branch("G4HitCellIdentifier", &m_g4hit_cellidentifier);
- m_tree->Branch("G4HitSamplingFraction",&m_g4hit_samplingfraction);
- m_tree->Branch("G4HitSampling", &m_g4hit_sampling);
- //And this looks like will fail since ROOT has no idea what the object is...
- //m_tree->Branch("MatchedCells", &m_matched_cells);
-
- //#########################
- /*
- m_tree->Branch("TTC_back_eta",&m_TTC_back_eta);
- m_tree->Branch("TTC_back_phi",&m_TTC_back_phi);
- m_tree->Branch("TTC_back_r",&m_TTC_back_r);
- m_tree->Branch("TTC_back_z",&m_TTC_back_z);
- m_tree->Branch("TTC_entrance_eta",&m_TTC_entrance_eta);
- m_tree->Branch("TTC_entrance_phi",&m_TTC_entrance_phi);
- m_tree->Branch("TTC_entrance_r",&m_TTC_entrance_r);
- m_tree->Branch("TTC_entrance_z",&m_TTC_entrance_z);
- m_tree->Branch("TTC_IDCaloBoundary_eta",&m_TTC_IDCaloBoundary_eta);
- m_tree->Branch("TTC_IDCaloBoundary_phi",&m_TTC_IDCaloBoundary_phi);
- m_tree->Branch("TTC_IDCaloBoundary_r",&m_TTC_IDCaloBoundary_r);
- m_tree->Branch("TTC_IDCaloBoundary_z",&m_TTC_IDCaloBoundary_z);
- m_tree->Branch("TTC_Angle3D",&m_TTC_Angle3D);
- m_tree->Branch("TTC_AngleEta",&m_TTC_AngleEta);
- */
+ m_oneeventcells = new FCS_matchedcellvector;
+ if(m_doAllCells){
+ m_tree->Branch("AllCells", &m_oneeventcells);
+ }
+
+ //write cells per layer
+ if(m_doLayers){
+ for (Int_t i = 0; i < MAX_LAYER; i++)
+ {
+ TString branchname = "Sampling_";
+ branchname += i;
+ m_layercells[i] = new FCS_matchedcellvector;
+ m_tree->Branch(branchname, &m_layercells[i]);
+ }
+ }
+
+ if(m_doLayerSums){
+ //write also energies per layer:
+ m_tree->Branch("cell_energy", &m_final_cell_energy);
+ m_tree->Branch("hit_energy", &m_final_hit_energy);
+ m_tree->Branch("g4hit_energy", &m_final_g4hit_energy);
+
+ //This is a duplicate of cell_energy[25]
+ m_tree->Branch("total_cell_energy", &m_total_cell_e);
+ m_tree->Branch("total_hit_energy", &m_total_hit_e);
+ m_tree->Branch("total_g4hit_energy", &m_total_g4hit_e);
+ }
+
m_tree->Branch("newTTC_back_eta",&m_newTTC_back_eta);
m_tree->Branch("newTTC_back_phi",&m_newTTC_back_phi);
m_tree->Branch("newTTC_back_r",&m_newTTC_back_r);
@@ -571,31 +615,6 @@ StatusCode ISF_HitAnalysis::initialize()
m_tree->Branch("newTTC_Angle3D",&m_newTTC_Angle3D);
m_tree->Branch("newTTC_AngleEta",&m_newTTC_AngleEta);
- //#########################
- //m_tree->Print();
- if (!m_hit_x || !m_hit_y || !m_hit_z || !m_hit_energy || !m_hit_time || !m_hit_identifier || !m_hit_cellidentifier || !m_islarbarrel || !m_islarendcap || !m_islarfcal || !m_islarhec || !m_istile || !m_hit_sampling || !m_hit_samplingfraction || !m_truth_energy || !m_truth_px || !m_truth_py || !m_truth_pz || !m_truth_pdg || !m_truth_barcode || !m_truth_vtxbarcode)
- {
- ATH_MSG_ERROR("Unable to create TTree branch correctly");
- return StatusCode::FAILURE;
- }
- if (!m_cell_identifier || !m_cell_energy || !m_cell_sampling || !m_g4hit_energy || !m_g4hit_time || !m_g4hit_identifier || !m_g4hit_cellidentifier || !m_g4hit_samplingfraction || !m_g4hit_sampling )
- {
- ATH_MSG_ERROR("Unable to create TTree branch correctly (cell or g4hit)");
- return StatusCode::FAILURE;
- }
- /*
- if (!m_TTC_back_eta || !m_TTC_back_phi || !m_TTC_back_r || !m_TTC_back_z || !m_TTC_entrance_eta || !m_TTC_entrance_phi || !m_TTC_entrance_r || !m_TTC_entrance_z || !m_TTC_IDCaloBoundary_eta || !m_TTC_IDCaloBoundary_phi || !m_TTC_IDCaloBoundary_r || !m_TTC_IDCaloBoundary_z || !m_TTC_Angle3D || !m_TTC_AngleEta)
- {
- ATH_MSG_ERROR("Unable to create TTree branch correctly (TTC variables)");
- return StatusCode::FAILURE;
- }
- */
- if (!m_newTTC_back_eta || !m_newTTC_back_phi || !m_newTTC_back_r || !m_newTTC_back_z || !m_newTTC_entrance_eta || !m_newTTC_entrance_phi || !m_newTTC_entrance_r || !m_newTTC_entrance_z || !m_newTTC_mid_eta || !m_newTTC_mid_phi || !m_newTTC_mid_r || !m_newTTC_mid_z || !m_newTTC_IDCaloBoundary_eta || !m_newTTC_IDCaloBoundary_phi || !m_newTTC_IDCaloBoundary_r || !m_newTTC_IDCaloBoundary_z || !m_newTTC_Angle3D || !m_newTTC_AngleEta)
- {
- ATH_MSG_ERROR("Unable to create TTree branch correctly (newTTC variables)");
- return StatusCode::FAILURE;
- }
-
}
return StatusCode::SUCCESS;
@@ -750,24 +769,23 @@ StatusCode ISF_HitAnalysis::execute()
m_g4hit_samplingfraction->clear();
//which fails for this one!!
//m_matched_cells->clear();
+ std::map<Long64_t, FCS_cell> cells; //read all objects and collect them by identifier (Long64_t)
+ std::map<Long64_t, std::vector<FCS_g4hit> > g4hits;
+ std::map<Long64_t, std::vector<FCS_hit> > hits;
- //##########################
- /*
- m_TTC_back_eta->clear();
- m_TTC_back_phi->clear();
- m_TTC_back_r->clear();
- m_TTC_back_z->clear();
- m_TTC_entrance_eta->clear();
- m_TTC_entrance_phi->clear();
- m_TTC_entrance_r->clear();
- m_TTC_entrance_z->clear();
- m_TTC_IDCaloBoundary_eta->clear();
- m_TTC_IDCaloBoundary_phi->clear();
- m_TTC_IDCaloBoundary_r->clear();
- m_TTC_IDCaloBoundary_z->clear();
- m_TTC_Angle3D->clear();
- m_TTC_AngleEta->clear();
- */
+ cells.clear();
+ g4hits.clear();
+ hits.clear();
+
+ FCS_cell one_cell; //note that this is not extra safe if I don't have a clear method!
+ FCS_g4hit one_g4hit;
+ FCS_hit one_hit;
+ FCS_matchedcell one_matchedcell;
+
+ m_oneeventcells->m_vector.clear();
+ m_final_g4hit_energy->clear();
+ m_final_hit_energy->clear();
+ m_final_cell_energy->clear();
m_newTTC_back_eta->clear();
m_newTTC_back_phi->clear();
@@ -1008,79 +1026,6 @@ StatusCode ISF_HitAnalysis::execute()
m_truth_pdg->push_back((*it)->pdg_id());
m_truth_barcode->push_back((*it)->barcode());
- //extrapolate only if the vertex is within the IDcaloboundary - margin
- //bool inside_ID=false;
-
- /*
-
- HepMC::GenVertex* pvtx = (*it)->production_vertex();
- if(pvtx)
- {
- truth.set_vertex(pvtx->position().x(),pvtx->position().y(),pvtx->position().z(),pvtx->position().t());
-
- m_truth_vtxbarcode->push_back(pvtx->barcode());
- double radius2=(pvtx->position().x()*pvtx->position().x())+(pvtx->position().y()*pvtx->position().y());
- ATH_MSG_VERBOSE("particle "<<particleIndex<<" Mom: "<<(*it)->momentum().e()<<"|"<<(*it)->momentum().px()<<"|"<<(*it)->momentum().py()<<"|"<<(*it)->momentum().pz()<<" vertex_x "<<pvtx->position().x()<<" y "<<pvtx->position().y()<<" z "<<pvtx->position().z()<<" r2 "<<radius2<<" VertexBarcode: "<<pvtx->barcode()<<" ParticleBarcode:"<<(*it)->barcode()<<" status: "<<(*it)->status());
- //std::cout<<"m_CaloBoundaryZ+m_calomargin "<<m_CaloBoundaryZ+m_calomargin<<" m_CaloBoundaryR+m_calomargin "<<m_CaloBoundaryR+m_calomargin<<std::endl;
-
- if(fabs(pvtx->position().z())<(m_CaloBoundaryZ+m_calomargin) && radius2<((m_CaloBoundaryR+m_calomargin)*(m_CaloBoundaryR+m_calomargin)))
- {
- //inside_ID=true;
- //if(inside_ID)
- // {
- ATH_MSG_DEBUG("*** Do extrapolation ***");
- extrapolate(*it,hitVector);
- ATH_MSG_DEBUG("*** Do extrapolation to ID ***");
- extrapolate_to_ID(*it,hitVector);
- ATH_MSG_DEBUG("Truth extrapolation done for "<<particleIndex<<" "<<(*it)->barcode());
- } //in this case no extrapolation possible
- else
- {
- //fill extrapolated angles with some weird values
- m_TTC_IDCaloBoundary_eta->push_back(-999.);
- m_TTC_IDCaloBoundary_phi->push_back(-999.);
- m_TTC_IDCaloBoundary_r->push_back(-999.);
- m_TTC_IDCaloBoundary_z->push_back(-999.);
- m_TTC_Angle3D->push_back(-999.);
- m_TTC_AngleEta->push_back(-999.);
- }
- } //if pvtx
- else //no pvtx
- {
- ATH_MSG_WARNING( "No vertex found for truth particle, no extrapolation");
- m_truth_vtxbarcode->push_back(-999999);
-
- //fill extrapolated angles with some weird values
- std::vector<float> eta_safe;
- std::vector<float> phi_safe;
- std::vector<float> r_safe;
- std::vector<float> z_safe;
- for(int sample=CaloCell_ID_FCS::FirstSample;sample<CaloCell_ID_FCS::MaxSample;++sample)
- {
- eta_safe.push_back(-999.0);
- phi_safe.push_back(-999.0);
- r_safe.push_back(-999.0);
- z_safe.push_back(-999.0);
- }
-
- m_TTC_back_eta->push_back(eta_safe);
- m_TTC_back_phi->push_back(phi_safe);
- m_TTC_back_r->push_back(r_safe);
- m_TTC_back_z->push_back(z_safe);
- m_TTC_entrance_eta->push_back(eta_safe);
- m_TTC_entrance_phi->push_back(phi_safe);
- m_TTC_entrance_r->push_back(r_safe);
- m_TTC_entrance_z->push_back(z_safe);
- m_TTC_IDCaloBoundary_eta->push_back(-999.);
- m_TTC_IDCaloBoundary_phi->push_back(-999.);
- m_TTC_IDCaloBoundary_r->push_back(-999.);
- m_TTC_IDCaloBoundary_z->push_back(-999.);
- m_TTC_Angle3D->push_back(-999.);
- m_TTC_AngleEta->push_back(-999.);
-
- } //no pvtx
- */
-
for(std::vector<Trk::HitInfo>::iterator it = hitVector->begin();it < hitVector->end();++it)
{
if((*it).trackParms)
@@ -1199,6 +1144,344 @@ StatusCode ISF_HitAnalysis::execute()
ATH_MSG_INFO( "Read "<<hitnumber<<" G4Hits from TileHitVec");
}
+
+ // CaloHitAna
+ ATH_MSG_DEBUG("CaloHitAna begin!");
+
+ //cells
+ for (unsigned int cell_i = 0; cell_i < m_cell_identifier->size(); cell_i++)
+ {
+ if (cells.find((*m_cell_identifier)[cell_i]) == cells.end()) //doesn't exist
+ {
+ one_cell.cell_identifier = (*m_cell_identifier)[cell_i];
+ one_cell.sampling = (*m_cell_sampling)[cell_i];
+ one_cell.energy = (*m_cell_energy)[cell_i];
+ one_cell.center_x = 0.0; //for now
+ one_cell.center_y = 0.0;
+ one_cell.center_z = 0.0;
+ cells.insert(std::pair<Long64_t, FCS_cell>(one_cell.cell_identifier, one_cell));
+ }
+ else
+ {
+ //there shouldn't be a cell with the same identifier in this event
+ ATH_MSG_DEBUG("ISF_HitAnalysis: Same cell???? ERROR");
+ }
+ }
+
+ // g4 hits
+ if(m_doG4Hits){
+ for (unsigned int g4hit_i = 0; g4hit_i < m_g4hit_identifier->size(); g4hit_i++)
+ {
+ if ((*m_g4hit_sampling)[g4hit_i] >= 0 && (*m_g4hit_sampling)[g4hit_i] <= 25 && (*m_g4hit_time)[g4hit_i] > m_TimingCut)
+ {
+ ATH_MSG_DEBUG("Ignoring G4hit, time too large: " << g4hit_i << " time: " << (*m_g4hit_time)[g4hit_i]);
+ continue;
+ }
+
+ if (g4hits.find((*m_g4hit_cellidentifier)[g4hit_i]) == g4hits.end())
+ {
+ //this G4 hit doesn't exist yet
+ one_g4hit.identifier = (*m_g4hit_identifier)[g4hit_i];
+ one_g4hit.cell_identifier = (*m_g4hit_cellidentifier)[g4hit_i];
+ one_g4hit.sampling = (*m_g4hit_sampling)[g4hit_i];
+ one_g4hit.hit_time = (*m_g4hit_time)[g4hit_i];
+ //one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i];
+ //scale the hit energy with the sampling fraction
+ if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20)
+ { //tile
+ //std::cout <<"Tile: "<<(*m_g4hit_energy)[g4hit_i]<<" "<<(*m_g4hit_samplingfraction)[g4hit_i]<<std::endl;
+ if ((*m_g4hit_samplingfraction)[g4hit_i])
+ {
+ one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] * (*m_g4hit_samplingfraction)[g4hit_i];
+ //std::cout <<"TileE: "<<one_g4hit.hit_energy<<std::endl;
+ }
+ else one_g4hit.hit_energy = 0.;
+ }
+ else
+ {
+ //std::cout <<"LAr: "<<(*m_g4hit_energy)[g4hit_i]<<" "<<(*m_g4hit_samplingfraction)[g4hit_i]<<std::endl;
+ one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] / (*m_g4hit_samplingfraction)[g4hit_i];
+ }
+ //one_g4hit.hit_sampfrac = (*m_g4hit_samplingfraction)[g4hit_i];
+ //new g4hit -> insert vector with 1 element
+ g4hits.insert(std::pair<Long64_t, std::vector<FCS_g4hit> >(one_g4hit.cell_identifier, std::vector<FCS_g4hit>(1, one_g4hit)));
+ }
+ else
+ {
+ //G4 hit exists in this identifier -> push_back new to the vector //FCS_g4hit one_g4hit;
+ one_g4hit.identifier = (*m_g4hit_identifier)[g4hit_i];
+ one_g4hit.cell_identifier = (*m_g4hit_cellidentifier)[g4hit_i];
+ one_g4hit.sampling = (*m_g4hit_sampling)[g4hit_i];
+ one_g4hit.hit_time = (*m_g4hit_time)[g4hit_i];
+ if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20)
+ { //tile
+ //std::cout <<"Tile2: "<<(*m_g4hit_energy)[g4hit_i]<<" "<<(*m_g4hit_samplingfraction)[g4hit_i]<<std::endl;
+ if ((*m_g4hit_samplingfraction)[g4hit_i])
+ {
+ one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] * (*m_g4hit_samplingfraction)[g4hit_i];
+ //std::cout <<"TileE2: "<<one_g4hit.hit_energy<<std::endl;
+ }
+ else one_g4hit.hit_energy = 0.;
+ }
+ else
+ {
+ //std::cout <<"LAr2: "<<(*m_g4hit_energy)[g4hit_i]<<" "<<(*m_g4hit_samplingfraction)[g4hit_i]<<std::endl;
+ one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] / (*m_g4hit_samplingfraction)[g4hit_i];
+ }
+ //one_g4hit.hit_sampfrac = (*m_g4hit_samplingfraction)[g4hit_i];
+ g4hits[(*m_g4hit_cellidentifier)[g4hit_i]].push_back(one_g4hit);
+ }
+ }
+ }
+
+ //hits
+ for (unsigned int hit_i = 0; hit_i < m_hit_identifier->size(); hit_i++)
+ {
+ if ((*m_hit_sampling)[hit_i] >= 0 && (*m_hit_sampling)[hit_i] <= 25 && (*m_hit_time)[hit_i] > m_TimingCut)
+ {
+ // if (m_Debug > 1)
+ ATH_MSG_DEBUG("Ignoring FCS hit, time too large: " << hit_i << " time: " << (*m_hit_time)[hit_i]);
+ continue;
+ }
+ if (hits.find((*m_hit_cellidentifier)[hit_i]) == hits.end())
+ {
+ //Detailed hit doesn't exist yet
+ one_hit.identifier = (*m_hit_identifier)[hit_i];
+ one_hit.cell_identifier = (*m_hit_cellidentifier)[hit_i];
+ one_hit.sampling = (*m_hit_sampling)[hit_i];
+
+ if (one_hit.sampling >= 12 && one_hit.sampling <= 20)
+ { //tile
+ if ((*m_hit_samplingfraction)[hit_i])
+ {
+ one_hit.hit_energy = (*m_hit_energy)[hit_i] * (*m_hit_samplingfraction)[hit_i];
+ }
+ else one_hit.hit_energy = 0.;
+ }
+ else
+ {
+ one_hit.hit_energy = (*m_hit_energy)[hit_i] / (*m_hit_samplingfraction)[hit_i];
+ }
+ //one_hit.hit_sampfrac = (*m_hit_samplingfraction)[hit_i];
+ one_hit.hit_time = (*m_hit_time)[hit_i];
+ one_hit.hit_x = (*m_hit_x)[hit_i];
+ one_hit.hit_y = (*m_hit_y)[hit_i];
+ one_hit.hit_z = (*m_hit_z)[hit_i];
+ hits.insert(std::pair<Long64_t, std::vector<FCS_hit> >(one_hit.cell_identifier, std::vector<FCS_hit>(1, one_hit)));
+ }
+ else
+ {
+ //Detailed hit exists in this identifier -> push_back new to the vector
+ one_hit.identifier = (*m_hit_identifier)[hit_i];
+ one_hit.cell_identifier = (*m_hit_cellidentifier)[hit_i];
+ one_hit.sampling = (*m_hit_sampling)[hit_i];
+ //one_hit.hit_energy = (*m_hit_energy)[hit_i];
+ if (one_hit.sampling >= 12 && one_hit.sampling <= 20)
+ { //tile
+ if ((*m_hit_samplingfraction)[hit_i])
+ {
+ one_hit.hit_energy = (*m_hit_energy)[hit_i] * (*m_hit_samplingfraction)[hit_i];
+ }
+ else one_hit.hit_energy = 0.;
+ }
+ else
+ {
+ one_hit.hit_energy = (*m_hit_energy)[hit_i] / (*m_hit_samplingfraction)[hit_i];
+ }
+ //one_hit.hit_sampfrac = (*m_hit_samplingfraction)[hit_i];
+ one_hit.hit_time = (*m_hit_time)[hit_i];
+ one_hit.hit_x = (*m_hit_x)[hit_i];
+ one_hit.hit_y = (*m_hit_y)[hit_i];
+ one_hit.hit_z = (*m_hit_z)[hit_i];
+ hits[(*m_hit_cellidentifier)[hit_i]].push_back(one_hit);
+ }
+ }
+
+ //Start matching:
+ Int_t iindex = 0;
+ for (std::map<Long64_t, FCS_cell>::iterator it = cells.begin(); it != cells.end(); )
+ {
+ iindex++;
+ // std::cout <<iindex<<std::endl;
+ one_matchedcell.clear(); //maybe not completely necessery, as we're not pushing_back into vectors
+ //set the cell part
+ one_matchedcell.cell = it->second;
+ //now look for FCS detailed hits in this cell
+ std::map<Long64_t, std::vector<FCS_hit> >::iterator it2 = hits.find(it->first);
+ if (it2 != hits.end())
+ {
+ //std::cout <<"FCS hits found in this cell"<<std::endl;
+ one_matchedcell.hit = it2->second;
+ hits.erase(it2); //remove it
+ }
+ else
+ {
+ //no hit found for this cell
+ one_matchedcell.hit.clear(); //important!
+ }
+ //now look for G4hits in this cell
+ std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first);
+ if (it3 != g4hits.end())
+ {
+ //std::cout <<"G4 hits found in this cell"<<std::endl;
+ one_matchedcell.g4hit = it3->second;
+ g4hits.erase(it3);
+ }
+ else
+ {
+ //no g4hit found for this cell
+ one_matchedcell.g4hit.clear();//important!
+ }
+ //std::cout <<"Erase cell"<<std::endl;
+ cells.erase(it++);
+ //std::cout <<"Insert matched object"<<std::endl;
+ //push_back matched cell for event jentry
+ m_oneeventcells->push_back(one_matchedcell);
+ //std::cout <<"Go next"<<std::endl;
+
+ }
+
+ //ok, cells should be empty, what about hits and g4hits?
+ //There could be G4hits/FCS hits for which we don't have a cell ->create a dummy empty cell with 0 energy, take the cell identifier from the hit
+ ATH_MSG_DEBUG("ISF_HitAnalysis Check after cells: " << cells.size() << " " << g4hits.size() << " " << hits.size());
+
+ for (std::map<Long64_t, std::vector<FCS_hit> >::iterator it = hits.begin(); it != hits.end();)
+ {
+ one_matchedcell.clear();
+ one_matchedcell.cell.cell_identifier = it->first;
+ //std::cout <<"This hit didn't exist in cell: "<<it->first<<std::endl;
+ if (it->second.size())
+ {
+ one_matchedcell.cell.sampling = (it->second)[0].sampling;
+ }
+ else
+ {
+ one_matchedcell.cell.sampling = -1; //
+ //ok, but you really shouldn't be here
+ ATH_MSG_DEBUG("ERROR: You shouldn't really be here");
+ }
+ one_matchedcell.cell.energy = 0.;
+ one_matchedcell.cell.center_x = 0.0;
+ one_matchedcell.cell.center_y = 0.0;
+ one_matchedcell.cell.center_z = 0.0;
+ one_matchedcell.hit = it->second;
+ std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first);
+ if (it3 != g4hits.end())
+ {
+ one_matchedcell.g4hit = it3->second;
+ g4hits.erase(it3);
+ }
+ else
+ {
+ //no g4hit found for this cell
+ one_matchedcell.g4hit.clear(); //important!
+ }
+ hits.erase(it++);
+ m_oneeventcells->push_back(one_matchedcell);
+
+ }
+
+ //ok, hits should be empty, what about g4hits?
+ ATH_MSG_DEBUG("ISF_HitAnalysis Check after hits: " << cells.size() << " " << g4hits.size() << " " << hits.size());
+ for (std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it = g4hits.begin(); it != g4hits.end();)
+ {
+ one_matchedcell.clear(); //maybe not so important
+ one_matchedcell.cell.cell_identifier = it->first;
+ if (it->second.size())
+ {
+ one_matchedcell.cell.sampling = (it->second)[0].sampling;
+ }
+ else
+ {
+ one_matchedcell.cell.sampling = -1; //
+ //not really
+ ATH_MSG_DEBUG("ERROR: You shouldn't really be here");
+ }
+ one_matchedcell.cell.energy = 0.;
+ one_matchedcell.cell.center_x = 0.0;
+ one_matchedcell.cell.center_y = 0.0;
+ one_matchedcell.cell.center_z = 0.0;
+ one_matchedcell.g4hit = it->second;
+ one_matchedcell.hit.clear(); //important!!
+ g4hits.erase(it++);
+ m_oneeventcells->push_back(one_matchedcell);
+ }
+
+ //Can fill the output tree already here:
+ m_total_cell_e = 0;
+ m_total_hit_e = 0;
+ m_total_g4hit_e = 0;
+
+ for (int j = 0; j < MAX_LAYER - 1; j++)
+ {
+ m_layercells[j]->m_vector = m_oneeventcells->GetLayer(j);
+ }
+
+ //this is for invalid cells
+ m_layercells[MAX_LAYER - 1]->m_vector = m_oneeventcells->GetLayer(-1);
+ for (int i = 0; i < MAX_LAYER; i++)
+ {
+ m_final_cell_energy->push_back(0.0); //zero for each event!
+ m_final_hit_energy->push_back(0.0);
+ m_final_g4hit_energy->push_back(0.0);
+
+ for (unsigned int cellindex = 0; cellindex < m_layercells[i]->size(); cellindex++)
+ {
+ if (i != MAX_LAYER - 1)
+ {
+ m_final_cell_energy->at(i) += m_layercells[i]->m_vector.at(cellindex).cell.energy;
+ m_total_cell_e += m_layercells[i]->m_vector.at(cellindex).cell.energy;
+ }
+ else
+ {
+ //don't add the energy in the invalid layer to the total energy (if there is any (shouldn't)
+ m_final_cell_energy->at(i) += m_layercells[i]->m_vector.at(cellindex).cell.energy; //this should be here anyway
+ }
+
+ //sum energy of all FCS detailed hits in this layer/cell
+ for (unsigned int j = 0; j < m_layercells[i]->m_vector.at(cellindex).hit.size(); j++)
+ {
+ if (i != MAX_LAYER - 1)
+ {
+ m_total_hit_e += m_layercells[i]->m_vector.at(cellindex).hit[j].hit_energy;
+ m_final_hit_energy->at(i) += m_layercells[i]->m_vector.at(cellindex).hit[j].hit_energy;
+ }
+ else
+ {
+ //again, don't add invalid layer energy to the sum
+ m_final_hit_energy->at(i) += m_layercells[i]->m_vector.at(cellindex).hit[j].hit_energy;
+ }
+ }
+
+ //sum energy of all G4 hits in this layer/cell
+ for (unsigned int j = 0; j < m_layercells[i]->m_vector.at(cellindex).g4hit.size(); j++)
+ {
+ if (i != MAX_LAYER - 1)
+ {
+ m_total_g4hit_e += m_layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy;
+ m_final_g4hit_energy->at(i) += m_layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy;
+ }
+ else
+ {
+ //don't add invalied layer energy to the sum
+ m_final_g4hit_energy->at(i) += m_layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy;
+ }
+ }
+ }
+ }
+
+ // push_back for total energy
+ m_final_cell_energy->push_back(0.0);
+ m_final_hit_energy->push_back(0.0);
+ m_final_g4hit_energy->push_back(0.0);
+
+ m_final_cell_energy->at(MAX_LAYER) = m_total_cell_e;
+ m_final_hit_energy->at(MAX_LAYER) = m_total_hit_e;
+ m_final_g4hit_energy->at(MAX_LAYER) = m_total_g4hit_e;
+
+
+
+
//Fill the tree and finish
if (m_tree) m_tree->Fill();
@@ -1348,609 +1631,6 @@ std::vector<Trk::HitInfo>* ISF_HitAnalysis::caloHits(const HepMC::GenParticle& p
return hitVector;
} //caloHits
-/*
-//#######################################################################
-void ISF_HitAnalysis::extrapolate(const HepMC::GenParticle* part,std::vector<Trk::HitInfo>* hitVector)
-{
- ATH_MSG_DEBUG("Start extrapolate()");
-
- m_ptruth_eta=part->momentum().eta();
- m_ptruth_phi=part->momentum().phi();
- m_ptruth_e =part->momentum().e();
- m_ptruth_pt =part->momentum().perp();
- m_ptruth_p =part->momentum().rho();
- m_pdgid=part->pdg_id();
-
- //////////////////////////////////////
- // Start calo extrapolation
- // First: get entry point into first calo sample
- //////////////////////////////////////
-
- get_calo_surface(hitVector);
-
- std::vector< std::vector<double> > eta_safe(3);
- std::vector< std::vector<double> > phi_safe(3);
- std::vector< std::vector<double> > r_safe(3);
- std::vector< std::vector<double> > z_safe(3);
- for(int subpos=CaloSubPos::SUBPOS_MID;subpos<=CaloSubPos::SUBPOS_EXT;++subpos)
- {
- eta_safe[subpos].resize(CaloCell_ID_FCS::MaxSample,-999.0);
- phi_safe[subpos].resize(CaloCell_ID_FCS::MaxSample,-999.0);
- r_safe[subpos].resize(CaloCell_ID_FCS::MaxSample,-999.0);
- z_safe[subpos].resize(CaloCell_ID_FCS::MaxSample,-999.0);
- }
-
- // only continue if inside the calo
- if( fabs(m_eta_calo_surf)<6 )
- {
- // now try to extrpolate to all calo layers, that contain energy
- ATH_MSG_DEBUG("Calo position for particle id "<<m_pdgid<<", trutheta= " << m_ptruth_eta <<", truthphi= "<<m_ptruth_phi<<", truthp="<<m_ptruth_p<<", truthpt="<<m_ptruth_pt);
- for(int sample=CaloCell_ID_FCS::FirstSample;sample<CaloCell_ID_FCS::MaxSample;++sample)
- {
- for(int subpos=CaloSubPos::SUBPOS_MID;subpos<=CaloSubPos::SUBPOS_EXT;++subpos)
- {
- m_letaCalo[sample][subpos]=-12345;
- m_lphiCalo[sample][subpos]=-12345;
- m_lrCalo[sample][subpos]=-12345;
- m_lzCalo[sample][subpos]=-12345;
- m_layerCaloOK[sample][subpos]=0;
- //ATH_MSG_DEBUG("============= Getting Calo position for sample "<<sample<<" E/Etot="<<p.E_layer[sample]<<" ; in : eta= " << m_ptruth_eta);
-
- if(get_calo_etaphi(hitVector,sample,subpos))
- {
- ATH_MSG_DEBUG( "Result in sample "<<sample<<"."<<subpos<<": eta="<<m_letaCalo[sample][subpos]<<" phi="<<m_lphiCalo[sample][subpos]<<" r="<<m_lrCalo[sample][subpos]<<" z="<<m_lzCalo[sample][subpos]<<" (ok="<<m_layerCaloOK[sample][subpos]<<")");
- eta_safe[subpos][sample]=m_letaCalo[sample][subpos];
- phi_safe[subpos][sample]=m_lphiCalo[sample][subpos];
- r_safe[subpos][sample]=m_lrCalo[sample][subpos];
- z_safe[subpos][sample]=m_lzCalo[sample][subpos];
- }
- else
- ATH_MSG_DEBUG( "Extrapolation to sample "<<sample<<" failed (ok="<<m_layerCaloOK[sample][subpos]<<")");
- }
- }
- } //inside calo
-
- m_TTC_back_eta->push_back(eta_safe[CaloSubPos::SUBPOS_EXT]);
- m_TTC_back_phi->push_back(phi_safe[CaloSubPos::SUBPOS_EXT]);
- m_TTC_back_r->push_back(r_safe[CaloSubPos::SUBPOS_EXT]);
- m_TTC_back_z->push_back(z_safe[CaloSubPos::SUBPOS_EXT]);
- m_TTC_entrance_eta->push_back(eta_safe[CaloSubPos::SUBPOS_ENT]);
- m_TTC_entrance_phi->push_back(phi_safe[CaloSubPos::SUBPOS_ENT]);
- m_TTC_entrance_r->push_back(r_safe[CaloSubPos::SUBPOS_ENT]);
- m_TTC_entrance_z->push_back(z_safe[CaloSubPos::SUBPOS_ENT]);
-
- ATH_MSG_DEBUG("End extrapolate()");
-} //extrapolate
-*/
-
-/*
-void ISF_HitAnalysis::extrapolate_to_ID(const HepMC::GenParticle* part,std::vector<Trk::HitInfo>* hitVector)
-{
- ATH_MSG_DEBUG("Start extrapolate_to_ID()");
-
- m_ptruth_eta=part->momentum().eta();
- m_ptruth_phi=part->momentum().phi();
- m_ptruth_e =part->momentum().e();
- m_ptruth_pt =part->momentum().perp();
- m_ptruth_p =part->momentum().rho();
- m_pdgid =part->pdg_id();
-
- Amg::Vector3D hitpos(0,0,0);
- Amg::Vector3D hitmom(0,0,0);
- if(rz_cylinder_get_calo_etaphi(hitVector,m_CaloBoundaryR,m_CaloBoundaryZ,hitpos,hitmom))
- {
- ATH_MSG_DEBUG("BOUNDARY ID-CALO eta="<<hitpos.eta()<<" phi="<<hitpos.phi()<<" r="<<hitpos.perp()<<" z="<<hitpos[Amg::z]<<" theta="<<hitpos.theta()<<" ; momentum eta="<<hitmom.eta()<<" phi="<<hitmom.phi()<<" theta="<<hitmom.theta());
- m_TTC_IDCaloBoundary_eta->push_back(hitpos.eta());
- m_TTC_IDCaloBoundary_phi->push_back(hitpos.phi());
- m_TTC_IDCaloBoundary_r->push_back(hitpos.perp());
- m_TTC_IDCaloBoundary_z->push_back(hitpos[Amg::z]);
- }
- else
- {
- ATH_MSG_DEBUG("Extrapolation to IDCaloBoundary failed");
- m_TTC_IDCaloBoundary_eta->push_back(-999.);
- m_TTC_IDCaloBoundary_phi->push_back(-999.);
- m_TTC_IDCaloBoundary_r->push_back(-999.);
- m_TTC_IDCaloBoundary_z->push_back(-999.);
- }
-
- TVector3 vec(hitpos[Amg::x],hitpos[Amg::y],hitpos[Amg::z]);
-
- //get the tangentvector on this interaction point:
- //GlobalMomentum* mom=params_on_surface_ID->TrackParameters::momentum().unit() ;
- //Trk::GlobalMomentum* trackmom=params_on_surface_ID->Trk::TrackParameters::momentum();
-
- if(hitmom.mag()>0)
- {
- //angle between vec and trackmom:
- TVector3 Trackmom(hitmom[Amg::x],hitmom[Amg::y],hitmom[Amg::z]);
- double angle3D=Trackmom.Angle(vec); //isn't this the same as TVector3 vec?
- ATH_MSG_DEBUG(" 3D ANGLE "<<angle3D);
-
- double angleEta=vec.Theta()-Trackmom.Theta();
- ATH_MSG_DEBUG(" ANGLE dTHEA"<<angleEta);
-
- m_TTC_Angle3D->push_back(angle3D);
- m_TTC_AngleEta->push_back(angleEta);
- }
- else
- {
- m_TTC_Angle3D->push_back(-999.);
- m_TTC_AngleEta->push_back(-999.);
- }
- ATH_MSG_DEBUG("End extrapolate_to_ID()");
-} //extrapolate_to_ID
-*/
-
-//UPDATED
-bool ISF_HitAnalysis::get_calo_etaphi(std::vector<Trk::HitInfo>* hitVector, int sample,int subpos)
-{
- m_layerCaloOK[sample][subpos]=false;
- m_letaCalo[sample][subpos]=m_eta_calo_surf;
- m_lphiCalo[sample][subpos]=m_phi_calo_surf;
- m_lrCalo[sample][subpos] = rpos(sample,m_eta_calo_surf,subpos);
- m_lzCalo[sample][subpos] = zpos(sample,m_eta_calo_surf,subpos);
- double distsamp =deta(sample,m_eta_calo_surf);
- double lrzpos =rzpos(sample,m_eta_calo_surf,subpos);
- //double hitdist=0;
- bool best_found=false;
- double best_target=0;
-
- std::vector<Trk::HitInfo>::iterator it = hitVector->begin();
- while ( it!= hitVector->end() && it->detID != (3000+sample) ) { it++;}
- //while ((*it).detID != (3000+sample) && it < hitVector->end() ) it++;
-
- if (it!=hitVector->end()) {
- Amg::Vector3D hitPos1 = (*it).trackParms->position();
- int sid1=(*it).detID;
- int sid2=-1;
- Amg::Vector3D hitPos;
- Amg::Vector3D hitPos2;
-
- std::vector<Trk::HitInfo>::iterator itnext = it;
- ++itnext;
- if(itnext!=hitVector->end()) {
- hitPos2 = (*itnext).trackParms->position();
- sid2=(*itnext).detID;
- double eta_avg=0.5*(hitPos1.eta()+hitPos2.eta());
- double t;
- if(isCaloBarrel(sample)) {
- double r=rpos(sample,eta_avg,subpos);
- double r1=hitPos1.perp();
- double r2=hitPos2.perp();
- t=(r-r1)/(r2-r1);
- best_target=r;
- } else {
- double z=zpos(sample,eta_avg,subpos);
- double z1=hitPos1[Amg::z];
- double z2=hitPos2[Amg::z];
- t=(z-z1)/(z2-z1);
- best_target=z;
- }
- hitPos=t*hitPos2+(1-t)*hitPos1;
- } else {
- hitPos=hitPos1;
- hitPos2=hitPos1;
- }
- double etaCalo = hitPos.eta();
- double phiCalo = hitPos.phi();
- m_letaCalo[sample][subpos]=etaCalo;
- m_lphiCalo[sample][subpos]=phiCalo;
- m_lrCalo[sample][subpos] = hitPos.perp();
- m_lzCalo[sample][subpos] = hitPos[Amg::z];
- //hitdist=hitPos.mag();
- m_layerCaloOK[sample][subpos]=true;
- lrzpos=rzpos(sample,etaCalo,subpos);
- distsamp=deta(sample,etaCalo);
- best_found=true;
- ATH_MSG_DEBUG(" extrapol with layer hit: id="<<sid1<<" -> "<<sid2<<
- " target r/z="<<best_target<<
- " r1="<<hitPos1.perp()<<" z1="<<hitPos1[Amg::z]<<" r2="<<hitPos2.perp()<<" z2="<<hitPos2[Amg::z]<<
- " re="<<m_lrCalo[sample][subpos]<<" ze="<<m_lzCalo[sample][subpos]
- );
- }
- if(!best_found) {
- it = hitVector->begin();
- double best_dist=0.5;
- bool best_inside=false;
- int best_id1,best_id2;
- Amg::Vector3D best_hitPos=(*it).trackParms->position();
- while (it < hitVector->end()-1 ) {
- Amg::Vector3D hitPos1 = (*it).trackParms->position();
- int sid1=(*it).detID;
- it++;
- Amg::Vector3D hitPos2 = (*it).trackParms->position();
- int sid2=(*it).detID;
- double eta_avg=0.5*(hitPos1.eta()+hitPos2.eta());
- double t;
- double tmp_target=0;
- if(isCaloBarrel(sample)) {
- double r=rpos(sample,eta_avg,subpos);
- double r1=hitPos1.perp();
- double r2=hitPos2.perp();
- t=(r-r1)/(r2-r1);
- tmp_target=r;
- } else {
- double z=zpos(sample,eta_avg,subpos);
- double z1=hitPos1[Amg::z];
- double z2=hitPos2[Amg::z];
- t=(z-z1)/(z2-z1);
- tmp_target=z;
- }
- Amg::Vector3D hitPos=t*hitPos2+(1-t)*hitPos1;
- double dist=TMath::Min(TMath::Abs(t-0),TMath::Abs(t-1));
- bool inside=false;
- if(t>=0 && t<=1) inside=true;
- if(!best_found || inside) {
- if(!best_inside || dist<best_dist) {
- best_dist=dist;
- best_hitPos=hitPos;
- best_inside=inside;
- best_found=true;
- best_id1=sid1;
- best_id2=sid2;
- best_target=tmp_target;
- }
- } else {
- if(!best_inside && dist<best_dist) {
- best_dist=dist;
- best_hitPos=hitPos;
- best_inside=inside;
- best_found=true;
- best_id1=sid1;
- best_id2=sid2;
- best_target=tmp_target;
- }
- }
- ATH_MSG_VERBOSE(" extrapol without layer hit: id="<<sid1<<" -> "<<sid2<<" dist="<<dist<<" mindist="<<best_dist<<
- " t="<<t<<" best_inside="<<best_inside<<" target r/z="<<tmp_target<<
- " r1="<<hitPos1.perp()<<" z1="<<hitPos1[Amg::z]<<" r2="<<hitPos2.perp()<<" z2="<<hitPos2[Amg::z]<<
- " re="<<hitPos.perp()<<" ze="<<hitPos[Amg::z]<<
- " rb="<<best_hitPos.perp()<<" zb="<<best_hitPos[Amg::z]
- );
- if(best_found) {
- m_letaCalo[sample][subpos]=best_hitPos.eta();
- m_lphiCalo[sample][subpos]=best_hitPos.phi();
- m_lrCalo[sample][subpos] = best_hitPos.perp();
- m_lzCalo[sample][subpos] = best_hitPos[Amg::z];
- //hitdist=best_hitPos.mag();
- lrzpos=rzpos(sample,m_letaCalo[sample][subpos],subpos);
- distsamp=deta(sample,m_letaCalo[sample][subpos]);
- m_layerCaloOK[sample][subpos]=true;
- }
- }
- if(best_found) {
- ATH_MSG_DEBUG(" extrapol without layer hit: id="<<best_id1<<" -> "<<best_id2<<" mindist="<<best_dist<<
- " best_inside="<<best_inside<<" target r/z="<<best_target<<
- " rb="<<best_hitPos.perp()<<" zb="<<best_hitPos[Amg::z]
- );
- }
- }
-
- if(isCaloBarrel(sample)) lrzpos*=cosh(m_letaCalo[sample][subpos]);
- else lrzpos= fabs(lrzpos/tanh(m_letaCalo[sample][subpos]));
-
- m_dCalo[sample][subpos]=lrzpos;
- m_distetaCaloBorder[sample][subpos]=distsamp;
-
- return m_layerCaloOK[sample][subpos];
-} //get_calo_etaphi
-
-//UPDATED
-bool ISF_HitAnalysis::rz_cylinder_get_calo_etaphi(std::vector<Trk::HitInfo>* hitVector, double cylR, double cylZ, Amg::Vector3D& pos, Amg::Vector3D& mom)
-{
- bool best_found=false;
- double best_dist=10000;
- bool best_inside=false;
- int best_id1,best_id2;
-
- std::vector<Trk::HitInfo>::iterator it = hitVector->begin();
- Amg::Vector3D best_hitPos=(*it).trackParms->position();
- for(int rz=0;rz<=1;++rz) {
- it = hitVector->begin();
- while (it < hitVector->end()-1 ) {
- Amg::Vector3D hitPos1 = (*it).trackParms->position();
- Amg::Vector3D hitMom1 = (*it).trackParms->momentum();
- int sid1=(*it).detID;
- it++;
- Amg::Vector3D hitPos2 = (*it).trackParms->position();
- Amg::Vector3D hitMom2 = (*it).trackParms->momentum();
- int sid2=(*it).detID;
-
- double t;
- if(rz==1) {
- double r=cylR;
- double r1=hitPos1.perp();
- double r2=hitPos2.perp();
- t=(r-r1)/(r2-r1);
- } else {
- double z=cylZ;
- double z1=hitPos1[Amg::z];
- double z2=hitPos2[Amg::z];
- t=(z-z1)/(z2-z1);
- }
- Amg::Vector3D hitPos=t*hitPos2+(1-t)*hitPos1;
-
- double dist=hitPos.mag();
- bool inside=false;
- if(t>=-0.001 && t<=1.001) inside=true;
-
- if(!best_found || inside) {
- if(!best_inside || dist<best_dist) {
- best_dist=dist;
- best_hitPos=hitPos;
- best_inside=inside;
- best_found=true;
- best_id1=sid1;
- best_id2=sid2;
- mom=t*hitMom2+(1-t)*hitMom1;
- }
- } else {
- if(!best_inside && dist<best_dist) {
- best_dist=dist;
- best_hitPos=hitPos;
- best_inside=inside;
- best_found=true;
- best_id1=sid1;
- best_id2=sid2;
- mom=t*hitMom2+(1-t)*hitMom1;
- }
- }
- ATH_MSG_DEBUG(" extrapol without layer hit: id="<<sid1<<" -> "<<sid2<<" dist="<<dist<<" bestdist="<<best_dist<<
- " t="<<t<<" inside="<<inside<<" best_inside="<<best_inside<<
- " r1="<<hitPos1.perp()<<" z1="<<hitPos1[Amg::z]<<" r2="<<hitPos2.perp()<<" z2="<<hitPos2[Amg::z]<<
- " re="<<hitPos.perp()<<" ze="<<hitPos[Amg::z]<<
- " rb="<<best_hitPos.perp()<<" zb="<<best_hitPos[Amg::z]
- );
- }
- }
-
- if(best_found) {
- ATH_MSG_DEBUG(" extrapol to r="<<cylR<<" z="<<cylZ<<": id="<<best_id1<<" -> "<<best_id2<<" dist="<<best_dist<<
- " best_inside="<<best_inside<<
- " rb="<<best_hitPos.perp()<<" zb="<<best_hitPos[Amg::z]
- );
- }
- pos=best_hitPos;
-
- return best_found;
-} //rz_cylinder_get_calo_etaphi
-
-//UPDATED
-/*
-bool ISF_HitAnalysis::get_calo_etaphi(std::vector<Trk::HitInfo>* hitVector, CaloCell_ID_FCS::CaloSample sample)
-{
- m_layerCaloOK[sample]=false;
- m_letaCalo[sample]=m_eta_calo_surf;
- m_lphiCalo[sample]=m_phi_calo_surf;
- double distsamp =deta(sample,m_eta_calo_surf);
- double rzmiddle =rzmid(sample,m_eta_calo_surf);
- double hitdist=0;
- bool best_found=false;
- double best_target=0;
-
- std::vector<Trk::HitInfo>::iterator it = hitVector->begin();
- while ( it!= hitVector->end() && it->detID != (3000+sample) ) { it++;}
- //while ((*it).detID != (3000+sample) && it < hitVector->end() ) it++;
-
- if (it!=hitVector->end()) {
- Amg::Vector3D hitPos1 = (*it).trackParms->position();
- int sid1=(*it).detID;
- int sid2=-1;
- Amg::Vector3D hitPos;
- Amg::Vector3D hitPos2;
-
- std::vector<Trk::HitInfo>::iterator itnext = it;
- ++itnext;
- if(itnext!=hitVector->end()) {
- hitPos2 = (*itnext).trackParms->position();
- sid2=(*itnext).detID;
- double eta_avg=0.5*(hitPos1.eta()+hitPos2.eta());
- double t;
- if(isCaloBarrel(sample)) {
- double r=rmid(sample,eta_avg);
- double r1=hitPos1.perp();
- double r2=hitPos2.perp();
- t=(r-r1)/(r2-r1);
- best_target=r;
- } else {
- double z=zmid(sample,eta_avg);
- double z1=hitPos1[Amg::z];
- double z2=hitPos2[Amg::z];
- t=(z-z1)/(z2-z1);
- best_target=z;
- }
- hitPos=t*hitPos2+(1-t)*hitPos1;
- } else {
- hitPos=hitPos1;
- hitPos2=hitPos1;
- }
- double etaCalo = hitPos.eta();
- double phiCalo = hitPos.phi();
- m_letaCalo[sample]=etaCalo;
- m_lphiCalo[sample]=phiCalo;
- hitdist=hitPos.mag();
- m_layerCaloOK[sample]=true;
- rzmiddle=rzmid((CaloCell_ID_FCS::CaloSample)sample,etaCalo);
- distsamp=deta((CaloCell_ID_FCS::CaloSample)sample,etaCalo);
- best_found=true;
- ATH_MSG_DEBUG(" extrapol with layer hit: id="<<sid1<<" -> "<<sid2<<
- " target r/z="<<best_target<<
- " r1="<<hitPos1.perp()<<" z1="<<hitPos1[Amg::z]<<" r2="<<hitPos2.perp()<<" z2="<<hitPos2[Amg::z]<<
- " re="<<hitPos.perp()<<" ze="<<hitPos[Amg::z]
- );
- }
- if(!best_found) {
- it = hitVector->begin();
- double best_dist=0.5;
- bool best_inside=false;
- int best_id1,best_id2;
- Amg::Vector3D best_hitPos=(*it).trackParms->position();
- while (it < hitVector->end()-1 ) {
- Amg::Vector3D hitPos1 = (*it).trackParms->position();
- int sid1=(*it).detID;
- it++;
- Amg::Vector3D hitPos2 = (*it).trackParms->position();
- int sid2=(*it).detID;
- double eta_avg=0.5*(hitPos1.eta()+hitPos2.eta());
- double t;
- double tmp_target=0;
- if(isCaloBarrel(sample)) {
- double r=rmid(sample,eta_avg);
- double r1=hitPos1.perp();
- double r2=hitPos2.perp();
- t=(r-r1)/(r2-r1);
- tmp_target=r;
- } else {
- double z=zmid(sample,eta_avg);
- double z1=hitPos1[Amg::z];
- double z2=hitPos2[Amg::z];
- t=(z-z1)/(z2-z1);
- tmp_target=z;
- }
- Amg::Vector3D hitPos=t*hitPos2+(1-t)*hitPos1;
- double dist=TMath::Min(TMath::Abs(t-0),TMath::Abs(t-1));
- bool inside=false;
- if(t>=0 && t<=1) inside=true;
- if(!best_found || inside) {
- if(!best_inside || dist<best_dist) {
- best_dist=dist;
- best_hitPos=hitPos;
- best_inside=inside;
- best_found=true;
- best_id1=sid1;
- best_id2=sid2;
- best_target=tmp_target;
- }
- } else {
- if(!best_inside && dist<best_dist) {
- best_dist=dist;
- best_hitPos=hitPos;
- best_inside=inside;
- best_found=true;
- best_id1=sid1;
- best_id2=sid2;
- best_target=tmp_target;
- }
- }
- ATH_MSG_VERBOSE(" extrapol without layer hit: id="<<sid1<<" -> "<<sid2<<" dist="<<dist<<" mindist="<<best_dist<<
- " t="<<t<<" best_inside="<<best_inside<<" target r/z="<<tmp_target<<
- " r1="<<hitPos1.perp()<<" z1="<<hitPos1[Amg::z]<<" r2="<<hitPos2.perp()<<" z2="<<hitPos2[Amg::z]<<
- " re="<<hitPos.perp()<<" ze="<<hitPos[Amg::z]<<
- " rb="<<best_hitPos.perp()<<" zb="<<best_hitPos[Amg::z]
- );
- if(best_found) {
- m_letaCalo[sample]=best_hitPos.eta();
- m_lphiCalo[sample]=best_hitPos.phi();
- hitdist=best_hitPos.mag();
- rzmiddle=rzmid((CaloCell_ID_FCS::CaloSample)sample,m_letaCalo[sample]);
- distsamp=deta((CaloCell_ID_FCS::CaloSample)sample,m_letaCalo[sample]);
- m_layerCaloOK[sample]=true;
- }
- }
- if(best_found) {
- ATH_MSG_DEBUG(" extrapol without layer hit: id="<<best_id1<<" -> "<<best_id2<<" mindist="<<best_dist<<
- " best_inside="<<best_inside<<" target r/z="<<best_target<<
- " rb="<<best_hitPos.perp()<<" zb="<<best_hitPos[Amg::z]
- );
- }
- }
-
- if(isCaloBarrel((CaloCell_ID_FCS::CaloSample)sample)) rzmiddle*=cosh(m_letaCalo[sample]);
- else rzmiddle= fabs(rzmiddle/tanh(m_letaCalo[sample]));
-
- m_dCalo[sample]=rzmiddle;
- m_distetaCaloBorder[sample]=distsamp;
-
- if(msgLvl(MSG::DEBUG)) {
- msg(MSG::DEBUG)<<" Final par TTC sample "<<(int)sample;
- if(m_layerCaloOK[sample]) msg()<<" (good)";
- else msg()<<" (bad)";
- msg()<<" eta=" << m_letaCalo[sample] << " phi=" << m_lphiCalo[sample] <<" m_dCalo="<<m_dCalo[sample]<<" dist(hit)="<<hitdist<< endmsg;
- }
-
- return m_layerCaloOK[sample];
-}
-*/
-
-//UPDATED
-bool ISF_HitAnalysis::get_calo_surface(std::vector<Trk::HitInfo>* hitVector)
-{
- ATH_MSG_DEBUG("Start get_calo_surface()");
- m_sample_calo_surf=CaloCell_ID_FCS::noSample;
- m_eta_calo_surf=-999;
- m_phi_calo_surf=-999;
- m_d_calo_surf=0;
- double min_calo_surf_dist=1000;
-
- for(unsigned int i=0;i<m_surfacelist.size();++i) {
- CaloCell_ID_FCS::CaloSample sample=m_surfacelist[i];
- std::vector<Trk::HitInfo>::iterator it = hitVector->begin();
- while (it != hitVector->end() && it->detID != (3000+sample) ) { it++;}
- if(it==hitVector->end()) continue;
- Amg::Vector3D hitPos = (*it).trackParms->position();
-
- //double offset = 0.;
- double etaCalo = hitPos.eta();
-
- ATH_MSG_DEBUG("test: entrance to calo surface : sample="<<sample<<" eta="<<etaCalo);
-
- if (fabs(etaCalo)<900) {
- double phiCalo = hitPos.phi();
- double distsamp =deta(sample,etaCalo);
-
- ATH_MSG_DEBUG(" phi="<<phiCalo<<" dist="<<distsamp);
- if(distsamp<min_calo_surf_dist && min_calo_surf_dist>=0)
- {
- //hitVector is ordered in r, so if first surface was hit, keep it
- m_eta_calo_surf=etaCalo;
- m_phi_calo_surf=phiCalo;
- min_calo_surf_dist=distsamp;
- m_sample_calo_surf=sample;
- m_d_calo_surf=rzent(sample,etaCalo);
- ATH_MSG_DEBUG(" r/z="<<m_d_calo_surf);
- if(isCaloBarrel(sample)) m_d_calo_surf*=cosh(etaCalo);
- else m_d_calo_surf= fabs(m_d_calo_surf/tanh(etaCalo));
- ATH_MSG_DEBUG(" d="<<m_d_calo_surf);
- if(distsamp<0)
- {
- break;
- }
- }
- }
- else
- {
- ATH_MSG_DEBUG(": eta > 900, not using this");
- }
- }
-
- if(m_sample_calo_surf==CaloCell_ID_FCS::noSample) {
- // first intersection with sensitive calo layer
- std::vector<Trk::HitInfo>::iterator it = hitVector->begin();
- while ( it < hitVector->end() && (*it).detID != 3 ) { it++;} // to be updated
- if (it==hitVector->end()) { // no calo intersection, abort
- return false;
- }
- Amg::Vector3D surface_hitPos = (*it).trackParms->position();
- m_eta_calo_surf=surface_hitPos.eta();
- m_phi_calo_surf=surface_hitPos.phi();
- m_d_calo_surf=surface_hitPos.mag();
-
- double pT=(*it).trackParms->momentum().perp();
- if(TMath::Abs(m_eta_calo_surf)>4.9 || pT<500 || (TMath::Abs(m_eta_calo_surf)>4 && pT<1000) ) {
- ATH_MSG_DEBUG("only entrance to calo entrance layer found, no surface : eta="<<m_eta_calo_surf<<" phi="<<m_phi_calo_surf<<" d="<<m_d_calo_surf<<" pT="<<pT);
- } else {
- ATH_MSG_WARNING("only entrance to calo entrance layer found, no surface : eta="<<m_eta_calo_surf<<" phi="<<m_phi_calo_surf<<" d="<<m_d_calo_surf<<" pT="<<pT);
- }
- } else {
- ATH_MSG_DEBUG("entrance to calo surface : sample="<<m_sample_calo_surf<<" eta="<<m_eta_calo_surf<<" phi="<<m_phi_calo_surf<<" deta="<<min_calo_surf_dist<<" dsurf="<<m_d_calo_surf);
- }
-
- ATH_MSG_DEBUG("End get_calo_surface()");
- return true;
-}
-
-
bool ISF_HitAnalysis::isCaloBarrel(int sample) const
{
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.C
deleted file mode 100644
index eb6473568b009f67b68286c7cab8078dc609802b..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.C
+++ /dev/null
@@ -1,529 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#define CaloHitAna_cxx
-#include "CaloHitAna.h"
-#include <TH2.h>
-#include <TStyle.h>
-#include <TCanvas.h>
-#include "FCS_Cell.h"
-#include <iostream>
-#include <map>
-#include "TSystem.h"
-
-using namespace std;
-
-void CaloHitAna::Loop()
-{
- //This Loop is supposed to read the output ntuple/tree
- //(vector of cells, vector of G4 hits, vector of FCS detailed hits, ...)
- //from the athena ISF_HitAnalysis algorithm (==ESD dumper)
- //and match them together + save structured output
- //NB: E = Ehit * SamplingFraction for Tile or E = Ehit / SamplingFraction (LAr) conversion
- //is done here
-
- if (fChain == 0) return;
- ProcInfo_t procinfo;
- Long64_t nentries = fChain->GetEntriesFast();
-
- if (m_max_nentries >= 0 && m_max_nentries < nentries) nentries = m_max_nentries;
-
- std::map<Long64_t, FCS_cell> cells; //read all objects and collect them by identifier (Long64_t)
- std::map<Long64_t, std::vector<FCS_g4hit> > g4hits;
- std::map<Long64_t, std::vector<FCS_hit> > hits;
-
- FCS_cell one_cell; //note that this is not extra safe if I don't have a clear method!
- FCS_g4hit one_g4hit;
- FCS_hit one_hit;
- FCS_matchedcell one_matchedcell;
- //FCS_truth one_truth;
-
- //From here: Loop over events:
- Long64_t nbytes = 0, nb = 0;
- std::cout << "before event loop" << std::endl;
- for (Long64_t jentry = 0; jentry < nentries; jentry++)
- {
- if (jentry % 100 == 0)
- std::cout << "jentry " << jentry << endl;
-
- Long64_t ientry = LoadTree(jentry);
- if (ientry < 0) break;
- nb = fChain->GetEntry(jentry); nbytes += nb;
- if (jentry % m_PrintOutFrequency == 0)
- {
- std::cout << "Processing entry " << jentry << std::endl;
- gSystem->GetProcInfo(&procinfo);
- std::cout << "Memory usage: " << procinfo.fMemResident << " " << procinfo.fMemVirtual << std::endl;
- }
-
- //Copy truth
- new_truthE->clear();
- new_truthPx->clear();
- new_truthPy->clear();
- new_truthPz->clear();
- new_truthPDG->clear();
- new_truthBarcode->clear();
- new_truthVtxBarcode->clear();
- //truthcollection->clear();
- m_newTTC_entrance_eta->clear();
- m_newTTC_entrance_phi->clear();
- m_newTTC_entrance_r->clear();
- m_newTTC_entrance_z->clear();
- m_newTTC_back_eta->clear();
- m_newTTC_back_phi->clear();
- m_newTTC_back_r->clear();
- m_newTTC_back_z->clear();
- m_newTTC_mid_eta->clear();
- m_newTTC_mid_phi->clear();
- m_newTTC_mid_r->clear();
- m_newTTC_mid_z->clear();
- m_newTTC_IDCaloBoundary_eta->clear();
- m_newTTC_IDCaloBoundary_phi->clear();
- m_newTTC_IDCaloBoundary_r->clear();
- m_newTTC_IDCaloBoundary_z->clear();
- m_newTTC_Angle3D->clear();
- m_newTTC_AngleEta->clear();
-
- oneeventcells->m_vector.clear();
-
- if (TruthE->size() != 1)
- std::cout << "Strange! TruthE->size() is " << TruthE->size() << ", but it should be 1. Please check input." << std::endl;
-
- for (unsigned int truth_i = 0; truth_i < TruthE->size(); truth_i++)
- {
- //std::cout <<truth_i<<" "<<TruthE->size()<<" "<<TruthPDG->at(truth_i)<<std::endl;
- //std::cout <<"TTC2: "<<(*TTC_entrance_eta)[truth_i].size()<<std::endl;//this one has 24 layers ->ok
-
- new_truthE->push_back(TruthE->at(truth_i));
- new_truthPx->push_back(TruthPx->at(truth_i));
- new_truthPy->push_back(TruthPy->at(truth_i));
- new_truthPz->push_back(TruthPz->at(truth_i));
- new_truthPDG->push_back(TruthPDG->at(truth_i));
- new_truthBarcode->push_back(TruthBarcode->at(truth_i));
- //new_truthVtxBarcode->push_back(TruthVtxBarcode->at(truth_i));
-
- m_newTTC_IDCaloBoundary_eta->push_back(newTTC_IDCaloBoundary_eta->at(truth_i));
- m_newTTC_IDCaloBoundary_phi->push_back(newTTC_IDCaloBoundary_phi->at(truth_i));
- m_newTTC_IDCaloBoundary_r ->push_back(newTTC_IDCaloBoundary_r->at(truth_i));
- m_newTTC_IDCaloBoundary_z ->push_back(newTTC_IDCaloBoundary_z->at(truth_i));
- m_newTTC_Angle3D ->push_back(newTTC_Angle3D->at(truth_i));
- m_newTTC_AngleEta ->push_back(newTTC_AngleEta->at(truth_i));
-
- //create temporary vectors
- vector<float> new_entrance_eta;
- vector<float> new_entrance_phi;
- vector<float> new_entrance_r;
- vector<float> new_entrance_z;
- vector<float> new_back_eta;
- vector<float> new_back_phi;
- vector<float> new_back_r;
- vector<float> new_back_z;
- vector<float> new_mid_eta;
- vector<float> new_mid_phi;
- vector<float> new_mid_r;
- vector<float> new_mid_z;
-
- for (unsigned int s = 0; s < 24; s++)
- {
- new_entrance_eta.push_back((newTTC_entrance_eta->at(truth_i))[s]);
- new_entrance_phi.push_back((newTTC_entrance_phi->at(truth_i))[s]);
- new_entrance_r.push_back((newTTC_entrance_r->at(truth_i))[s]);
- new_entrance_z.push_back((newTTC_entrance_z->at(truth_i))[s]);
- new_back_eta.push_back((newTTC_back_eta->at(truth_i))[s]);
- new_back_phi.push_back((newTTC_back_phi->at(truth_i))[s]);
- new_back_r.push_back((newTTC_back_r->at(truth_i))[s]);
- new_back_z.push_back((newTTC_back_z->at(truth_i))[s]);
- new_mid_eta.push_back((newTTC_mid_eta->at(truth_i))[s]);
- new_mid_phi.push_back((newTTC_mid_phi->at(truth_i))[s]);
- new_mid_r.push_back((newTTC_mid_r->at(truth_i))[s]);
- new_mid_z.push_back((newTTC_mid_z->at(truth_i))[s]);
- }
-
- //push back temporary vectors
- m_newTTC_entrance_eta->push_back(new_entrance_eta);
- m_newTTC_entrance_phi->push_back(new_entrance_phi);
- m_newTTC_entrance_r->push_back(new_entrance_r);
- m_newTTC_entrance_z->push_back(new_entrance_z);
- m_newTTC_back_eta->push_back(new_back_eta);
- m_newTTC_back_phi->push_back(new_back_phi);
- m_newTTC_back_r->push_back(new_back_r);
- m_newTTC_back_z->push_back(new_back_z);
- m_newTTC_mid_eta->push_back(new_mid_eta);
- m_newTTC_mid_phi->push_back(new_mid_phi);
- m_newTTC_mid_r->push_back(new_mid_r);
- m_newTTC_mid_z->push_back(new_mid_z);
-
-
- /*
- one_truth.SetPxPyPzE((*TruthPx)[truth_i], (*TruthPy)[truth_i], (*TruthPz)[truth_i], (*TruthE)[truth_i]);
- one_truth.pdgid = TruthPDG->at(truth_i);
- one_truth.barcode = TruthBarcode->at(truth_i);
- one_truth.vtxbarcode = TruthVtxBarcode->at(truth_i);
- truthcollection->push_back(one_truth);
- */
-
- }
-
- //std::cout<<"check. after truth block"<<std::endl;
-
- //Now matching between cells, G4hits and G4detailed hits
- //sort cells by identifier:
- //clear first the containers for this event:
- cells.clear();
- g4hits.clear();
- hits.clear();
- //std::cout <<"Check Original size: Cells: "<<CellIdentifier->size()<<" G4Hits: "<<G4HitCellIdentifier->size()<<" FCSHits: "<<HitCellIdentifier->size()<<std::endl;
- if (m_Debug > 1) std::cout << "Reading cells...";
-
- for (unsigned int cell_i = 0; cell_i < CellIdentifier->size(); cell_i++)
- {
- if (cells.find((*CellIdentifier)[cell_i]) == cells.end()) //doesn't exist
- {
- one_cell.cell_identifier = (*CellIdentifier)[cell_i];
- one_cell.sampling = (*CellSampling)[cell_i];
- one_cell.energy = (*CellE)[cell_i];
- one_cell.center_x = 0.0; //for now
- one_cell.center_y = 0.0;
- one_cell.center_z = 0.0;
- cells.insert(std::pair<Long64_t, FCS_cell>(one_cell.cell_identifier, one_cell));
- }
- else
- {
- //there shouldn't be a cell with the same identifier in this event
- std::cout << "ISF_HitAnalysis: Same cell???? ERROR" << std::endl;
- }
- }
-
- if (m_Debug > 1) std::cout << " Done" << std::endl;
-
- if (m_Debug > 1) std::cout << "Reading G4hits";
-
- if (m_do_g4_hits)
- {
- for (unsigned int g4hit_i = 0; g4hit_i < G4HitIdentifier->size(); g4hit_i++)
- {
- if ((*G4HitSampling)[g4hit_i] >= 0 && (*G4HitSampling)[g4hit_i] <= 25 && (*G4HitT)[g4hit_i] > m_TimingCut)
- {
- if (m_Debug > 1) std::cout << "Ignoring G4hit, time too large: " << g4hit_i << " time: " << (*G4HitT)[g4hit_i] << std::endl;
- continue;
- }
-
- if (g4hits.find((*G4HitCellIdentifier)[g4hit_i]) == g4hits.end())
- {
- //this G4 hit doesn't exist yet
- one_g4hit.identifier = (*G4HitIdentifier)[g4hit_i];
- one_g4hit.cell_identifier = (*G4HitCellIdentifier)[g4hit_i];
- one_g4hit.sampling = (*G4HitSampling)[g4hit_i];
- one_g4hit.hit_time = (*G4HitT)[g4hit_i];
- //one_g4hit.hit_energy = (*G4HitE)[g4hit_i];
- //scale the hit energy with the sampling fraction
- if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20)
- { //tile
- //std::cout <<"Tile: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl;
- if ((*G4HitSamplingFraction)[g4hit_i])
- {
- one_g4hit.hit_energy = (*G4HitE)[g4hit_i] * (*G4HitSamplingFraction)[g4hit_i];
- //std::cout <<"TileE: "<<one_g4hit.hit_energy<<std::endl;
- }
- else one_g4hit.hit_energy = 0.;
- }
- else
- {
- //std::cout <<"LAr: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl;
- one_g4hit.hit_energy = (*G4HitE)[g4hit_i] / (*G4HitSamplingFraction)[g4hit_i];
- }
- //one_g4hit.hit_sampfrac = (*G4HitSamplingFraction)[g4hit_i];
- //new g4hit -> insert vector with 1 element
- g4hits.insert(std::pair<Long64_t, std::vector<FCS_g4hit> >(one_g4hit.cell_identifier, std::vector<FCS_g4hit>(1, one_g4hit)));
- }
- else
- {
- //G4 hit exists in this identifier -> push_back new to the vector //FCS_g4hit one_g4hit;
- one_g4hit.identifier = (*G4HitIdentifier)[g4hit_i];
- one_g4hit.cell_identifier = (*G4HitCellIdentifier)[g4hit_i];
- one_g4hit.sampling = (*G4HitSampling)[g4hit_i];
- one_g4hit.hit_time = (*G4HitT)[g4hit_i];
- if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20)
- { //tile
- //std::cout <<"Tile2: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl;
- if ((*G4HitSamplingFraction)[g4hit_i])
- {
- one_g4hit.hit_energy = (*G4HitE)[g4hit_i] * (*G4HitSamplingFraction)[g4hit_i];
- //std::cout <<"TileE2: "<<one_g4hit.hit_energy<<std::endl;
- }
- else one_g4hit.hit_energy = 0.;
- }
- else
- {
- //std::cout <<"LAr2: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl;
- one_g4hit.hit_energy = (*G4HitE)[g4hit_i] / (*G4HitSamplingFraction)[g4hit_i];
- }
- //one_g4hit.hit_sampfrac = (*G4HitSamplingFraction)[g4hit_i];
- g4hits[(*G4HitCellIdentifier)[g4hit_i]].push_back(one_g4hit);
- }
- }
- }
-
- if (m_Debug > 1) std::cout << " Done" << std::endl;
-
- if (m_Debug > 1) std::cout << "Reading detailed FCS hits " << HitIdentifier->size() << std::endl;
-
- for (unsigned int hit_i = 0; hit_i < HitIdentifier->size(); hit_i++)
- {
- if ((*HitSampling)[hit_i] >= 0 && (*HitSampling)[hit_i] <= 25 && (*HitT)[hit_i] > m_TimingCut)
- {
- if (m_Debug > 1)
- std::cout << "Ignoring FCS hit, time too large: " << hit_i << " time: " << (*HitT)[hit_i] << std::endl;
- continue;
- }
- if (hits.find((*HitCellIdentifier)[hit_i]) == hits.end())
- {
- //Detailed hit doesn't exist yet
- one_hit.identifier = (*HitIdentifier)[hit_i];
- one_hit.cell_identifier = (*HitCellIdentifier)[hit_i];
- one_hit.sampling = (*HitSampling)[hit_i];
-
- if (one_hit.sampling >= 12 && one_hit.sampling <= 20)
- { //tile
- if ((*HitSamplingFraction)[hit_i])
- {
- one_hit.hit_energy = (*HitE)[hit_i] * (*HitSamplingFraction)[hit_i];
- }
- else one_hit.hit_energy = 0.;
- }
- else
- {
- one_hit.hit_energy = (*HitE)[hit_i] / (*HitSamplingFraction)[hit_i];
- }
- //one_hit.hit_sampfrac = (*HitSamplingFraction)[hit_i];
- one_hit.hit_time = (*HitT)[hit_i];
- one_hit.hit_x = (*HitX)[hit_i];
- one_hit.hit_y = (*HitY)[hit_i];
- one_hit.hit_z = (*HitZ)[hit_i];
- hits.insert(std::pair<Long64_t, std::vector<FCS_hit> >(one_hit.cell_identifier, std::vector<FCS_hit>(1, one_hit)));
- }
- else
- {
- //Detailed hit exists in this identifier -> push_back new to the vector
- one_hit.identifier = (*HitIdentifier)[hit_i];
- one_hit.cell_identifier = (*HitCellIdentifier)[hit_i];
- one_hit.sampling = (*HitSampling)[hit_i];
- //one_hit.hit_energy = (*HitE)[hit_i];
- if (one_hit.sampling >= 12 && one_hit.sampling <= 20)
- { //tile
- if ((*HitSamplingFraction)[hit_i])
- {
- one_hit.hit_energy = (*HitE)[hit_i] * (*HitSamplingFraction)[hit_i];
- }
- else one_hit.hit_energy = 0.;
- }
- else
- {
- one_hit.hit_energy = (*HitE)[hit_i] / (*HitSamplingFraction)[hit_i];
- }
- //one_hit.hit_sampfrac = (*HitSamplingFraction)[hit_i];
- one_hit.hit_time = (*HitT)[hit_i];
- one_hit.hit_x = (*HitX)[hit_i];
- one_hit.hit_y = (*HitY)[hit_i];
- one_hit.hit_z = (*HitZ)[hit_i];
- hits[(*HitCellIdentifier)[hit_i]].push_back(one_hit);
- }
- }
-
- if (m_Debug > 1) std::cout << " Done" << std::endl;
-
- if (m_Debug > 1) std::cout << "ISF_HitAnalysis Check: Cells: " << cells.size() << " G4hits: " << g4hits.size() << " FCS detailed hits: " << hits.size() << std::endl;
-
- //Start matching:
- Int_t iindex = 0;
- for (std::map<Long64_t, FCS_cell>::iterator it = cells.begin(); it != cells.end(); )
- {
- iindex++;
- // std::cout <<iindex<<std::endl;
- one_matchedcell.clear(); //maybe not completely necessery, as we're not pushing_back into vectors
- //set the cell part
- one_matchedcell.cell = it->second;
- //now look for FCS detailed hits in this cell
- std::map<Long64_t, std::vector<FCS_hit> >::iterator it2 = hits.find(it->first);
- if (it2 != hits.end())
- {
- //std::cout <<"FCS hits found in this cell"<<std::endl;
- one_matchedcell.hit = it2->second;
- hits.erase(it2); //remove it
- }
- else
- {
- //no hit found for this cell
- one_matchedcell.hit.clear(); //important!
- }
- //now look for G4hits in this cell
- std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first);
- if (it3 != g4hits.end())
- {
- //std::cout <<"G4 hits found in this cell"<<std::endl;
- one_matchedcell.g4hit = it3->second;
- g4hits.erase(it3);
- }
- else
- {
- //no g4hit found for this cell
- one_matchedcell.g4hit.clear();//important!
- }
- //std::cout <<"Erase cell"<<std::endl;
- cells.erase(it++);
- //std::cout <<"Insert matched object"<<std::endl;
- //push_back matched cell for event jentry
- oneeventcells->push_back(one_matchedcell);
- //std::cout <<"Go next"<<std::endl;
-
- }
-
- //ok, cells should be empty, what about hits and g4hits?
- //There could be G4hits/FCS hits for which we don't have a cell ->create a dummy empty cell with 0 energy, take the cell identifier from the hit
- if (m_Debug > 1) std::cout << "ISF_HitAnalysis Check after cells: " << cells.size() << " " << g4hits.size() << " " << hits.size() << std::endl;
-
- for (std::map<Long64_t, std::vector<FCS_hit> >::iterator it = hits.begin(); it != hits.end();)
- {
- one_matchedcell.clear();
- one_matchedcell.cell.cell_identifier = it->first;
- //std::cout <<"This hit didn't exist in cell: "<<it->first<<std::endl;
- if (it->second.size())
- {
- one_matchedcell.cell.sampling = (it->second)[0].sampling;
- }
- else
- {
- one_matchedcell.cell.sampling = -1; //
- //ok, but you really shouldn't be here
- std::cout << "ERROR: You shouldn't really be here" << std::endl;
- }
- one_matchedcell.cell.energy = 0.;
- one_matchedcell.cell.center_x = 0.0;
- one_matchedcell.cell.center_y = 0.0;
- one_matchedcell.cell.center_z = 0.0;
- one_matchedcell.hit = it->second;
- std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first);
- if (it3 != g4hits.end())
- {
- one_matchedcell.g4hit = it3->second;
- g4hits.erase(it3);
- }
- else
- {
- //no g4hit found for this cell
- one_matchedcell.g4hit.clear(); //important!
- }
- hits.erase(it++);
- oneeventcells->push_back(one_matchedcell);
-
- }
-
- //ok, hits should be empty, what about g4hits?
- if (m_Debug > 1 )std::cout << "ISF_HitAnalysis Check after hits: " << cells.size() << " " << g4hits.size() << " " << hits.size() << std::endl;
- for (std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it = g4hits.begin(); it != g4hits.end();)
- {
- one_matchedcell.clear(); //maybe not so important
- one_matchedcell.cell.cell_identifier = it->first;
- if (it->second.size())
- {
- one_matchedcell.cell.sampling = (it->second)[0].sampling;
- }
- else
- {
- one_matchedcell.cell.sampling = -1; //
- //not really
- std::cout << "ERROR: You shouldn't really be here" << std::endl;
- }
- one_matchedcell.cell.energy = 0.;
- one_matchedcell.cell.center_x = 0.0;
- one_matchedcell.cell.center_y = 0.0;
- one_matchedcell.cell.center_z = 0.0;
- one_matchedcell.g4hit = it->second;
- one_matchedcell.hit.clear(); //important!!
- g4hits.erase(it++);
- oneeventcells->push_back(one_matchedcell);
- }
-
- if (m_Debug > 1) std::cout << "ISF_HitAnalysis Check after g4hits: " << cells.size() << " " << g4hits.size() << " " << hits.size() << std::endl;
- //Final size for this event
- if (m_Debug > 1) std::cout << "ISF_HitAnalysis Matched cells size: " << oneeventcells->size() << std::endl;
-
- //Can fill the output tree already here:
- total_cell_e = 0;
- total_hit_e = 0;
- total_g4hit_e = 0;
-
- for (int j = 0; j < MAX_LAYER - 1; j++)
- {
- layercells[j]->m_vector = oneeventcells->GetLayer(j);
- }
-
- //this is for invalid cells
- layercells[MAX_LAYER - 1]->m_vector = oneeventcells->GetLayer(-1);
-
- for (int i = 0; i < MAX_LAYER; i++)
- {
- cell_energy->at(i) = 0.0; //zero for each event!
- hit_energy->at(i) = 0.0;
- g4hit_energy->at(i) = 0.0;
-
- for (unsigned int cellindex = 0; cellindex < layercells[i]->size(); cellindex++)
- {
- if (i != MAX_LAYER - 1)
- {
- cell_energy->at(i) += layercells[i]->m_vector.at(cellindex).cell.energy;
- total_cell_e += layercells[i]->m_vector.at(cellindex).cell.energy;
- }
- else
- {
- //don't add the energy in the invalid layer to the total energy (if there is any (shouldn't)
- cell_energy->at(i) += layercells[i]->m_vector.at(cellindex).cell.energy; //this should be here anyway
- }
-
- //sum energy of all FCS detailed hits in this layer/cell
- for (unsigned int j = 0; j < layercells[i]->m_vector.at(cellindex).hit.size(); j++)
- {
- if (i != MAX_LAYER - 1)
- {
- total_hit_e += layercells[i]->m_vector.at(cellindex).hit[j].hit_energy;
- hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).hit[j].hit_energy;
- }
- else
- {
- //again, don't add invalid layer energy to the sum
- hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).hit[j].hit_energy;
- }
- }
-
- //sum energy of all G4 hits in this layer/cell
- for (unsigned int j = 0; j < layercells[i]->m_vector.at(cellindex).g4hit.size(); j++)
- {
- if (i != MAX_LAYER - 1)
- {
- total_g4hit_e += layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy;
- g4hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy;
- }
- else
- {
- //don't add invalied layer energy to the sum
- g4hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy;
- }
- }
- }
- }
-
- cell_energy->at(MAX_LAYER) = total_cell_e;
- hit_energy->at(MAX_LAYER) = total_hit_e;
- g4hit_energy->at(MAX_LAYER) = total_g4hit_e;
-
- //Uncomment this to get energy cross check
- if (jentry % m_PrintOutFrequency == 0) std::cout << "Energy check: event: " << jentry << " cell: " << total_cell_e << " g4hits: " << total_g4hit_e << " hits: " << total_hit_e << std::endl;
- //tree gets filled
- m_OutputTree->Fill();
-
- }//loop over entries
-
- m_OutputTree->Write();
- m_Output->Close();
-};
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.h
deleted file mode 100644
index 7af0d141ca9955f74a1e7fa65eed93320ee4ecf1..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.h
+++ /dev/null
@@ -1,521 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-//////////////////////////////////////////////////////////
-// This class has been automatically generated on
-// Fri May 2 18:25:08 2014 by ROOT version 5.34/05
-// from TTree CaloHitAna/CaloHitAna
-// found on file: ISF_HitAnalysispion_eta1_020514.root
-//////////////////////////////////////////////////////////
-
-#ifndef CaloHitAna_h
-#define CaloHitAna_h
-
-#include <TROOT.h>
-#include <TChain.h>
-#include <TFile.h>
-
-// Header file for the classes stored in the TTree if any.
-#include <vector>
-#include <FCS_Cell.h>
-
-#include <iostream>
-
-class CaloHitAna
-{
-public :
- TTree *fChain; //!pointer to the analyzed TTree or TChain
- Int_t fCurrent; //!current Tree number in a TChain
- TString fFilename;
-
- //output structure
- TFile *m_Output;
- TString m_OutputName;
- TTree *m_OutputTree;
- std::vector<Int_t> m_Settings;
- Float_t m_TimingCut;
- Int_t m_Debug;
- Int_t m_PrintOutFrequency;
- Int_t m_max_nentries;
- bool m_do_g4_hits = false;
-
- const static int MAX_LAYER = 25;
-
- // Declaration of leaf types
- std::vector<float> *HitX;
- std::vector<float> *HitY;
- std::vector<float> *HitZ;
- std::vector<float> *HitE;
- std::vector<float> *HitT;
- std::vector<Long64_t> *HitIdentifier;
- std::vector<Long64_t> *HitCellIdentifier;
- std::vector<bool> *HitIsLArBarrel;
- std::vector<bool> *HitIsLArEndCap;
- std::vector<bool> *HitIsHEC;
- std::vector<bool> *HitIsFCAL;
- std::vector<bool> *HitIsTile;
- std::vector<int> *HitSampling;
- std::vector<float> *HitSamplingFraction;
- std::vector<float> *TruthE;
- std::vector<float> *TruthPx;
- std::vector<float> *TruthPy;
- std::vector<float> *TruthPz;
- std::vector<int> *TruthPDG;
- std::vector<int> *TruthBarcode;
- std::vector<int> *TruthVtxBarcode;
- std::vector<Long64_t> *CellIdentifier;
- std::vector<float> *CellE;
- std::vector<int> *CellSampling;
- std::vector<float> *G4HitE;
- std::vector<float> *G4HitT;
- std::vector<Long64_t> *G4HitIdentifier;
- std::vector<Long64_t> *G4HitCellIdentifier;
- std::vector<float> *G4HitSamplingFraction;
- std::vector<int> *G4HitSampling;
-
- //input (from is hit analysis ntuple):
- std::vector<std::vector<float> >* newTTC_entrance_eta;
- std::vector<std::vector<float> >* newTTC_entrance_phi;
- std::vector<std::vector<float> >* newTTC_entrance_r;
- std::vector<std::vector<float> >* newTTC_entrance_z;
- std::vector<std::vector<float> >* newTTC_back_eta;
- std::vector<std::vector<float> >* newTTC_back_phi;
- std::vector<std::vector<float> >* newTTC_back_r;
- std::vector<std::vector<float> >* newTTC_back_z;
- std::vector<std::vector<float> >* newTTC_mid_eta;
- std::vector<std::vector<float> >* newTTC_mid_phi;
- std::vector<std::vector<float> >* newTTC_mid_r;
- std::vector<std::vector<float> >* newTTC_mid_z;
- std::vector<float>* newTTC_IDCaloBoundary_eta;
- std::vector<float>* newTTC_IDCaloBoundary_phi;
- std::vector<float>* newTTC_IDCaloBoundary_r;
- std::vector<float>* newTTC_IDCaloBoundary_z;
- std::vector<float>* newTTC_Angle3D;
- std::vector<float>* newTTC_AngleEta;
-
- //output:
- std::vector<std::vector<float> >* m_newTTC_entrance_eta = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_entrance_phi = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_entrance_r = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_entrance_z = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_back_eta = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_back_phi = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_back_r = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_back_z = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_mid_eta = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_mid_phi = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_mid_r = new std::vector<std::vector<float>>;
- std::vector<std::vector<float> >* m_newTTC_mid_z = new std::vector<std::vector<float>>;
- std::vector<float>* m_newTTC_IDCaloBoundary_eta = new std::vector<float>;
- std::vector<float>* m_newTTC_IDCaloBoundary_phi = new std::vector<float>;
- std::vector<float>* m_newTTC_IDCaloBoundary_r = new std::vector<float>;
- std::vector<float>* m_newTTC_IDCaloBoundary_z = new std::vector<float>;
- std::vector<float>* m_newTTC_Angle3D = new std::vector<float>;
- std::vector<float>* m_newTTC_AngleEta = new std::vector<float>;
-
- FCS_matchedcellvector* oneeventcells = new FCS_matchedcellvector; //these are all matched cells in a single event
- FCS_matchedcellvector* layercells[MAX_LAYER]; //these are all matched cells in a given layer in a given event
-
- std::vector<FCS_truth>* truthcollection = new std::vector<FCS_truth>;
-
- Float_t total_cell_e = 0;
- Float_t total_hit_e = 0;
- Float_t total_g4hit_e = 0;
-
- std::vector<Float_t>* cell_energy = new std::vector<Float_t>(MAX_LAYER + 1);
- std::vector<Float_t>* hit_energy = new std::vector<Float_t>(MAX_LAYER + 1);
- std::vector<Float_t>* g4hit_energy = new std::vector<Float_t>(MAX_LAYER + 1);
-
- std::vector<Float_t>* new_truthE = new std::vector<Float_t>;
- std::vector<Float_t>* new_truthPx = new std::vector<Float_t>;
- std::vector<Float_t>* new_truthPy = new std::vector<Float_t>;
- std::vector<Float_t>* new_truthPz = new std::vector<Float_t>;
- std::vector<int>* new_truthBarcode = new std::vector<int>;
- std::vector<int>* new_truthPDG = new std::vector<int>;
- std::vector<int>* new_truthVtxBarcode = new std::vector<int>;
-
- // List of branches
- TBranch *b_HitX; //!
- TBranch *b_HitY; //!
- TBranch *b_HitZ; //!
- TBranch *b_HitE; //!
- TBranch *b_HitT; //!
- TBranch *b_HitIdentifier; //!
- TBranch *b_HitCellIdentifier; //!
- TBranch *b_HitIsLArBarrel; //!
- TBranch *b_HitIsLArEndCap; //!
- TBranch *b_HitIsHEC; //!
- TBranch *b_HitIsFCAL; //!
- TBranch *b_HitIsTile; //!
- TBranch *b_HitSampling; //!
- TBranch *b_HitSamplingFraction; //!
- TBranch *b_TruthE; //!
- TBranch *b_TruthPx; //!
- TBranch *b_TruthPy; //!
- TBranch *b_TruthPz; //!
- TBranch *b_TruthPDG; //!
- TBranch *b_TruthBarcode; //!
- TBranch *b_TruthVtxBarcode; //!
- TBranch *b_CellIdentifier; //!
- TBranch *b_CellE; //!
- TBranch *b_CellSampling; //!
- TBranch *b_G4HitE; //!
- TBranch *b_G4HitT; //!
- TBranch *b_G4HitIdentifier; //!
- TBranch *b_G4HitCellIdentifier; //!
- TBranch *b_G4HitSamplingFraction; //!
- TBranch *b_G4HitSampling; //!
- /*
- TBranch *b_TTC_back_eta; //!
- TBranch *b_TTC_back_phi; //!
- TBranch *b_TTC_back_r; //!
- TBranch *b_TTC_back_z; //!
- TBranch *b_TTC_entrance_eta; //!
- TBranch *b_TTC_entrance_phi; //!
- TBranch *b_TTC_entrance_r; //!
- TBranch *b_TTC_entrance_z; //!
- TBranch *b_TTC_IDCaloBoundary_eta; //!
- TBranch *b_TTC_IDCaloBoundary_phi; //!
- TBranch *b_TTC_IDCaloBoundary_r; //!
- TBranch *b_TTC_IDCaloBoundary_z; //!
- TBranch *b_TTC_Angle3D; //!
- TBranch *b_TTC_AngleEta; //!
- */
-
- TBranch *b_newTTC_back_eta; //!
- TBranch *b_newTTC_back_phi; //!
- TBranch *b_newTTC_back_r; //!
- TBranch *b_newTTC_back_z; //!
- TBranch *b_newTTC_entrance_eta; //!
- TBranch *b_newTTC_entrance_phi; //!
- TBranch *b_newTTC_entrance_r; //!
- TBranch *b_newTTC_entrance_z; //!
- TBranch *b_newTTC_mid_eta; //!
- TBranch *b_newTTC_mid_phi; //!
- TBranch *b_newTTC_mid_r; //!
- TBranch *b_newTTC_mid_z; //!
- TBranch *b_newTTC_IDCaloBoundary_eta; //!
- TBranch *b_newTTC_IDCaloBoundary_phi; //!
- TBranch *b_newTTC_IDCaloBoundary_r; //!
- TBranch *b_newTTC_IDCaloBoundary_z; //!
- TBranch *b_newTTC_Angle3D; //!
- TBranch *b_newTTC_AngleEta; //!
-
-
- CaloHitAna(TString filename = "ISF_HitAnalysispion_eta1.root", TString outputname = "output1.root", std::vector<Int_t> settings = std::vector<Int_t>(), Float_t timingcut = 999999., Int_t debug = 0, TTree *tree = 0);
- virtual ~CaloHitAna();
- virtual Int_t Cut(Long64_t entry);
- virtual Int_t GetEntry(Long64_t entry);
- virtual Long64_t LoadTree(Long64_t entry);
- virtual void Init(TTree *tree);
- virtual void InitOutTree();
- virtual void Loop();
- virtual Bool_t Notify();
- virtual void Show(Long64_t entry = -1);
-
- void SetDebug(Int_t debug = 0) {m_Debug = debug;};
- void SetTimingCut(Int_t timingcut = 999999) {m_TimingCut = timingcut;};
- void SetDoAllCells(Int_t doit = 0) {if (m_Settings.size() < 1) m_Settings.resize(1, 0); m_Settings[0] = doit;};
- void SetDoLayers(Int_t doit = 0) {if (m_Settings.size() < 2) m_Settings.resize(2, 0); m_Settings[1] = doit;};
- void SetDoLayerSums(Int_t doit = 0) {if (m_Settings.size() < 3) m_Settings.resize(3, 0); m_Settings[2] = doit;};
-};
-
-#endif
-
-#ifdef CaloHitAna_cxx
-CaloHitAna::CaloHitAna(TString filename, TString outputname, std::vector<Int_t> settings, Float_t timingcut, Int_t debug, TTree *tree) : fChain(0)
-{
- fFilename = filename;
- // if parameter tree is not specified (or zero), connect the file
- // used to generate this class and read the Tree.
- if (tree == 0) {
- TFile *f = TFile::Open(filename, "READ");
- TString dirname = filename;
- dirname += ":/ISF_HitAnalysis";
- TDirectory * dir = (TDirectory*)f->Get(dirname);
- dir->GetObject("CaloHitAna", tree);
- }
- Init(tree);
- m_Settings = settings;
- m_Debug = debug;
- m_TimingCut = timingcut;
- m_OutputName = outputname;
- m_PrintOutFrequency = 100;
- m_max_nentries = -1;
- m_Output = new TFile(m_OutputName, "RECREATE");
- m_OutputTree = new TTree("FCS_ParametrizationInput", "Output_Matched_cell_Tree");
- InitOutTree();
-}
-
-CaloHitAna::~CaloHitAna()
-{
- // if (m_OutputTree) delete m_OutputTree;
- // if (m_Output) delete m_Output;
- if (!fChain) return;
- delete fChain->GetCurrentFile();
-}
-
-Int_t CaloHitAna::GetEntry(Long64_t entry)
-{
- // Read contents of entry.
- if (!fChain) return 0;
- return fChain->GetEntry(entry);
-}
-Long64_t CaloHitAna::LoadTree(Long64_t entry)
-{
- // Set the environment to read one entry
- if (!fChain) return -5;
- Long64_t centry = fChain->LoadTree(entry);
- if (centry < 0) return centry;
- if (fChain->GetTreeNumber() != fCurrent) {
- fCurrent = fChain->GetTreeNumber();
- Notify();
- }
- return centry;
-}
-
-void CaloHitAna::InitOutTree()
-{
- m_OutputTree->Branch("TruthE", &new_truthE);
- m_OutputTree->Branch("TruthPx", &new_truthPx);
- m_OutputTree->Branch("TruthPy", &new_truthPy);
- m_OutputTree->Branch("TruthPz", &new_truthPz);
- m_OutputTree->Branch("TruthPDG", &new_truthPDG);
- m_OutputTree->Branch("TruthBarcode", &new_truthBarcode);
- m_OutputTree->Branch("TruthVtxBarcode", &new_truthVtxBarcode); //this is duplicate of what is in the truth collection, will be good to remove/hide at some point
-
- m_OutputTree->Branch("newTTC_back_eta", &m_newTTC_back_eta);
- m_OutputTree->Branch("newTTC_back_phi", &m_newTTC_back_phi);
- m_OutputTree->Branch("newTTC_back_r", &m_newTTC_back_r);
- m_OutputTree->Branch("newTTC_back_z", &m_newTTC_back_z);
- m_OutputTree->Branch("newTTC_entrance_eta", &m_newTTC_entrance_eta);
- m_OutputTree->Branch("newTTC_entrance_phi", &m_newTTC_entrance_phi);
- m_OutputTree->Branch("newTTC_entrance_r", &m_newTTC_entrance_r);
- m_OutputTree->Branch("newTTC_entrance_z", &m_newTTC_entrance_z);
- m_OutputTree->Branch("newTTC_mid_eta", &m_newTTC_mid_eta);
- m_OutputTree->Branch("newTTC_mid_phi", &m_newTTC_mid_phi);
- m_OutputTree->Branch("newTTC_mid_r", &m_newTTC_mid_r);
- m_OutputTree->Branch("newTTC_mid_z", &m_newTTC_mid_z);
- m_OutputTree->Branch("newTTC_IDCaloBoundary_eta", &m_newTTC_IDCaloBoundary_eta);
- m_OutputTree->Branch("newTTC_IDCaloBoundary_phi", &m_newTTC_IDCaloBoundary_phi);
- m_OutputTree->Branch("newTTC_IDCaloBoundary_r", &m_newTTC_IDCaloBoundary_r);
- m_OutputTree->Branch("newTTC_IDCaloBoundary_z", &m_newTTC_IDCaloBoundary_z);
- m_OutputTree->Branch("newTTC_Angle3D", &m_newTTC_Angle3D);
- m_OutputTree->Branch("newTTC_AngleEta", &m_newTTC_AngleEta);
-
- //create branches in the output tree according to the settings vector
- if (! m_Settings.size() || m_Settings[0] == 1)
- {
- //Write all FCS_matchedcells
- m_OutputTree->Branch("AllCells", &oneeventcells);
- }
- if (m_Settings.size() >= 2 && m_Settings[1] == 1)
- {
- //write cells per layer
- for (Int_t i = 0; i < MAX_LAYER; i++)
- {
- TString branchname = "Sampling_";
- branchname += i;
- // std::cout<< "fail?" << std::endl;
- layercells[i] = new FCS_matchedcellvector;
- m_OutputTree->Branch(branchname, &layercells[i]);
- }
- }
- if (m_Settings.size() >= 3 && m_Settings[2] == 1)
- {
- //write also energies per layer:
- m_OutputTree->Branch("cell_energy", &cell_energy);
- m_OutputTree->Branch("hit_energy", &hit_energy);
- m_OutputTree->Branch("g4hit_energy", &g4hit_energy);
-
- //This is a duplicate of cell_energy[25]
- m_OutputTree->Branch("total_cell_energy", &total_cell_e);
- m_OutputTree->Branch("total_hit_energy", &total_hit_e);
- m_OutputTree->Branch("total_g4hit_energy", &total_g4hit_e);
- }
- // Enable/Disable recording of g4hits
- if (m_Settings.size() >= 4 && m_Settings[3] == 1) {
- std::cout << "do g4hits!" << std::endl;
- m_do_g4_hits = true;
- }
-}
-
-void CaloHitAna::Init(TTree *tree)
-{
- // The Init() function is called when the selector needs to initialize
- // a new tree or chain. Typically here the branch addresses and branch
- // pointers of the tree will be set.
- // It is normally not necessary to make changes to the generated
- // code, but the routine can be extended by the user if needed.
- // Init() will be called many times when running on PROOF
- // (once per file to be processed).
-
- // Set object pointer
- HitX = 0;
- HitY = 0;
- HitZ = 0;
- HitE = 0;
- HitT = 0;
- HitIdentifier = 0;
- HitCellIdentifier = 0;
- HitIsLArBarrel = 0;
- HitIsLArEndCap = 0;
- HitIsHEC = 0;
- HitIsFCAL = 0;
- HitIsTile = 0;
- HitSampling = 0;
- HitSamplingFraction = 0;
- TruthE = 0;
- TruthPx = 0;
- TruthPy = 0;
- TruthPz = 0;
- TruthPDG = 0;
- TruthBarcode = 0;
- TruthVtxBarcode = 0;
- CellIdentifier = 0;
- CellE = 0;
- CellSampling = 0;
- G4HitE = 0;
- G4HitT = 0;
- G4HitIdentifier = 0;
- G4HitCellIdentifier = 0;
- G4HitSamplingFraction = 0;
- G4HitSampling = 0;
- /*
- TTC_back_eta = 0;
- TTC_back_phi = 0;
- TTC_back_r = 0;
- TTC_back_z = 0;
- TTC_entrance_eta = 0;
- TTC_entrance_phi = 0;
- TTC_entrance_r = 0;
- TTC_entrance_z = 0;
- TTC_IDCaloBoundary_eta = 0;
- TTC_IDCaloBoundary_phi = 0;
- TTC_IDCaloBoundary_r = 0;
- TTC_IDCaloBoundary_z = 0;
- TTC_Angle3D = 0;
- TTC_AngleEta = 0;
- */
- newTTC_back_eta = 0;
- newTTC_back_phi = 0;
- newTTC_back_r = 0;
- newTTC_back_z = 0;
- newTTC_entrance_eta = 0;
- newTTC_entrance_phi = 0;
- newTTC_entrance_r = 0;
- newTTC_entrance_z = 0;
- newTTC_mid_eta = 0;
- newTTC_mid_phi = 0;
- newTTC_mid_r = 0;
- newTTC_mid_z = 0;
- newTTC_IDCaloBoundary_eta = 0;
- newTTC_IDCaloBoundary_phi = 0;
- newTTC_IDCaloBoundary_r = 0;
- newTTC_IDCaloBoundary_z = 0;
- newTTC_Angle3D = 0;
- newTTC_AngleEta = 0;
-
-
- // Set branch addresses and branch pointers
- if (!tree) return;
- fChain = tree;
- fCurrent = -1;
- fChain->SetMakeClass(1);
-
- fChain->SetBranchAddress("HitX", &HitX, &b_HitX);
- fChain->SetBranchAddress("HitY", &HitY, &b_HitY);
- fChain->SetBranchAddress("HitZ", &HitZ, &b_HitZ);
- fChain->SetBranchAddress("HitE", &HitE, &b_HitE);
- fChain->SetBranchAddress("HitT", &HitT, &b_HitT);
- fChain->SetBranchAddress("HitIdentifier", &HitIdentifier, &b_HitIdentifier);
- fChain->SetBranchAddress("HitCellIdentifier", &HitCellIdentifier, &b_HitCellIdentifier);
- fChain->SetBranchAddress("HitIsLArBarrel", &HitIsLArBarrel, &b_HitIsLArBarrel);
- fChain->SetBranchAddress("HitIsLArEndCap", &HitIsLArEndCap, &b_HitIsLArEndCap);
- fChain->SetBranchAddress("HitIsHEC", &HitIsHEC, &b_HitIsHEC);
- fChain->SetBranchAddress("HitIsFCAL", &HitIsFCAL, &b_HitIsFCAL);
- fChain->SetBranchAddress("HitIsTile", &HitIsTile, &b_HitIsTile);
- fChain->SetBranchAddress("HitSampling", &HitSampling, &b_HitSampling);
- fChain->SetBranchAddress("HitSamplingFraction", &HitSamplingFraction, &b_HitSamplingFraction);
- fChain->SetBranchAddress("TruthE", &TruthE, &b_TruthE);
- fChain->SetBranchAddress("TruthPx", &TruthPx, &b_TruthPx);
- fChain->SetBranchAddress("TruthPy", &TruthPy, &b_TruthPy);
- fChain->SetBranchAddress("TruthPz", &TruthPz, &b_TruthPz);
- fChain->SetBranchAddress("TruthPDG", &TruthPDG, &b_TruthPDG);
- fChain->SetBranchAddress("TruthBarcode", &TruthBarcode, &b_TruthBarcode);
- fChain->SetBranchAddress("TruthVtxBarcode", &TruthVtxBarcode, &b_TruthVtxBarcode);
- fChain->SetBranchAddress("CellIdentifier", &CellIdentifier, &b_CellIdentifier);
- fChain->SetBranchAddress("CellE", &CellE, &b_CellE);
- fChain->SetBranchAddress("CellSampling", &CellSampling, &b_CellSampling);
- fChain->SetBranchAddress("G4HitE", &G4HitE, &b_G4HitE);
- fChain->SetBranchAddress("G4HitT", &G4HitT, &b_G4HitT);
- fChain->SetBranchAddress("G4HitIdentifier", &G4HitIdentifier, &b_G4HitIdentifier);
- fChain->SetBranchAddress("G4HitCellIdentifier", &G4HitCellIdentifier, &b_G4HitCellIdentifier);
- fChain->SetBranchAddress("G4HitSamplingFraction", &G4HitSamplingFraction, &b_G4HitSamplingFraction);
- fChain->SetBranchAddress("G4HitSampling", &G4HitSampling, &b_G4HitSampling);
- /*
- fChain->SetBranchAddress("TTC_back_eta", &TTC_back_eta, &b_TTC_back_eta);
- fChain->SetBranchAddress("TTC_back_phi", &TTC_back_phi, &b_TTC_back_phi);
- fChain->SetBranchAddress("TTC_back_r", &TTC_back_r, &b_TTC_back_r);
- fChain->SetBranchAddress("TTC_back_z", &TTC_back_z, &b_TTC_back_z);
- fChain->SetBranchAddress("TTC_entrance_eta", &TTC_entrance_eta, &b_TTC_entrance_eta);
- fChain->SetBranchAddress("TTC_entrance_phi", &TTC_entrance_phi, &b_TTC_entrance_phi);
- fChain->SetBranchAddress("TTC_entrance_r", &TTC_entrance_r, &b_TTC_entrance_r);
- fChain->SetBranchAddress("TTC_entrance_z", &TTC_entrance_z, &b_TTC_entrance_z);
- fChain->SetBranchAddress("TTC_IDCaloBoundary_eta", &TTC_IDCaloBoundary_eta, &b_TTC_IDCaloBoundary_eta);
- fChain->SetBranchAddress("TTC_IDCaloBoundary_phi", &TTC_IDCaloBoundary_phi, &b_TTC_IDCaloBoundary_phi);
- fChain->SetBranchAddress("TTC_IDCaloBoundary_r", &TTC_IDCaloBoundary_r, &b_TTC_IDCaloBoundary_r);
- fChain->SetBranchAddress("TTC_IDCaloBoundary_z", &TTC_IDCaloBoundary_z, &b_TTC_IDCaloBoundary_z);
- fChain->SetBranchAddress("TTC_Angle3D", &TTC_Angle3D, &b_TTC_Angle3D);
- fChain->SetBranchAddress("TTC_AngleEta", &TTC_AngleEta, &b_TTC_AngleEta);
- */
- fChain->SetBranchAddress("newTTC_back_eta", &newTTC_back_eta, &b_newTTC_back_eta);
- fChain->SetBranchAddress("newTTC_back_phi", &newTTC_back_phi, &b_newTTC_back_phi);
- fChain->SetBranchAddress("newTTC_back_r", &newTTC_back_r, &b_newTTC_back_r);
- fChain->SetBranchAddress("newTTC_back_z", &newTTC_back_z, &b_newTTC_back_z);
- fChain->SetBranchAddress("newTTC_entrance_eta", &newTTC_entrance_eta, &b_newTTC_entrance_eta);
- fChain->SetBranchAddress("newTTC_entrance_phi", &newTTC_entrance_phi, &b_newTTC_entrance_phi);
- fChain->SetBranchAddress("newTTC_entrance_r", &newTTC_entrance_r, &b_newTTC_entrance_r);
- fChain->SetBranchAddress("newTTC_entrance_z", &newTTC_entrance_z, &b_newTTC_entrance_z);
- fChain->SetBranchAddress("newTTC_mid_eta", &newTTC_mid_eta, &b_newTTC_mid_eta);
- fChain->SetBranchAddress("newTTC_mid_phi", &newTTC_mid_phi, &b_newTTC_mid_phi);
- fChain->SetBranchAddress("newTTC_mid_r", &newTTC_mid_r, &b_newTTC_mid_r);
- fChain->SetBranchAddress("newTTC_mid_z", &newTTC_mid_z, &b_newTTC_mid_z);
- fChain->SetBranchAddress("newTTC_IDCaloBoundary_eta", &newTTC_IDCaloBoundary_eta, &b_newTTC_IDCaloBoundary_eta);
- fChain->SetBranchAddress("newTTC_IDCaloBoundary_phi", &newTTC_IDCaloBoundary_phi, &b_newTTC_IDCaloBoundary_phi);
- fChain->SetBranchAddress("newTTC_IDCaloBoundary_r", &newTTC_IDCaloBoundary_r, &b_newTTC_IDCaloBoundary_r);
- fChain->SetBranchAddress("newTTC_IDCaloBoundary_z", &newTTC_IDCaloBoundary_z, &b_newTTC_IDCaloBoundary_z);
- fChain->SetBranchAddress("newTTC_Angle3D", &newTTC_Angle3D, &b_newTTC_Angle3D);
- fChain->SetBranchAddress("newTTC_AngleEta", &newTTC_AngleEta, &b_newTTC_AngleEta);
-
- Notify();
-}
-
-Bool_t CaloHitAna::Notify()
-{
- // The Notify() function is called when a new file is opened. This
- // can be either for a new TTree in a TChain or when when a new TTree
- // is started when using PROOF. It is normally not necessary to make changes
- // to the generated code, but the routine can be extended by the
- // user if needed. The return value is currently not used.
-
- return kTRUE;
-}
-
-void CaloHitAna::Show(Long64_t entry)
-{
- // Print contents of entry.
- // If entry is not specified, print current entry
- if (!fChain) return;
- fChain->Show(entry);
-}
-Int_t CaloHitAna::Cut(Long64_t entry)
-{
- // This function may be called from Loop.
- // returns 1 if entry is accepted.
- // returns -1 otherwise.
- std::cout << entry << std::endl;
- return 1;
-}
-#endif // #ifdef CaloHitAna_cxx
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna2.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna2.C
deleted file mode 100644
index 5e4c61209328e12f2d1d8d14f4c0391d67d4c9a3..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna2.C
+++ /dev/null
@@ -1,171 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#include "FCS_Cell.h"
-#include "TFile.h"
-#include "TTree.h"
-#include "TString.h"
-#include <string>
-#include <sstream>
-#include <iostream>
-#include "TSystem.h"
-#include "TString.h"
-#include "TFile.h"
-#include <stdlib.h>
-#include "TLorentzVector.h"
-#include "TH2D.h"
-
-void CaloHitAna2(TString filename="output1.root", Int_t layer=0)
-{
- gSystem->Load("FCS_Cell_h.so");
- TFile *fin = TFile::Open(filename,"READ");
- if (!fin) abort();
- TTree *t = (TTree*) fin->Get("FCS_ParametrizationInput");
- if (!t) abort();
- TString foutname = filename;
- TString outname=TString::Format("_output_%i.root",layer);
- foutname.ReplaceAll(".root",outname);
-
- TFile *fout = new TFile(foutname, "RECREATE");
- TH2D *h1 = new TH2D("h1","total_cell_energy",100,-5,5,600,0,60000);
- TH2D *h2 = new TH2D("h2","total_hit_energy",100,-5,5,600,0,60000);
- TH2D *h3 = new TH2D("h3","total_g4hit_energy",100,-5,5,600,0,60000);
- TH2D *hh1 = new TH2D("hh1","nhits",100,-5,5,1000,0,10000);
- TH2D *hh2 = new TH2D("hh2","ng4hits",100,-5,5,1000,0,10000);
-
- TH2D *hdiff1 = new TH2D("hdiff1","cell-hit_energy",100,-5,5,200,-10000,10000);
- TH2D *hdiff2 = new TH2D("hdiff2","cell-g4hit_energy",100,-5,5,200,-10000,10000);
-
- TString hdiffL1name; hdiffL1name.Form("hdiff1_%d",layer);
- TString hdiffL2name; hdiffL2name.Form("hdiff2_%d",layer);
-
- TH2D *hdiffL1 = new TH2D(hdiffL1name,hdiffL1name,100,-5,5,200,-10000,10000);
- TH2D *hdiffL2 = new TH2D(hdiffL2name,hdiffL2name,100,-5,5,200,-10000,10000);
-
- h1->Sumw2();
- h2->Sumw2();
- h3->Sumw2();
- hh1->Sumw2();
- hh2->Sumw2();
- hdiff1->Sumw2();
- hdiff2->Sumw2();
- hdiffL1->Sumw2();
- hdiffL2->Sumw2();
-
- FCS_matchedcellvector *vec=0;
- TString name="Sampling_";
- name+=layer;
- if (layer == -1)
- {
- name="AllCells";
- }
- t->SetBranchAddress(name,&vec);
- Float_t total_cell_energy;
- Float_t total_hit_energy;
- Float_t total_g4hit_energy;
- t->SetBranchAddress("total_cell_energy", &total_cell_energy);
- t->SetBranchAddress("total_hit_energy", &total_hit_energy);
- t->SetBranchAddress("total_g4hit_energy",&total_g4hit_energy);
- std::vector<float> *truth_px=0;
- std::vector<float> *truth_py = 0;
- std::vector<float> *truth_pz = 0;
- std::vector<float> *truth_e = 0;
-
-
- t->SetBranchAddress("TruthE",&truth_e);
- t->SetBranchAddress("TruthPx",&truth_px);
- t->SetBranchAddress("TruthPy",&truth_py);
- t->SetBranchAddress("TruthPz",&truth_pz);
-
- std::cout <<"Entries: "<<t->GetEntries()<<std::endl;
- for (Long64_t i=0; i<t->GetEntries(); i++)
- {
- t->GetEntry(i);
- if (i%100 ==0) std::cout <<"Event: "<<i<<" layer "<<layer<< " has: "<<vec->size()<<" cells"<<std::endl;
- UInt_t totalhits=0;
- UInt_t totalg4hits=0;
- Float_t totalhitE = 0.; //these are per layer!
- Float_t totalcellE = 0.;
- Float_t totalg4hitE = 0.0;
- for (UInt_t j=0; j<(*vec).size(); j++)
- {
- totalhits+=(*vec)[j].hit.size();
- totalg4hits+=(*vec)[j].g4hit.size();
- for (UInt_t k=0; k<(*vec)[j].hit.size(); k++)
- {
- totalhitE+=(*vec)[j].hit[k].hit_energy;
- }
- for (UInt_t k=0; k<(*vec)[j].g4hit.size(); k++)
- {
- totalg4hitE+=(*vec)[j].g4hit[k].hit_energy;
- }
- totalcellE+=(*vec)[j].cell.energy;
- }
- TLorentzVector pion;
- /*
- if (truth_e->size()>1)
- {
- std::cout <<"Truth: "<<truth_e->size()<<" "<<(*truth_px)[0]<<std::endl;
- }
- else std::cout <<truth_e->size()<<std::endl;
- */
- if (truth_e->size()>=1)
- {
- pion.SetPxPyPzE((*truth_px)[0], (*truth_py)[0], (*truth_pz)[0], (*truth_e)[0]);
- }
- else continue;
- /*
- std::cout <<std::endl;
- pion.Print();
- std::cout <<std::endl;
- */
- //std::cout <<"Pion rap: "<<pion.Rapidity()<<std::endl;
- h1->Fill(pion.Rapidity(), total_cell_energy);
- h2->Fill(pion.Rapidity(), total_hit_energy);
- h3->Fill(pion.Rapidity(), total_g4hit_energy);
- hh1->Fill(pion.Rapidity(), totalhits);
- hh2->Fill(pion.Rapidity(), totalg4hits);
- hdiff1->Fill(pion.Rapidity(), total_cell_energy-total_hit_energy);
- hdiff2->Fill(pion.Rapidity(), total_cell_energy-total_g4hit_energy);
- hdiffL1->Fill(pion.Rapidity(), totalcellE-totalhitE);
- hdiffL2->Fill(pion.Rapidity(), totalcellE-totalg4hitE);
- if (i%100 ==0) std::cout <<" and "<<totalhits<<" hits and "<<totalg4hits<<" g4hits"<<std::endl;
- }
- fout->cd();
- h1->Write();
- h2->Write();
- h3->Write();
- hh1->Write();
- hh2->Write();
- hdiff1->Write();
- hdiff2->Write();
- hdiffL1->Write();
- hdiffL2->Write();
- fout->Close();
-
-};
-
-int main(int argc, char *argv[])
-{
- TString filename="output1.root";
- Int_t layer =2;
- if (argc <=1)
- {
- std::cout <<"Usage: ./CaloHitAna2 inputfile sampling "<<std::endl;
- exit(1);
- }
-
- if (argc>=2)
- {
- filename=argv[1];
- }
- if (argc>=3)
- {
- TString param = argv[2];
- layer = param.Atoi();
- }
- std::cout <<"Parameters: inputfile: "<<filename<<" sampling: "<<layer<<std::endl;
- CaloHitAna2(filename, layer);
- return 0;
-}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna3.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna3.C
deleted file mode 100755
index 15b8dee9f55480efacffb5ac1272890e9e62cdb4..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna3.C
+++ /dev/null
@@ -1,202 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#include "FCS_Cell.h"
-#include "TFile.h"
-#include "TTree.h"
-#include "TString.h"
-#include <string>
-#include <sstream>
-#include <iostream>
-#include "TSystem.h"
-#include "TString.h"
-#include "TFile.h"
-#include <stdlib.h>
-#include "TLorentzVector.h"
-#include "TH2D.h"
-
-//another example how to read matched ntuples
-
-void CaloHitAna3(TString filename="output1.root", Int_t llayer=0, Long64_t max_entries=-1, float timing_cut=99999)
-{
- gSystem->Load("./FCS_Cell_h.so");
- TFile *fin = TFile::Open(filename,"READ");
- if (!fin) abort();
- TTree *t = (TTree*) fin->Get("FCS_ParametrizationInput");
- if (!t) abort();
- TString foutname = filename;
- TString outname="_output_";
- outname+=llayer;
- outname+=".root";
- foutname.ReplaceAll(".root",outname);
-
- TFile *fout = new TFile(foutname, "RECREATE");
-
-
- FCS_matchedcellvector *vec;
- TString name="Sampling_";
- name+=llayer;
- if (llayer == -1)
- {
- name="AllCells";
- }
- std::cout <<"Reading branch: "<<name<<" for layer: "<<llayer<<std::endl;
- t->SetBranchAddress(name,&vec);
-
- std::vector<float> *truth_px=0;
- std::vector<float> *truth_py = 0;
- std::vector<float> *truth_pz = 0;
- std::vector<float> *truth_e = 0;
-
-
- t->SetBranchAddress("TruthE",&truth_e);
- t->SetBranchAddress("TruthPx",&truth_px);
- t->SetBranchAddress("TruthPy",&truth_py);
- t->SetBranchAddress("TruthPz",&truth_pz);
-
- std::cout <<"Entries: "<<t->GetEntries()<<std::endl;
- if (max_entries == -1 )
- {
- max_entries = t->GetEntries();
- }
-
- for (Long64_t i=0; i<max_entries; i++)
- {
- //std::cout <<llayer<<std::endl;
- if (i%100 ==0) std::cout<<std::endl;
- t->GetEntry(i);
- //if (i%100 ==0) std::cout <<"Event: "<<i<<" layer "<<layer<< " "<<name<<" has: "<<vec->size()<<" cells"<<std::endl;
- //sort cells and hits
- vec->sort();
- //trim timing
- if (timing_cut<999999.)
- {
- if (i%100==0) std::cout <<"Trimming: orig size: "<<vec->size()<<" cells"<<std::endl;
- vec->time_trim(timing_cut);
- if (i%100==0) std::cout <<"Remaining non empty cells: "<<vec->size()<<std::endl;
- }
- else
- {
- if (i==0) std::cout<<"No time trimming"<<std::endl;
- }
-
- if (i%100 ==0)
- {
- std::cout <<"======================================"<<std::endl;
- std::cout <<"= EVENT "<<i<<" ="<<std::endl;
- std::cout <<"======================================"<<std::endl;
- std::cout <<"Layer: "<<llayer<<" : "<<vec->size()<< " cells"<<std::endl;
-
- float totalhitelayer = 0.;
- float totalcellelayer = 0.;
- for (UInt_t j=0; j<vec->size(); j++)
- {
- FCS_matchedcell cell = (*vec)[j];
- FCS_cell ccell = cell.cell;
- float ccelle = ccell.energy;
- UInt_t hitsize = cell.hit.size();
- std::cout <<" Cell "<<j<<" E: "<<ccell.energy<<" id: "<<ccell.cell_identifier<<" hits: "<<hitsize<<std::endl;
- float tothite=0.;
- totalcellelayer+=ccell.energy;
- for (UInt_t k=0; k<hitsize; k++)
- {
- FCS_hit hit = cell.hit[k];
- float hite = hit.hit_energy;
- float hitt = hit.hit_time;
- // float hitt = (*vec)[j].hit[k].hit_time;
- tothite+=hite;
- totalhitelayer+=hite;
- if (i%100 ==0) std::cout <<" Hit "<<k<<" HitE: "<<hite<<" Time: "<<hitt<<std::endl;
- if (i%100 ==0 && hitt>1000.) std::cout <<"HERE"<<std::endl;
- }
-
- float frac = 0.;
- //long emptycell =0;
- if (fabs(ccelle)>0.)
- {
- frac = tothite/ccelle;
- }
- // else emptycell++;
- std::cout <<" Total hitE in cell: "<<tothite<<" / cellE: "<<ccelle<<" ( "<<frac<<" ) "<<std::endl;//<<" empty cells: "<<emptycell<<std::endl;
-
- }
- // std::cout <<"----"<<std::endl;
- //float frac2 = 0.;
- //if (fabs(totalcellelayer)>0.) frac2= totalhitelayer/totalcellelayer;
- std::cout <<"Total check: layer: SUMhit: "<<totalhitelayer<<" / SUMcell: "<<totalcellelayer<<std::endl;//" ( "<<frac2<<" ) "<<std::endl;
- std::cout <<"======================================"<<std::endl;
- std::cout <<std::endl;
- }
-
- /*
- UInt_t totalhits=0;
- UInt_t totalg4hits=0;
- Float_t totalhitE = 0.; //these are per layer!
- Float_t totalcellE = 0.;
- Float_t totalg4hitE = 0.0;
- for (UInt_t j=0; j<(*vec).size(); j++)
- {
- if (i%100 ==0)
- {
- std::cout <<"Hit size: "<<(*vec)[j].hit.size()<<" g4hit size: "<<(*vec)[j].g4hit.size()<<std::endl;
- std::cout <<"Reading in cell: "<<j<<std::endl;
- }
- totalhits+=(*vec)[j].hit.size();
- totalg4hits+=(*vec)[j].g4hit.size();
- for (UInt_t k=0; k<(*vec)[j].hit.size(); k++)
- {
- if (i%100 ==0) std::cout <<"Reading hit: "<<k<<" e: "<<(*vec)[j].hit[k].hit_energy<<" tot: "<<totalhitE<<std::endl;
- totalhitE+=(*vec)[j].hit[k].hit_energy;
- }
- for (UInt_t k=0; k<(*vec)[j].g4hit.size(); k++)
- {
- if (i%100 ==0) std::cout <<"Reading hit: "<<k<<" e: "<<(*vec)[j].g4hit[k].hit_energy<<" tot: "<<totalg4hitE<<std::endl;
- totalg4hitE+=(*vec)[j].g4hit[k].hit_energy;
- }
- totalcellE+=(*vec)[j].cell.energy;
- }
- if (i%100 ==0) std::cout <<"Test? "<<totalhits<<" hits and "<<totalg4hits<<" g4hits with C: "<<total_cell_energy<<" / "<<totalcellE<<" H: "<<total_hit_energy<<" / "<<totalhitE<<" G4: "<<total_g4hit_energy<<" / "<<totalg4hitE<<std::endl;
- */
- TLorentzVector pion;
- if (truth_e->size()>=1)
- {
- pion.SetPxPyPzE((*truth_px)[0], (*truth_py)[0], (*truth_pz)[0], (*truth_e)[0]);
- }
- else continue;
- }
- fout->cd();
- fout->Close();
-
-};
-
-int main(int argc, char *argv[])
-{
- TString filename="output1.root";
- Int_t layer =2;
- float timing_cut=999999;
- Long64_t max_entries = -1;
- std::cout <<"Usage: ./CaloHitAna3 inputfile sampling max_entries timing_cut"<<std::endl;
- if (argc>=2)
- {
- filename=argv[1];
- }
- if (argc>=3)
- {
- TString param = argv[2];
- layer = param.Atoi();
- }
- if (argc>=4)
- {
- TString param = argv[3];
- max_entries = param.Atoi();
- }
- if (argc>=5)
- {
- TString param = argv[4];
- timing_cut = param.Atof();
- }
- std::cout <<"Parameters: inputfile: "<<filename<<" sampling: "<<layer<<" max entries: "<<max_entries<<" timing cut: "<<timing_cut<<std::endl;
- CaloHitAna3(filename, layer,max_entries,timing_cut);
- return 0;
-}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/Linkdef.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/Linkdef.h
deleted file mode 100644
index 73dcb053904cc3e6d3e7d77c3119c6a95d55f1df..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/Linkdef.h
+++ /dev/null
@@ -1,18 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifdef __MAKECINT__
-#pragma link C++ struct FCS_cell+;
-#pragma link C++ struct FCS_hit+;
-#pragma link C++ struct FCS_g4hit+;
-#pragma link C++ struct std::vector<FCS_hit>+;
-#pragma link C++ struct std::vector<FCS_g4hit>+;
-#pragma link C++ struct FCS_matchedcell+;
-#pragma link C++ struct std::vector<FCS_matchedcell>+;
-#pragma link C++ struct FCS_matchedcellvector+;
-#pragma link C++ class std::vector<Float_t>+;
-#pragma link C++ struct FCS_truth+;
-#pragma link C++ struct std::vector<FCS_truth>+;
-#pragma link C++ class std::vector<std::vector<float> >+;
-#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/Makefile b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/Makefile
deleted file mode 100755
index c409a954cd88a33b32ddfdae674b52d21ae93537..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/Makefile
+++ /dev/null
@@ -1,26 +0,0 @@
-if [ -f FCS_CellDict.cxx ];
-then
- rm FCS_CellDict.cxx
- rm FCS_CellDict.h
-fi
-rootcint FCS_CellDict.cxx -c FCS_Cell.h Linkdef.h
-
-#compile
-g++ -g -c run.C FCS_CellDict.cxx -I $(root-config --incdir) -I . -fPIC $CPPEXPFLAGS -std=c++11
-
-#link executable
-g++ -g -o CaloHitAna run.o FCS_CellDict.o -L $(root-config --libdir) -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lpthread -Wl,-rpath,$(root-config --libdir) -lm -ldl -fPIC $CPPEXPFLAGS -std=c++11
-
-#this creates a shared library for ROOT (which can be loaded by hand)
-g++ -shared -o FCS_Cell_h.so FCS_CellDict.o -L $(root-config --libdir) -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lpthread -Wl,-rpath,$(root-config --libdir) -lm -ldl -fPIC $CPPEXPFLAGS -std=c++11
-
-#example
-g++ -g -c CaloHitAna2.C -I $(root-config --incdir) -I . -fPIC $CPPEXPFLAGS -std=c++11
-
-g++ -g -o CaloHitAna2 CaloHitAna2.o FCS_CellDict.o -L $(root-config --libdir) -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lpthread -Wl,-rpath,$(root-config --libdir) -lm -ldl -lstdc++ -fPIC $CPPEXPFLAGS -std=c++11
-
-#another example
-g++ -g -c CaloHitAna3.C FCS_CellDict.cxx -Wall -I $(root-config --incdir) -I . -fPIC $CPPEXPFLAGS -Wall -std=c++11
-
-g++ -g -o CaloHitAna3 CaloHitAna3.o FCS_CellDict.o -L $(root-config --libdir) -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lpthread -Wl,-rpath,$(root-config --libdir) -lm -ldl -lstdc++ -fPIC $CPPEXPFLAGS -Wall -std=c++11
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/make_caloHitAna.sh b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/make_caloHitAna.sh
deleted file mode 100644
index 0f6fbdc60a0b5518fa7349547daa4745dc375a8a..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/make_caloHitAna.sh
+++ /dev/null
@@ -1,4 +0,0 @@
-g++ -g -c run.C FCS_CellDict.cxx -I $(root-config --incdir) -I . -fPIC $CPPEXPFLAGS -std=c++11
-g++ -g -o CaloHitAna run.o FCS_CellDict.o -L $(root-config --libdir) -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lpthread -Wl,-rpath,$(root-config --libdir) -lm -ldl -fPIC $CPPEXPFLAGS -std=c++11
-
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run.C
deleted file mode 100644
index 36e0d327799b7b294dd721368e0685ff70cf65bd..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run.C
+++ /dev/null
@@ -1,91 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#include "CaloHitAna.C"
-#include <iostream>
-#include <stdlib.h>
-void run(TString filename="ISF_HitAnalysispion_eta1.root",TString outputname="output1.root", std::vector<Int_t> settings =std::vector<Int_t>(), Float_t TimingCut=999999., Int_t debug = 0)
-{
- //CaloHitAna t;
- CaloHitAna t(filename,outputname, settings, TimingCut,debug,0);
- t.Loop();
- //No need to call finish now
-}
-
-
-int main(int argc, char *argv[])
-{
- TString filename="ISF_HitAnalysispion_eta1.root";
- TString outputname="output1.root";
- if (argc <=1)
- {
- std::cout <<"Usage: ./CaloHitAna input.root output.root save_all_cells save_cells_per_sampling save_energies_per_sampling hit_timing_cut"<<std::endl;
- exit(1);
- }
- if (argc >=2)
- {
- filename = argv[1];
- }
- std::cout <<"Processing "<<filename<<std::endl;
-
- if (argc >=3)
- {
- outputname = argv[2];
- }
- std::cout <<"Saving output to "<<outputname<<std::endl;
- std::vector<Int_t> settings(0);
- Float_t TimingCut = 999999.;
- Int_t Debug = 0;
- if (argc >=4)
- {
- for (Int_t i=3; i<argc; i++)
- {
- TString parameter = argv[i];
- if (! (parameter.Contains("tc") || parameter.Contains("debug")))
- {
- Int_t pari = parameter.Atoi();
- settings.push_back(pari);
- }
- else if (parameter.Contains("tc"))
- {
- parameter.ReplaceAll("tc=","");
- TimingCut = parameter.Atof();
- }
- else if (parameter.Contains("debug"))
- {
- parameter.ReplaceAll("debug=","");
- Debug = parameter.Atoi();
- std::cout <<Debug<<std::endl;
- }
- }
- }
- else
- {
- settings.push_back(1);
- }
- std::cout <<"Settings: ";
- for (UInt_t i=0; i<settings.size(); i++)
- {
- std::cout <<settings[i]<<" ";
- }
- std::cout <<std::endl;
- std::cout <<"Timing cut: "<<TimingCut<<std::endl;
- std::cout <<"Debug: "<<Debug<<std::endl;
- TFile *ftest = TFile::Open(filename,"READ");
- if (!ftest)
- {
- std::cout <<"File "<<filename<<"doesn't exist!"<<std::endl;
- exit(1);
- }
- else
- {
- ftest->Close();
- }
- delete ftest;
- std::cout <<std::endl;
- std::cout <<"Running"<<std::endl;
- run(filename, outputname, settings, TimingCut,Debug);
-
- return 0;
-}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/simplerun.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/simplerun.C
deleted file mode 100644
index 40fff4c2b4ec81459afda87ce1c517dd272415c8..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/simplerun.C
+++ /dev/null
@@ -1,20 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-void simplerun(TString filename="ISF_calo__211__E10000_10000__eta20_20_Evts0-1000_z03350.pool.root",TString outputname="output1.root")
-{
- gROOT->LoadMacro("FCS_Cell.h+");
- gROOT->LoadMacro("CaloHitAna.C+");
- //CaloHitAna t;
- CaloHitAna t(filename,outputname);
- t.SetDoLayers(1);
- t.SetDoLayerSums(1);
- t.m_PrintOutFrequency=1;
- t.m_max_nentries=5;
- t.m_Debug=2;
- t.Loop();
- //No need to call finish now
-}
-
-