From 4e05e76e080dc510451d7aa22fed9d8c0e386d4a Mon Sep 17 00:00:00 2001
From: Petr Jacka <petr.jacka@cern.ch>
Date: Fri, 1 Sep 2017 11:03:06 +0200
Subject: [PATCH] Cleaned mess in CaloGeometry classes
---
.../CMakeLists.txt | 3 +-
.../CaloDetDescrElement.h | 268 ----
.../CaloGeoDetDescrElement.h | 249 ----
.../CaloGeoGeometry.h | 200 ---
.../CaloGeometry.h | 97 +-
.../CaloGeometryFromFile.h | 22 -
.../FCAL_ChannelMap.h | 286 ----
.../ICaloGeometry.h | 1 -
.../ISF_HitAnalysis.h | 17 +
.../ISF_FastCaloSimParametrizationConfig.py | 1 -
...tCaloSimParametrization_DigiPostInclude.py | 2 +-
...oSimParametrization_DigiTilePostInclude.py | 18 +-
...aloSimParametrization_SimPreInclude_1mm.py | 5 +
.../src/CaloGeoGeometry.cxx | 1220 -----------------
.../src/CaloGeometry.cxx | 493 ++-----
.../src/CaloGeometryFromCaloDDM.cxx | 24 +-
.../src/CaloGeometryFromCaloDDM.h | 21 -
.../src/CaloGeometryFromFile.cxx | 373 -----
.../src/FCAL_ChannelMap.cxx | 488 -------
.../src/FastCaloSimParamAlg.cxx | 4 +-
.../src/ISF_HitAnalysis.cxx | 2 +-
.../tools/CaloGeometryFromFile.cxx | 144 +-
.../tools/init_geo.C | 8 +-
.../tools/run_geo.C | 84 +-
.../ISF_FastCaloSimServices/CMakeLists.txt | 8 +-
.../python/AdditionalConfig.py | 61 +-
.../python/ISF_FastCaloSimServicesConfig.py | 23 +-
.../python/ISF_FastCaloSimServicesConfigDb.py | 4 +-
.../src/FastCaloSimSvcV2.cxx | 1173 +++++++++++-----
.../src/FastCaloSimSvcV2.h | 206 ++-
30 files changed, 1302 insertions(+), 4203 deletions(-)
delete mode 100755 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloDetDescrElement.h
delete mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h
delete mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoGeometry.h
delete mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometryFromFile.h
delete mode 100755 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/FCAL_ChannelMap.h
create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude_1mm.py
delete mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeoGeometry.cxx
delete mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.h
delete mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromFile.cxx
delete mode 100755 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FCAL_ChannelMap.cxx
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
index 0ff7df8d978c..5640dd4b4695 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
@@ -16,6 +16,7 @@ atlas_depends_on_subdirs( PUBLIC
DetectorDescription/Identifier
GaudiKernel
LArCalorimeter/LArElecCalib
+ LArCalorimeter/LArGeoModel
Simulation/Barcode/BarcodeEvent
Simulation/ISF/ISF_Core/ISF_Interfaces
Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent
@@ -63,7 +64,7 @@ atlas_add_root_dictionary( ISF_FastCaloSimParametrizationLib
EXTERNAL_PACKAGES ROOT HepPDT XercesC CLHEP HepMC Geant4 )
atlas_add_library( ISF_FastCaloSimParametrizationLib
- Root/*.cxx src/CaloGeoGeometry.cxx src/CaloGeometryFromFile.cxx
+ Root/*.cxx src/CaloGeometryFromCaloDDM.cxx src/CaloGeometryLookup.cxx
${ISF_FastCaloSimParametrizationLibDictSource}
PUBLIC_HEADERS ISF_FastCaloSimParametrization
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${HEPPDT_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloDetDescrElement.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloDetDescrElement.h
deleted file mode 100755
index 4a4ecf6b2f15..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloDetDescrElement.h
+++ /dev/null
@@ -1,268 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef CALOGEODETECTORELEMENT_H
-#define CALOGEODETECTORELEMENT_H
-
-class CaloGeometry;
-
-#include "Identifier/Identifier.h"
-#include <iostream>
-#include <iomanip>
-#include <cmath>
-
-class CaloGeoDetDescrElement
-{
- friend class CaloGeometry;
- public:
- CaloGeoDetDescrElement() {
- m_identify = 0;
- m_hash_id = 0;
- m_calosample = 0;
- m_eta = 0;
- m_phi = 0;
- m_deta = 0;
- m_dphi = 0;
- m_r = 0;
- m_eta_raw = 0;
- m_phi_raw = 0;
- m_r_raw = 0;
- m_dr = 0;
- m_x = 0;
- m_y = 0;
- m_z = 0;
- m_x_raw = 0;
- m_y_raw = 0;
- m_z_raw = 0;
- m_dx = 0;
- m_dy = 0;
- m_dz = 0;
- };
-
- /** @brief virtual destructor
- */
- virtual ~CaloGeoDetDescrElement() {};
-
- /** @brief cell eta
- * @copydoc CaloGeoDetDescrElement::m_eta
- */
- float eta() const;
- /** @brief cell phi
- * @copydoc CaloGeoDetDescrElement::m_phi
- */
- float phi() const;
- /** @brief cell r
- * @copydoc CaloGeoDetDescrElement::m_r
- */
- float r() const;
- /** @brief cell eta_raw
- * @copydoc CaloGeoDetDescrElement::m_eta_raw
- */
- float eta_raw() const;
- /** @brief cell phi_raw
- * @copydoc CaloGeoDetDescrElement::m_phi_raw
- */
- float phi_raw() const;
- /** @brief cell r_raw
- * @copydoc CaloGeoDetDescrElement::m_r_raw
- */
- float r_raw() const;
- /** @brief cell dphi
- * @copydoc CaloGeoDetDescrElement::m_dphi
- */
- float dphi() const;
- /** @brief cell deta
- * @copydoc CaloGeoDetDescrElement::m_deta
- */
- float deta() const;
- /** @brief cell dr
- * @copydoc CaloGeoDetDescrElement::m_dr
- */
- float dr() const;
-
- /** @brief cell x
- * @copydoc CaloGeoDetDescrElement::m_x
- */
- float x() const;
- /** @brief cell y
- * @copydoc CaloGeoDetDescrElement::m_y
- */
- float y() const;
- /** @brief cell z
- * @copydoc CaloGeoDetDescrElement::m_z
- */
- float z() const;
- /** @brief cell x_raw
- * @copydoc CaloGeoDetDescrElement::m_x_raw
- */
- float x_raw() const;
- /** @brief cell y_raw
- * @copydoc CaloGeoDetDescrElement::m_y_raw
- */
- float y_raw() const;
- /** @brief cell z_raw
- * @copydoc CaloGeoDetDescrElement::m_z_raw
- */
- float z_raw() const;
- /** @brief cell dx
- * @copydoc CaloGeoDetDescrElement::m_dx
- */
- float dx() const;
- /** @brief cell dy
- * @copydoc CaloGeoDetDescrElement::m_dy
- */
- float dy() const;
- /** @brief cell dz
- * @copydoc CaloGeoDetDescrElement::m_dz
- */
- float dz() const;
-
- /** @brief cell identifier
- */
- Identifier identify() const;
-
- unsigned long long calo_hash() const;
-
- int getSampling() const ;
-
- private:
- //ACH protected:
- //
- long long m_identify;
- long long m_hash_id;
-
- int m_calosample;
-
- /** @brief cylindric coordinates : eta
- */
- float m_eta;
- /** @brief cylindric coordinates : phi
- */
- float m_phi;
-
- /** @brief this one is cached for algorithm working in transverse Energy
- */
- float m_sinTh;
- /** @brief this one is cached for algorithm working in transverse Energy
- */
- float m_cosTh;
-
- /** @brief cylindric coordinates : delta eta
- */
- float m_deta;
- /** @brief cylindric coordinates : delta phi
- */
- float m_dphi;
-
- /** @brief cylindric coordinates : r
- */
-
- float m_volume;
-
- /** @brief cache to allow fast px py pz computation
- */
- float m_sinPhi;
-
- /** @brief cache to allow fast px py pz computation
- */
- float m_cosPhi;
-
- /** @brief cylindric coordinates : r
- */
- float m_r;
- /** @brief cylindric coordinates : eta_raw
- */
- float m_eta_raw;
- /** @brief cylindric coordinates : phi_raw
- */
- float m_phi_raw;
- /** @brief cylindric coordinates : r_raw
- */
- float m_r_raw;
- /** @brief cylindric coordinates : delta r
- */
- float m_dr;
-
- /** @brief cartesian coordinates : X
- */
- float m_x;
- /** @brief cartesian coordinates : Y
- */
- float m_y;
- /** @brief cartesian coordinates : Z
- */
- float m_z;
- /** @brief cartesian coordinates : X raw
- */
- float m_x_raw;
- /** @brief cartesian coordinates : Y raw
- */
- float m_y_raw;
- /** @brief cartesian coordinates : Z raw
- */
- float m_z_raw;
- /** @brief cartesian coordinates : delta X
- */
- float m_dx;
- /** @brief cartesian coordinates : delta Y
- */
- float m_dy;
- /** @brief cartesian coordinates : delta Z
- */
- float m_dz;
-
-};
-
-inline Identifier CaloGeoDetDescrElement::identify() const
-{
- Identifier id((unsigned long long) m_identify);
- return id;
-}
-
-inline unsigned long long CaloGeoDetDescrElement::calo_hash() const
-{
- return m_hash_id;
-}
-
-inline int CaloGeoDetDescrElement::getSampling() const
-{ return m_calosample;}
-inline float CaloGeoDetDescrElement::eta() const
-{ return m_eta;}
-inline float CaloGeoDetDescrElement::phi() const
-{ return m_phi;}
-inline float CaloGeoDetDescrElement::r() const
-{ return m_r;}
-inline float CaloGeoDetDescrElement::eta_raw() const
-{ return m_eta_raw;}
-inline float CaloGeoDetDescrElement::phi_raw() const
-{ return m_phi_raw;}
-inline float CaloGeoDetDescrElement::r_raw() const
-{ return m_r_raw;}
-inline float CaloGeoDetDescrElement::deta() const
-{ return m_deta;}
-inline float CaloGeoDetDescrElement::dphi() const
-{ return m_dphi;}
-inline float CaloGeoDetDescrElement::dr() const
-{ return m_dr;}
-
-inline float CaloGeoDetDescrElement::x() const
-{ return m_x;}
-inline float CaloGeoDetDescrElement::y() const
-{ return m_y;}
-inline float CaloGeoDetDescrElement::z() const
-{ return m_z;}
-inline float CaloGeoDetDescrElement::x_raw() const
-{ return m_x_raw;}
-inline float CaloGeoDetDescrElement::y_raw() const
-{ return m_y_raw;}
-inline float CaloGeoDetDescrElement::z_raw() const
-{ return m_z_raw;}
-inline float CaloGeoDetDescrElement::dx() const
-{ return m_dx;}
-inline float CaloGeoDetDescrElement::dy() const
-{ return m_dy;}
-inline float CaloGeoDetDescrElement::dz() const
-{ return m_dz;}
-
-#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h
deleted file mode 100644
index 20cc3cd3c9c3..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h
+++ /dev/null
@@ -1,249 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef ISF_FASTCALOSIMPARAMETRIZATION_CALOGEODETDESCRELEMENT_H
-#define ISF_FASTCALOSIMPARAMETRIZATION_CALOGEODETDESCRELEMENT_H
-
-class CaloGeometry;
-
-#include "Identifier/Identifier.h"
-#include <iostream>
-#include <iomanip>
-#include <cmath>
-
-class CaloGeoDetDescrElement
-{
- friend class CaloGeometry;
- public:
- CaloGeoDetDescrElement() {
- m_identify = 0;
- m_hash_id = 0;
- m_calosample = 0;
- m_eta = 0;
- m_phi = 0;
- m_deta = 0;
- m_dphi = 0;
- m_r = 0;
- m_eta_raw = 0;
- m_phi_raw = 0;
- m_r_raw = 0;
- m_dr = 0;
- m_x = 0;
- m_y = 0;
- m_z = 0;
- m_x_raw = 0;
- m_y_raw = 0;
- m_z_raw = 0;
- m_dx = 0;
- m_dy = 0;
- m_dz = 0;
- };
-
- /** @brief virtual destructor
- */
- virtual ~CaloGeoDetDescrElement() {};
-
- /** @brief cell eta
- */
- float eta() const;
- /** @brief cell phi
- */
- float phi() const;
- /** @brief cell r
- */
- float r() const;
- /** @brief cell eta_raw
- */
- float eta_raw() const;
- /** @brief cell phi_raw
- */
- float phi_raw() const;
- /** @brief cell r_raw
- */
- float r_raw() const;
- /** @brief cell dphi
- */
- float dphi() const;
- /** @brief cell deta
- */
- float deta() const;
- /** @brief cell dr
- */
- float dr() const;
-
- /** @brief cell x
- */
- float x() const;
- /** @brief cell y
- */
- float y() const;
- /** @brief cell z
- */
- float z() const;
- /** @brief cell x_raw
- */
- float x_raw() const;
- /** @brief cell y_raw
- */
- float y_raw() const;
- /** @brief cell z_raw
- */
- float z_raw() const;
- /** @brief cell dx
- */
- float dx() const;
- /** @brief cell dy
- */
- float dy() const;
- /** @brief cell dz
- */
- float dz() const;
-
- /** @brief cell identifier
- */
- Identifier identify() const;
-
- unsigned long long calo_hash() const;
-
- int getSampling() const ;
-
- //ACH protected:
- //
- long long m_identify;
- long long m_hash_id;
-
- int m_calosample;
-
- /** @brief cylindric coordinates : eta
- */
- float m_eta;
- /** @brief cylindric coordinates : phi
- */
- float m_phi;
-
- /** @brief this one is cached for algorithm working in transverse Energy
- */
- float m_sinTh;
- /** @brief this one is cached for algorithm working in transverse Energy
- */
- float m_cosTh;
-
- /** @brief cylindric coordinates : delta eta
- */
- float m_deta;
- /** @brief cylindric coordinates : delta phi
- */
- float m_dphi;
-
- /** @brief cylindric coordinates : r
- */
-
- float m_volume;
-
- /** @brief cache to allow fast px py pz computation
- */
- float m_sinPhi;
-
- /** @brief cache to allow fast px py pz computation
- */
- float m_cosPhi;
-
- /** @brief cylindric coordinates : r
- */
- float m_r;
- /** @brief cylindric coordinates : eta_raw
- */
- float m_eta_raw;
- /** @brief cylindric coordinates : phi_raw
- */
- float m_phi_raw;
- /** @brief cylindric coordinates : r_raw
- */
- float m_r_raw;
- /** @brief cylindric coordinates : delta r
- */
- float m_dr;
-
- /** @brief cartesian coordinates : X
- */
- float m_x;
- /** @brief cartesian coordinates : Y
- */
- float m_y;
- /** @brief cartesian coordinates : Z
- */
- float m_z;
- /** @brief cartesian coordinates : X raw
- */
- float m_x_raw;
- /** @brief cartesian coordinates : Y raw
- */
- float m_y_raw;
- /** @brief cartesian coordinates : Z raw
- */
- float m_z_raw;
- /** @brief cartesian coordinates : delta X
- */
- float m_dx;
- /** @brief cartesian coordinates : delta Y
- */
- float m_dy;
- /** @brief cartesian coordinates : delta Z
- */
- float m_dz;
-
-};
-
-inline Identifier CaloGeoDetDescrElement::identify() const
-{
- Identifier id((unsigned long long) m_identify);
- return id;
-}
-
-inline unsigned long long CaloGeoDetDescrElement::calo_hash() const
-{
- return m_hash_id;
-}
-
-inline int CaloGeoDetDescrElement::getSampling() const
-{ return m_calosample;}
-inline float CaloGeoDetDescrElement::eta() const
-{ return m_eta;}
-inline float CaloGeoDetDescrElement::phi() const
-{ return m_phi;}
-inline float CaloGeoDetDescrElement::r() const
-{ return m_r;}
-inline float CaloGeoDetDescrElement::eta_raw() const
-{ return m_eta_raw;}
-inline float CaloGeoDetDescrElement::phi_raw() const
-{ return m_phi_raw;}
-inline float CaloGeoDetDescrElement::r_raw() const
-{ return m_r_raw;}
-inline float CaloGeoDetDescrElement::deta() const
-{ return m_deta;}
-inline float CaloGeoDetDescrElement::dphi() const
-{ return m_dphi;}
-inline float CaloGeoDetDescrElement::dr() const
-{ return m_dr;}
-
-inline float CaloGeoDetDescrElement::x() const
-{ return m_x;}
-inline float CaloGeoDetDescrElement::y() const
-{ return m_y;}
-inline float CaloGeoDetDescrElement::z() const
-{ return m_z;}
-inline float CaloGeoDetDescrElement::x_raw() const
-{ return m_x_raw;}
-inline float CaloGeoDetDescrElement::y_raw() const
-{ return m_y_raw;}
-inline float CaloGeoDetDescrElement::z_raw() const
-{ return m_z_raw;}
-inline float CaloGeoDetDescrElement::dx() const
-{ return m_dx;}
-inline float CaloGeoDetDescrElement::dy() const
-{ return m_dy;}
-inline float CaloGeoDetDescrElement::dz() const
-{ return m_dz;}
-
-#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoGeometry.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoGeometry.h
deleted file mode 100644
index 19adbd677982..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeoGeometry.h
+++ /dev/null
@@ -1,200 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef ISF_FASTCALOSIMPARAMETRIZATION_CALOGEOGEOMETRY_H
-#define ISF_FASTCALOSIMPARAMETRIZATION_CALOGEOGEOMETRY_H
-
-#include <TMath.h>
-
-#include <vector>
-#include <map>
-#include <iostream>
-
-#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
-#include "Identifier/Identifier.h"
-//#include "ISF_FastCaloSimParametrization/ICaloGeometry.h"
-#include "ISF_FastCaloSimParametrization/MeanAndRMS.h"
-#include "ISF_FastCaloSimParametrization/FSmap.h"
-#include "ISF_FastCaloSimParametrization/FCAL_ChannelMap.h"
-
-class CaloGeoDetDescrElement;
-class TCanvas;
-class TGraphErrors;
-
-typedef std::map< Identifier , const CaloGeoDetDescrElement* > t_cellmap;
-typedef std::map< double , const CaloGeoDetDescrElement* > t_eta_cellmap;
-
-class CaloGeoGeometryLookup {
- public:
- CaloGeoGeometryLookup(int ind=0);
- virtual ~CaloGeoGeometryLookup();
-
- bool IsCompatible(const CaloGeoDetDescrElement* cell);
- void add(const CaloGeoDetDescrElement* cell);
- t_cellmap::size_type size() const {return m_cells.size();};
- int index() const {return m_index;};
- void set_index(int ind) {m_index=ind;};
- void post_process();
- bool has_overlap(CaloGeoGeometryLookup* ref);
- void merge_into_ref(CaloGeoGeometryLookup* ref);
- //void CalculateTransformation();
-
- float mineta() const {return m_mineta;};
- float maxeta() const {return m_maxeta;};
- float minphi() const {return m_minphi;};
- float maxphi() const {return m_maxphi;};
-
- float mineta_raw() const {return m_mineta_raw;};
- float maxeta_raw() const {return m_maxeta_raw;};
- float minphi_raw() const {return m_minphi_raw;};
- float maxphi_raw() const {return m_maxphi_raw;};
-
- float minx() const {return m_mineta;};
- float maxx() const {return m_maxeta;};
- float miny() const {return m_minphi;};
- float maxy() const {return m_maxphi;};
-
- float minx_raw() const {return m_mineta_raw;};
- float maxx_raw() const {return m_maxeta_raw;};
- float miny_raw() const {return m_minphi_raw;};
- float maxy_raw() const {return m_maxphi_raw;};
-
- const MeanAndRMS& deta() {return m_deta;};
- const MeanAndRMS& dphi() {return m_dphi;};
- float mindeta() {return m_mindeta;};
- float mindphi() {return m_mindphi;};
- const MeanAndRMS& dx() {return m_deta;};
- const MeanAndRMS& dy() {return m_dphi;};
- float mindx() {return m_mindeta;};
- float mindy() {return m_mindphi;};
-
- const MeanAndRMS& eta_correction() {return m_eta_correction;};
- const MeanAndRMS& phi_correction() {return m_phi_correction;};
- const MeanAndRMS& x_correction() {return m_eta_correction;};
- const MeanAndRMS& y_correction() {return m_phi_correction;};
-
- int cell_grid_eta() const {return m_cell_grid_eta;};
- int cell_grid_phi() const {return m_cell_grid_phi;};
- void set_xy_grid_adjustment_factor(float factor) {m_xy_grid_adjustment_factor=factor;};
-
- virtual const CaloGeoDetDescrElement* getDDE(float eta,float phi,float* distance=0,int* steps=0);
-
- protected:
- float neta_double() {return (maxeta_raw()-mineta_raw())/deta().mean();};
- float nphi_double() {return (maxphi_raw()-minphi_raw())/dphi().mean();};
- Int_t neta() {return TMath::Nint( neta_double() );};
- Int_t nphi() {return TMath::Nint( nphi_double() );};
-
- //FCal is not sufficiently regular to use submaps with regular mapping
- float nx_double() {return (maxx_raw()-minx_raw())/mindx();};
- float ny_double() {return (maxy_raw()-miny_raw())/mindy();};
- Int_t nx() {return TMath::Nint(TMath::Ceil( nx_double() ));};
- Int_t ny() {return TMath::Nint(TMath::Ceil( ny_double() ));};
-
- float m_xy_grid_adjustment_factor;
-
- int raw_eta_position_to_index(float eta_raw) const {return TMath::Floor((eta_raw-mineta_raw())/m_deta_double);};
- int raw_phi_position_to_index(float phi_raw) const {return TMath::Floor((phi_raw-minphi_raw())/m_dphi_double);};
- bool index_range_adjust(int& ieta,int& iphi);
- float calculate_distance_eta_phi(const CaloGeoDetDescrElement* DDE,float eta,float phi,float& dist_eta0,float& dist_phi0);
-
- int m_index;
- t_cellmap m_cells;
- std::vector< std::vector< const CaloGeoDetDescrElement* > > m_cell_grid;
- int m_cell_grid_eta,m_cell_grid_phi;
- float m_mineta,m_maxeta,m_minphi,m_maxphi;
- float m_mineta_raw,m_maxeta_raw,m_minphi_raw,m_maxphi_raw;
- float m_mineta_correction,m_maxeta_correction,m_minphi_correction,m_maxphi_correction;
-
- MeanAndRMS m_deta,m_dphi,m_eta_correction,m_phi_correction;
- float m_mindeta,m_mindphi;
- float m_deta_double,m_dphi_double;
-};
-
-
-class CaloGeoGeometry //: virtual public ICaloGeometry {
-{
- public :
- static const int MAX_SAMPLING;
-
- static Identifier m_debug_identify;
- static bool m_debug;
-
- CaloGeoGeometry();
- virtual ~CaloGeoGeometry();
-
- virtual bool PostProcessGeometry();
-
- virtual void Validate();
-
- virtual const CaloGeoDetDescrElement* getDDE(Identifier identify);
- virtual const CaloGeoDetDescrElement* getDDE(int sampling, Identifier identify);
-
- virtual const CaloGeoDetDescrElement* getDDE(int sampling,float eta,float phi,float* distance=0,int* steps=0);
- virtual const CaloGeoDetDescrElement* getFCalDDE(int sampling,float eta,float phi,float z);
-
- double deta(int sample,double eta) const;
- void minmaxeta(int sample,double eta,double& mineta,double& maxeta) const;
- double rzmid(int sample,double eta) const;
- double rzent(int sample,double eta) const;
- double rzext(int sample,double eta) const;
- double rmid(int sample,double eta) const;
- double rent(int sample,double eta) const;
- double rext(int sample,double eta) const;
- double zmid(int sample,double eta) const;
- double zent(int sample,double eta) const;
- double zext(int sample,double eta) const;
- double rpos(int sample,double eta,int subpos = CaloSubPos::SUBPOS_MID) const;
- double zpos(int sample,double eta,int subpos = CaloSubPos::SUBPOS_MID) const;
- double rzpos(int sample,double eta,int subpos = CaloSubPos::SUBPOS_MID) const;
- bool isCaloBarrel(int sample) const {return m_isCaloBarrel[sample];};
- static std::string SamplingName(int sample);
-
- TGraphErrors* GetGraph(unsigned int sample) const {return m_graph_layers[sample];};
- void SetDoGraphs(bool dographs=true) {m_dographs=dographs;};
- bool DoGraphs() const {return m_dographs;};
-
- TCanvas* DrawGeoForPhi0();
- FCAL_ChannelMap* GetFCAL_ChannelMap(){return &m_FCal_ChannelMap;}
-
- protected:
- virtual void addcell(const CaloGeoDetDescrElement* cell);
-
- virtual void post_process(int layer);
-
- virtual void PrintMapInfo(int i, int j);
-
- virtual void InitRZmaps();
-
- t_cellmap m_cells;
- std::vector< t_cellmap > m_cells_in_sampling;
- std::vector< t_eta_cellmap > m_cells_in_sampling_for_phi0;
- std::vector< std::vector< CaloGeoGeometryLookup* > > m_cells_in_regions;
-
- std::vector< bool > m_isCaloBarrel;
- std::vector< double > m_min_eta_sample[2]; //[side][calosample]
- std::vector< double > m_max_eta_sample[2]; //[side][calosample]
- std::vector< FSmap< double , double > > m_rmid_map[2]; //[side][calosample]
- std::vector< FSmap< double , double > > m_zmid_map[2]; //[side][calosample]
- std::vector< FSmap< double , double > > m_rent_map[2]; //[side][calosample]
- std::vector< FSmap< double , double > > m_zent_map[2]; //[side][calosample]
- std::vector< FSmap< double , double > > m_rext_map[2]; //[side][calosample]
- std::vector< FSmap< double , double > > m_zext_map[2]; //[side][calosample]
-
- bool m_dographs;
- std::vector< TGraphErrors* > m_graph_layers;
- FCAL_ChannelMap m_FCal_ChannelMap;
- /*
- double m_min_eta_sample[2][MAX_SAMPLING]; //[side][calosample]
- double m_max_eta_sample[2][MAX_SAMPLING]; //[side][calosample]
- FSmap< double , double > m_rmid_map[2][MAX_SAMPLING]; //[side][calosample]
- FSmap< double , double > m_zmid_map[2][MAX_SAMPLING]; //[side][calosample]
- FSmap< double , double > m_rent_map[2][MAX_SAMPLING]; //[side][calosample]
- FSmap< double , double > m_zent_map[2][MAX_SAMPLING]; //[side][calosample]
- FSmap< double , double > m_rext_map[2][MAX_SAMPLING]; //[side][calosample]
- FSmap< double , double > m_zext_map[2][MAX_SAMPLING]; //[side][calosample]
- */
-};
-
-#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h
index 3ba69de1560d..b0adf9aae342 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometry.h
@@ -12,104 +12,15 @@
#include <iostream>
#include "ISF_FastCaloSimParametrization/ICaloGeometry.h"
+#include "ISF_FastCaloSimParametrization/CaloGeometryLookup.h"
#include "ISF_FastCaloSimParametrization/MeanAndRMS.h"
#include "ISF_FastCaloSimParametrization/FSmap.h"
-#include "ISF_FastCaloSimParametrization/FCAL_ChannelMap.h"
+#include "LArReadoutGeometry/FCAL_ChannelMap.h"
class CaloDetDescrElement;
class TCanvas;
class TGraphErrors;
-typedef std::map< Identifier , const CaloDetDescrElement* > t_cellmap;
-typedef std::map< double , const CaloDetDescrElement* > t_eta_cellmap;
-
-class CaloGeometryLookup {
- public:
- CaloGeometryLookup(int ind=0);
- virtual ~CaloGeometryLookup();
-
- bool IsCompatible(const CaloDetDescrElement* cell);
- void add(const CaloDetDescrElement* cell);
- t_cellmap::size_type size() const {return m_cells.size();};
- int index() const {return m_index;};
- void set_index(int ind) {m_index=ind;};
- void post_process();
- bool has_overlap(CaloGeometryLookup* ref);
- void merge_into_ref(CaloGeometryLookup* ref);
- //void CalculateTransformation();
-
- float mineta() const {return m_mineta;};
- float maxeta() const {return m_maxeta;};
- float minphi() const {return m_minphi;};
- float maxphi() const {return m_maxphi;};
-
- float mineta_raw() const {return m_mineta_raw;};
- float maxeta_raw() const {return m_maxeta_raw;};
- float minphi_raw() const {return m_minphi_raw;};
- float maxphi_raw() const {return m_maxphi_raw;};
-
- float minx() const {return m_mineta;};
- float maxx() const {return m_maxeta;};
- float miny() const {return m_minphi;};
- float maxy() const {return m_maxphi;};
-
- float minx_raw() const {return m_mineta_raw;};
- float maxx_raw() const {return m_maxeta_raw;};
- float miny_raw() const {return m_minphi_raw;};
- float maxy_raw() const {return m_maxphi_raw;};
-
- const MeanAndRMS& deta() {return m_deta;};
- const MeanAndRMS& dphi() {return m_dphi;};
- float mindeta() {return m_mindeta;};
- float mindphi() {return m_mindphi;};
- const MeanAndRMS& dx() {return m_deta;};
- const MeanAndRMS& dy() {return m_dphi;};
- float mindx() {return m_mindeta;};
- float mindy() {return m_mindphi;};
-
- const MeanAndRMS& eta_correction() {return m_eta_correction;};
- const MeanAndRMS& phi_correction() {return m_phi_correction;};
- const MeanAndRMS& x_correction() {return m_eta_correction;};
- const MeanAndRMS& y_correction() {return m_phi_correction;};
-
- int cell_grid_eta() const {return m_cell_grid_eta;};
- int cell_grid_phi() const {return m_cell_grid_phi;};
- void set_xy_grid_adjustment_factor(float factor) {m_xy_grid_adjustment_factor=factor;};
-
- virtual const CaloDetDescrElement* getDDE(float eta,float phi,float* distance=0,int* steps=0);
-
- protected:
- float neta_double() {return (maxeta_raw()-mineta_raw())/deta().mean();};
- float nphi_double() {return (maxphi_raw()-minphi_raw())/dphi().mean();};
- Int_t neta() {return TMath::Nint( neta_double() );};
- Int_t nphi() {return TMath::Nint( nphi_double() );};
-
- //FCal is not sufficiently regular to use submaps with regular mapping
- float nx_double() {return (maxx_raw()-minx_raw())/mindx();};
- float ny_double() {return (maxy_raw()-miny_raw())/mindy();};
- Int_t nx() {return TMath::Nint(TMath::Ceil( nx_double() ));};
- Int_t ny() {return TMath::Nint(TMath::Ceil( ny_double() ));};
-
- float m_xy_grid_adjustment_factor;
-
- int raw_eta_position_to_index(float eta_raw) const {return TMath::Floor((eta_raw-mineta_raw())/m_deta_double);};
- int raw_phi_position_to_index(float phi_raw) const {return TMath::Floor((phi_raw-minphi_raw())/m_dphi_double);};
- bool index_range_adjust(int& ieta,int& iphi);
- float calculate_distance_eta_phi(const CaloDetDescrElement* DDE,float eta,float phi,float& dist_eta0,float& dist_phi0);
-
- int m_index;
- t_cellmap m_cells;
- std::vector< std::vector< const CaloDetDescrElement* > > m_cell_grid;
- int m_cell_grid_eta,m_cell_grid_phi;
- float m_mineta,m_maxeta,m_minphi,m_maxphi;
- float m_mineta_raw,m_maxeta_raw,m_minphi_raw,m_maxphi_raw;
- float m_mineta_correction,m_maxeta_correction,m_minphi_correction,m_maxphi_correction;
-
- MeanAndRMS m_deta,m_dphi,m_eta_correction,m_phi_correction;
- float m_mindeta,m_mindphi;
- float m_deta_double,m_dphi_double;
-};
-
class CaloGeometry : virtual public ICaloGeometry {
public :
static const int MAX_SAMPLING;
@@ -153,6 +64,7 @@ class CaloGeometry : virtual public ICaloGeometry {
TCanvas* DrawGeoForPhi0();
FCAL_ChannelMap* GetFCAL_ChannelMap(){return &m_FCal_ChannelMap;}
+ virtual void LoadFCalGeometryFromFiles(TString filename1,TString filename2,TString filename3); // Initialize m_FCal_ChannelMap
protected:
virtual void addcell(const CaloDetDescrElement* cell);
@@ -180,7 +92,8 @@ class CaloGeometry : virtual public ICaloGeometry {
bool m_dographs;
std::vector< TGraphErrors* > m_graph_layers;
- FCAL_ChannelMap m_FCal_ChannelMap;
+ FCAL_ChannelMap m_FCal_ChannelMap; // for hit-to-cell assignment in FCal
+
/*
double m_min_eta_sample[2][MAX_SAMPLING]; //[side][calosample]
double m_max_eta_sample[2][MAX_SAMPLING]; //[side][calosample]
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometryFromFile.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometryFromFile.h
deleted file mode 100644
index 5238fa4d78db..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/CaloGeometryFromFile.h
+++ /dev/null
@@ -1,22 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef CaloGeometryFromFile_h
-#define CaloGeometryFromFile_h
-
-#include "ISF_FastCaloSimParametrization/CaloGeoGeometry.h"
-
-
-class CaloGeometryFromFile:public CaloGeoGeometry {
-public :
- CaloGeometryFromFile();
- virtual ~CaloGeometryFromFile();
-
- virtual bool LoadGeometryFromFile(TString filename,TString treename);
- virtual bool LoadFCalGeometryFromFiles(TString filename1,TString filename2,TString filename3);
- void DrawFCalGraph(int isam,int color);
-};
-
-#endif
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/FCAL_ChannelMap.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/FCAL_ChannelMap.h
deleted file mode 100755
index 3a7cb07efbbf..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/FCAL_ChannelMap.h
+++ /dev/null
@@ -1,286 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-// ***************************************************************************
-// Liquid Argon FCAL detector description package
-// -----------------------------------------
-//
-//
-// 10-Sep-2000 S.Simion Handling of the FCAL read-out identifiers
-// Jan-2001 R.Sobie Modify for persistency
-// Feb-2002 R.Sobie Use same FCAL geometry files as simulation
-// ***************************************************************************
-
-#ifndef FCAL_CHANNELMAP_H
-#define FCAL_CHANNELMAP_H
-
-
-//<<<<<< INCLUDES >>>>>>
-
-
-#include <math.h>
-#include <vector>
-#include <functional>
-#include <map>
-#include <string>
-#include <stdexcept>
-
-
-/** This class contains the tube and tile maps for the FCAL <br>
- * A tile is of a set of FCAL tubes
- */
-class FCAL_ChannelMap
-{
-public:
-
-
-
- //
- // Typedefs:
- //
- typedef unsigned int tileName_t;
- typedef unsigned int tubeID_t;
-
- //------------ TubePosition: X,Y Position of an FCAL Tube---- //
- // //
- class TubePosition { //
- public: //
- TubePosition(); //
- TubePosition(tileName_t name, float x, float y, std::string hvFT);//
- tileName_t get_tileName() const; //
- float x() const; //
- float y() const; //
- std::string getHVft() const; //Gabe
- private: //
- tileName_t m_tileName; //
- float m_x; //
- float m_y; //
- std::string m_hvFT; // Gabe
- }; //
- //-------------------------------------------------------------//
-
-
-
- /** Constructors: */
- FCAL_ChannelMap( int itemp);
-
- typedef std::map<tubeID_t, TubePosition > tubeMap_t;
- typedef tubeMap_t::size_type tubemap_sizetype;
- typedef tubeMap_t::value_type tubemap_valuetype;
- typedef tubeMap_t::const_iterator tubemap_const_iterator;
-
- /** tubeMap access functions */
- tubemap_const_iterator tubemap_begin (int isam) const;
- tubemap_const_iterator tubemap_end (int isam) const;
- tubemap_sizetype tubemap_size (int isam) const;
-
-
- // Get the tube by copy number:
- tubemap_const_iterator getTubeByCopyNumber ( int isam, int copyNo) const;
-
-
- /** ---- For the new LArFCAL_ID Identifier
- */
-
- bool getTileID(int isam,
- float x,
- float y,
- int& eta,
- int& phi) const throw (std::range_error);
-
- /** For reconstruction, decoding of tile identifiers */
- float x(int isam,
- int eta,
- int phi) const throw(std::range_error) ;
- float y(int isam,
- int eta,
- int phi) const throw(std::range_error) ;
-
- void tileSize(int sam, int eta, int phi,
- float& dx, float& dy) const throw(std::range_error) ;
-
- void tileSize(int isam, int ntubes, float& dx, float& dy) const;
-
-
- /** print functions */
- void print_tubemap (int isam) const;
-
- /** set and get for the inversion flags
- */
-
- bool invert_x() const;
- bool invert_xy() const ;
- void set_invert_x(bool );
- void set_invert_xy(bool );
-
-
- // Fill this. The information comes from the database.
- void add_tube (const std::string & tileName, int mod, int id, int i, int j, double xCm, double yCm);//original
- void add_tube (const std::string & tileName, int mod, int id, int i, int j, double xCm, double yCm, std::string hvFT);//29-03-07 include HV
-
-
- // Finish the job. Create the tile map.
- void finish();
-
-
- class TilePosition {
- public:
- TilePosition();
- TilePosition(float x, float y, int ntubes);
- float x() const;
- float y() const;
- unsigned int ntubes() const;
- private:
- float m_x;
- float m_y;
- unsigned int m_ntubes;
- };
-
- /** TileMap */
- typedef std::map<tileName_t, TilePosition > tileMap_t;
- typedef tileMap_t::size_type tileMap_sizetype;
- typedef tileMap_t::value_type tileMap_valuetype;
- typedef tileMap_t::const_iterator tileMap_const_iterator;
-
- // Methods to iterate over the tile map:
- tileMap_const_iterator begin(int isam) const {return m_tileMap[isam-1].begin();}
- tileMap_const_iterator end(int isam) const {return m_tileMap[isam-1].end(); };
-
-
-private:
- /** Geometrical parameters here, in CLHEP::cm please to be compatible with G3 */
- static const double m_tubeSpacing[];
- double m_tubeDx[3];
- double m_tubeDy[3];
- double m_tileDx[3];
- double m_tileDy[3];
- mutable bool m_invert_x; // Some geometry need x inverted
- mutable bool m_invert_xy; // Some geometry need xy crossed
-
- tileMap_t m_tileMap[3];
- void create_tileMap( int isam );
-
- /** TubeMap */
- tubeMap_t m_tubeMap[3];
- std::vector<tubemap_const_iterator> m_tubeIndex[3];
-};
-
-
-// Tube Position functions
-
-inline
-FCAL_ChannelMap::TubePosition::TubePosition()
- :
- m_tileName(0),
- m_x(0),
- m_y(0),
- m_hvFT("")
-{}
-
-inline
-FCAL_ChannelMap::TubePosition::TubePosition(tileName_t name, float x, float y, std::string hvFT)
- :
- m_tileName(name),
- m_x(x),
- m_y(y),
- m_hvFT(hvFT)
-{}
-
-
-inline FCAL_ChannelMap::tileName_t
-FCAL_ChannelMap::TubePosition::get_tileName() const
-{
- return m_tileName;
-}
-
-
-inline float
-FCAL_ChannelMap::TubePosition::x() const
-{
- return m_x;
-}
-
-inline float
-FCAL_ChannelMap::TubePosition::y() const
-{
- return m_y;
-}
-
-
-inline std::string
-FCAL_ChannelMap::TubePosition::getHVft() const
-{
- return m_hvFT;
-}
-
-
-// Tile Position functions
-
-inline
-FCAL_ChannelMap::TilePosition::TilePosition()
- :
- m_x(0),
- m_y(0),
- m_ntubes(0)
-{}
-
-inline
-FCAL_ChannelMap::TilePosition::TilePosition(float x, float y, int ntubes)
- :
- m_x(x),
- m_y(y),
- m_ntubes(ntubes)
-{}
-
-
-inline float
-FCAL_ChannelMap::TilePosition::x() const
-{
- return m_x;
-}
-
-inline float
-FCAL_ChannelMap::TilePosition::y() const
-{
- return m_y;
-}
-
-inline unsigned int
-FCAL_ChannelMap::TilePosition::ntubes() const
-{
- return m_ntubes;
-}
-
-inline
-bool FCAL_ChannelMap::invert_x() const {
-return m_invert_x;
-}
-
-inline
-bool FCAL_ChannelMap::invert_xy() const {
-return m_invert_xy;
-}
-
-inline
-void FCAL_ChannelMap::set_invert_x(bool flag ){
-m_invert_x = flag;
-return ;
-}
-
-inline
-void FCAL_ChannelMap::set_invert_xy(bool flag ){
-m_invert_xy = flag;
-return ;
-}
-
-//#include "CLIDSvc/CLASS_DEF.h"
-//CLASS_DEF(FCAL_ChannelMap, 74242524, 1)
-
-#endif // LARDETDESCR_FCAL_CHANNELMAP_H
-
-
-
-
-
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ICaloGeometry.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ICaloGeometry.h
index 7ac58a1d9e6b..d302c71e1de7 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ICaloGeometry.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/ICaloGeometry.h
@@ -8,7 +8,6 @@
#include "Identifier/Identifier.h"
#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
#include "CaloDetDescr/CaloDetDescrElement.h"
-#include "ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h"
class CaloDetDescrElement;
class ICaloGeometry {
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 924341e2462b..e4c345953e7d 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
@@ -160,6 +160,23 @@ class ISF_HitAnalysis : public AthAlgorithm {
double m_ptruth_p;
int m_pdgid;
+ /*
+ std::vector<std::vector<float> >* m_TTC_entrance_eta;
+ std::vector<std::vector<float> >* m_TTC_entrance_phi;
+ std::vector<std::vector<float> >* m_TTC_entrance_r;
+ std::vector<std::vector<float> >* m_TTC_entrance_z;
+ std::vector<std::vector<float> >* m_TTC_back_eta;
+ std::vector<std::vector<float> >* m_TTC_back_phi;
+ std::vector<std::vector<float> >* m_TTC_back_r;
+ std::vector<std::vector<float> >* m_TTC_back_z;
+ std::vector<float>* m_TTC_IDCaloBoundary_eta;
+ std::vector<float>* m_TTC_IDCaloBoundary_phi;
+ std::vector<float>* m_TTC_IDCaloBoundary_r;
+ std::vector<float>* m_TTC_IDCaloBoundary_z;
+ std::vector<float>* m_TTC_Angle3D;
+ std::vector<float>* m_TTC_AngleEta;
+ */
+
std::vector<std::vector<float> >* m_newTTC_entrance_eta;
std::vector<std::vector<float> >* m_newTTC_entrance_phi;
std::vector<std::vector<float> >* m_newTTC_entrance_r;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/python/ISF_FastCaloSimParametrizationConfig.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/python/ISF_FastCaloSimParametrizationConfig.py
index e69e3b183b57..76281ea3bedc 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/python/ISF_FastCaloSimParametrizationConfig.py
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/python/ISF_FastCaloSimParametrizationConfig.py
@@ -1,6 +1,5 @@
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-
"""
Tools configurations for ISF_FastCaloSimParametrization
"""
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiPostInclude.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiPostInclude.py
index 3491c5354028..55cef205defc 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiPostInclude.py
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiPostInclude.py
@@ -1,3 +1,3 @@
-streamRDO.ItemList+=["ISF_FCS_Parametrization::FCS_StepInfoCollection#MergedEventSteps","LArHitContainer#*","TileHitVector#*"];
+streamRDO.ItemList+=["ISF_FCS_Parametrization::FCS_StepInfoCollection#ZHMergedEventSteps","LArHitContainer#*","TileHitVector#*"];
from AthenaCommon.CfgGetter import getPublicTool
getPublicTool("LArPileUpTool").CrossTalk=False
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiTilePostInclude.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiTilePostInclude.py
index 64d3c37f2e08..81eef137a23c 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiTilePostInclude.py
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_DigiTilePostInclude.py
@@ -1,6 +1,7 @@
from AthenaCommon.AppMgr import ServiceMgr
-
-ToolSvc.TileHitVecToCntTool.HitTimeFlag=1
+ServiceMgr.MessageSvc.OutputLevel= VERBOSE
+#rel 19 is configuring this via a tool!!
+ToolSvc.TileHitVecToCntTool.HitTimeFlag=1
ServiceMgr.TileInfoLoader.NPhElec=0
ServiceMgr.TileInfoLoader.NPhElecA=0
@@ -8,14 +9,13 @@ ServiceMgr.TileInfoLoader.NPhElecC=0
ServiceMgr.TileInfoLoader.NPhElecD=0
ServiceMgr.TileInfoLoader.NPhElecE=0
-topSequence.TileDigitsMaker.IntegerDigits=True
-
###########
print ServiceMgr.TileInfoLoader
print topSequence.TileDigitsMaker
print ToolSvc.TileHitVecToCntTool
-#topSequence.TileDigitsMaker.OutputLevel=1
-#TileHitVecToCntTool.OutputLevel=1
-#ToolSvc.TileRawChannelBuilderOptATLAS.OutputLevel=1
-#ServiceMgr.MessageSvc.OutputLevel= VERBOSE
-###########
+topSequence.TileDigitsMaker.OutputLevel=1
+topSequence.TileDigitsMaker.IntegerDigits=True
+TileHitVecToCntTool.OutputLevel=1
+ToolSvc.TileRawChannelBuilderOptATLAS.OutputLevel=1
+###########
+
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude_1mm.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude_1mm.py
new file mode 100644
index 000000000000..3965eb8bb970
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude_1mm.py
@@ -0,0 +1,5 @@
+from G4AtlasApps.SimFlags import simFlags
+simFlags.OptionalUserActionList.addAction('G4UA::FastCaloSimParamActionTool',['BeginOfEvent','EndOfEvent','BeginOfRun','EndOfRun','Step'])
+simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"shift_lar_subhit",True)
+simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"shorten_lar_step",True)
+simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"MaxRadiusLAr",1)
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeoGeometry.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeoGeometry.cxx
deleted file mode 100644
index 92d40e7dc165..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeoGeometry.cxx
+++ /dev/null
@@ -1,1220 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#include "ISF_FastCaloSimParametrization/CaloGeoGeometry.h"
-#include <TTree.h>
-#include <TVector2.h>
-#include <TRandom.h>
-#include <TCanvas.h>
-#include <TH2D.h>
-#include <TGraphErrors.h>
-#include <TVector3.h>
-#include <TLegend.h>
-
-//#include "CaloDetDescr/CaloGeoDetDescrElement.h"
-#include "ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h"
-#include "CaloGeoHelpers/CaloSampling.h"
-#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
-//#include "TMVA/Tools.h"
-//#include "TMVA/Factory.h"
-
-using namespace std;
-
-const int CaloGeoGeometry::MAX_SAMPLING = CaloCell_ID_FCS::MaxSample; //number of calorimeter layers/samplings
-
-Identifier CaloGeoGeometry::m_debug_identify;
-bool CaloGeoGeometry::m_debug=false;
-
-CaloGeoGeometryLookup::CaloGeoGeometryLookup(int ind):m_xy_grid_adjustment_factor(0.75),m_index(ind)
-{
- m_mineta=+10000;
- m_maxeta=-10000;
- m_minphi=+10000;
- m_maxphi=-10000;
-
- m_mineta_raw=+10000;
- m_maxeta_raw=-10000;
- m_minphi_raw=+10000;
- m_maxphi_raw=-10000;
-
- m_mindeta=10000;
- m_mindphi=10000;
-
- m_mineta_correction=+10000;
- m_maxeta_correction=-10000;
- m_minphi_correction=+10000;
- m_maxphi_correction=-10000;
-
- m_cell_grid_eta=0.;
- m_cell_grid_phi=0.;
- m_deta_double =0.;
- m_dphi_double =0.;
-}
-
-CaloGeoGeometryLookup::~CaloGeoGeometryLookup()
-{
-}
-
-bool CaloGeoGeometryLookup::has_overlap(CaloGeoGeometryLookup* ref)
-{
- if(m_cells.size()==0) return false;
- for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
- const CaloGeoDetDescrElement* cell=ic->second;
- if(ref->IsCompatible(cell)) return true;
- }
- return false;
-}
-
-void CaloGeoGeometryLookup::merge_into_ref(CaloGeoGeometryLookup* ref)
-{
- for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
- const CaloGeoDetDescrElement* cell=ic->second;
- ref->add(cell);
- }
-}
-
-
-bool CaloGeoGeometryLookup::IsCompatible(const CaloGeoDetDescrElement* cell)
-{
- if(m_cells.size()==0) return true;
- t_cellmap::iterator ic=m_cells.begin();
- const CaloGeoDetDescrElement* refcell=ic->second;
- int sampling=refcell->getSampling();
- if(cell->getSampling()!=sampling) return false;
- if(cell->eta_raw()*refcell->eta_raw()<0) return false;
- if(sampling<21) { // keep only cells reasonable close in eta to each other
- if(TMath::Abs(cell->deta()-refcell->deta())>0.001) return false;
- if(TMath::Abs(cell->dphi()-refcell->dphi())>0.001) return false;
-
- if( cell->eta()<mineta()-2.1*cell->deta() ) return false;
- if( cell->eta()>maxeta()+2.1*cell->deta() ) return false;
-
- double delta_eta=(cell->eta_raw()-refcell->eta_raw())/cell->deta();
- double delta_phi=(cell->phi_raw()-refcell->phi_raw())/cell->dphi();
- if(TMath::Abs(delta_eta-TMath::Nint(delta_eta)) > 0.01) return false;
- if(TMath::Abs(delta_phi-TMath::Nint(delta_phi)) > 0.01) return false;
-
- //cout<<"Compatible: s="<<cell->getSampling()<<"~"<<refcell->getSampling()<<"; ";
- //cout<<"eta="<<cell->eta_raw()<<","<<refcell->eta_raw()<<"; ";
- //cout<<"phi="<<cell->phi_raw()<<","<<refcell->phi_raw()<<"; ";
- //cout<<"deta="<<cell->deta()<<"~"<<refcell->deta()<<"; ";
- //cout<<"dphi="<<cell->dphi()<<"~"<<refcell->dphi()<<";";
- //cout<<"mineta="<<mineta()<<", maxeta="<<maxeta()<<"; ";
- //cout<<endl;
-
- } else {
- //FCal is not sufficiently regular to use submaps with regular mapping
- return true;
-
-
- //if(TMath::Abs(cell->dx()-refcell->dx())>0.001) return false;
- //if(TMath::Abs(cell->dy()-refcell->dy())>0.001) return false;
- //if( cell->x()<minx()-20*cell->dx() ) return false;
- //if( cell->x()>maxx()+20*cell->dx() ) return false;
- //if( cell->y()<miny()-20*cell->dy() ) return false;
- //if( cell->y()>maxy()+20*cell->dy() ) return false;
-
- //double delta_x=(cell->x_raw()-refcell->x_raw())/cell->dx();
- //double delta_y=(cell->y_raw()-refcell->y_raw())/cell->dy();
- //if(TMath::Abs(delta_x-TMath::Nint(delta_x)) > 0.01) return false;
- //if(TMath::Abs(delta_y-TMath::Nint(delta_y)) > 0.01) return false;
-
- }
-
- return true;
-}
-
-void CaloGeoGeometryLookup::add(const CaloGeoDetDescrElement* cell)
-{
- if(cell->getSampling()<21) {
- m_deta.add(cell->deta());
- m_dphi.add(cell->dphi());
- m_mindeta=TMath::Min(cell->deta(),m_mindeta);
- m_mindphi=TMath::Min(cell->dphi(),m_mindphi);
-
- m_eta_correction.add(cell->eta_raw()-cell->eta());
- m_phi_correction.add(cell->phi_raw()-cell->phi());
- m_mineta_correction=TMath::Min(cell->eta_raw()-cell->eta() , m_mineta_correction);
- m_maxeta_correction=TMath::Max(cell->eta_raw()-cell->eta() , m_maxeta_correction);
- m_minphi_correction=TMath::Min(cell->phi_raw()-cell->phi() , m_minphi_correction);
- m_maxphi_correction=TMath::Max(cell->phi_raw()-cell->phi() , m_maxphi_correction);
-
- m_mineta=TMath::Min(cell->eta()-cell->deta()/2,m_mineta);
- m_maxeta=TMath::Max(cell->eta()+cell->deta()/2,m_maxeta);
- m_minphi=TMath::Min(cell->phi()-cell->dphi()/2,m_minphi);
- m_maxphi=TMath::Max(cell->phi()+cell->dphi()/2,m_maxphi);
-
- m_mineta_raw=TMath::Min(cell->eta_raw()-cell->deta()/2,m_mineta_raw);
- m_maxeta_raw=TMath::Max(cell->eta_raw()+cell->deta()/2,m_maxeta_raw);
- m_minphi_raw=TMath::Min(cell->phi_raw()-cell->dphi()/2,m_minphi_raw);
- m_maxphi_raw=TMath::Max(cell->phi_raw()+cell->dphi()/2,m_maxphi_raw);
- } else {
- m_deta.add(cell->dx());
- m_dphi.add(cell->dy());
- m_mindeta=TMath::Min(cell->dx(),m_mindeta);
- m_mindphi=TMath::Min(cell->dy(),m_mindphi);
-
- m_eta_correction.add(cell->x_raw()-cell->x());
- m_phi_correction.add(cell->y_raw()-cell->y());
- m_mineta_correction=TMath::Min(cell->x_raw()-cell->x() , m_mineta_correction);
- m_maxeta_correction=TMath::Max(cell->x_raw()-cell->x() , m_maxeta_correction);
- m_minphi_correction=TMath::Min(cell->y_raw()-cell->y() , m_minphi_correction);
- m_maxphi_correction=TMath::Max(cell->y_raw()-cell->y() , m_maxphi_correction);
-
- m_mineta=TMath::Min(cell->x()-cell->dx()/2,m_mineta);
- m_maxeta=TMath::Max(cell->x()+cell->dx()/2,m_maxeta);
- m_minphi=TMath::Min(cell->y()-cell->dy()/2,m_minphi);
- m_maxphi=TMath::Max(cell->y()+cell->dy()/2,m_maxphi);
-
- m_mineta_raw=TMath::Min(cell->x_raw()-cell->dx()/2,m_mineta_raw);
- m_maxeta_raw=TMath::Max(cell->x_raw()+cell->dx()/2,m_maxeta_raw);
- m_minphi_raw=TMath::Min(cell->y_raw()-cell->dy()/2,m_minphi_raw);
- m_maxphi_raw=TMath::Max(cell->y_raw()+cell->dy()/2,m_maxphi_raw);
- }
- m_cells[cell->identify()]=cell;
-}
-
-bool CaloGeoGeometryLookup::index_range_adjust(int& ieta,int& iphi)
-{
- while(iphi<0) {iphi+=m_cell_grid_phi;};
- while(iphi>=m_cell_grid_phi) {iphi-=m_cell_grid_phi;};
- if(ieta<0) {
- ieta=0;
- return false;
- }
- if(ieta>=m_cell_grid_eta) {
- ieta=m_cell_grid_eta-1;
- return false;
- }
- return true;
-}
-
-void CaloGeoGeometryLookup::post_process()
-{
- if(size()==0) return;
- t_cellmap::iterator ic=m_cells.begin();
- const CaloGeoDetDescrElement* refcell=ic->second;
- int sampling=refcell->getSampling();
- if(sampling<21) {
- double rneta=neta_double()-neta();
- double rnphi=nphi_double()-nphi();
- if(TMath::Abs(rneta)>0.05 || TMath::Abs(rnphi)>0.05) {
- cout<<"ERROR: can't map cells into a regular grid for sampling "<<sampling<<", map index "<<index()<<endl;
- return;
- }
- m_cell_grid_eta=neta();
- m_cell_grid_phi=nphi();
- m_deta_double=deta();
- m_dphi_double=dphi();
- } else {
- //double rnx=nx_double()-nx();
- //double rny=ny_double()-ny();
- //if(TMath::Abs(rnx)>0.05 || TMath::Abs(rny)>0.05) {
- // cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": mapping cells into a regular grid, although cells don't fit well"<<endl;
- //}
-
- m_cell_grid_eta=TMath::Nint(TMath::Ceil(nx_double()/m_xy_grid_adjustment_factor));
- m_cell_grid_phi=TMath::Nint(TMath::Ceil(ny_double()/m_xy_grid_adjustment_factor));
- m_deta_double=mindx()*m_xy_grid_adjustment_factor;
- m_dphi_double=mindy()*m_xy_grid_adjustment_factor;
- }
-
- m_cell_grid.resize(m_cell_grid_eta);
- for(int ieta=0;ieta<m_cell_grid_eta;++ieta) {
- m_cell_grid[ieta].resize(m_cell_grid_phi,0);
- }
-
- for(ic=m_cells.begin();ic!=m_cells.end();++ic) {
- refcell=ic->second;
- int ieta,iphi;
- if(sampling<21) {
- ieta=raw_eta_position_to_index(refcell->eta_raw());
- iphi=raw_phi_position_to_index(refcell->phi_raw());
- } else {
- ieta=raw_eta_position_to_index(refcell->x_raw());
- iphi=raw_phi_position_to_index(refcell->y_raw());
- }
- if(index_range_adjust(ieta,iphi)) {
- if(m_cell_grid[ieta][iphi]) {
- cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": Already cell found at pos ("<<ieta<<","<<iphi<<"): id="
- <<m_cell_grid[ieta][iphi]->identify()<<"->"<<refcell->identify()<<endl;
- cout<<" x="<<m_cells[m_cell_grid[ieta][iphi]->identify()]->x_raw()<<" -> "<<refcell->x_raw();
- cout<<" , y="<<m_cells[m_cell_grid[ieta][iphi]->identify()]->y_raw()<<" -> "<<refcell->y_raw();
- cout<<" mindx="<<m_deta_double<<" mindy="<<m_dphi_double<<endl;
- cout<<endl;
- }
- m_cell_grid[ieta][iphi]=refcell;
- } else {
- cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": Cell at pos ("<<ieta<<","<<iphi<<"): id="
- <<refcell->identify()<<" outside eta range"<<endl;
- }
- }
-
- int ncells=0;
- int nempty=0;
- for(int ieta=0;ieta<m_cell_grid_eta;++ieta) {
- for(int iphi=0;iphi<m_cell_grid_phi;++iphi) {
- if(!m_cell_grid[ieta][iphi]) {
- ++nempty;
- //cout<<"Sampling "<<sampling<<"_"<<index()<<": No cell at pos ("<<ieta<<","<<iphi<<")"<<endl;
- } else {
- ++ncells;
- }
- }
- }
- // cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": "<<ncells<<"/"<<size()<<" cells filled, "<<nempty<<" empty grid positions deta="<<m_deta_double<<" dphi="<<m_dphi_double<<endl;
-}
-
-float CaloGeoGeometryLookup::calculate_distance_eta_phi(const CaloGeoDetDescrElement* DDE,float eta,float phi,float& dist_eta0,float& dist_phi0)
-{
- dist_eta0=(eta - DDE->eta())/m_deta_double;
- dist_phi0=(TVector2::Phi_mpi_pi(phi - DDE->phi()))/m_dphi_double;
- float abs_dist_eta0=TMath::Abs(dist_eta0);
- float abs_dist_phi0=TMath::Abs(dist_phi0);
- return TMath::Max(abs_dist_eta0,abs_dist_phi0)-0.5;
-}
-
-
-const CaloGeoDetDescrElement* CaloGeoGeometryLookup::getDDE(float eta,float phi,float* distance,int* steps)
-{
- float dist;
- const CaloGeoDetDescrElement* bestDDE=0;
- if(!distance) distance=&dist;
- *distance=+10000000;
- int intsteps=0;
- if(!steps) steps=&intsteps;
-
- float best_eta_corr=m_eta_correction;
- float best_phi_corr=m_phi_correction;
-
- float raw_eta=eta+best_eta_corr;
- float raw_phi=phi+best_phi_corr;
-
- int ieta=raw_eta_position_to_index(raw_eta);
- int iphi=raw_phi_position_to_index(raw_phi);
- index_range_adjust(ieta,iphi);
-
- const CaloGeoDetDescrElement* newDDE=m_cell_grid[ieta][iphi];
- float bestdist=+10000000;
- ++(*steps);
- int nsearch=0;
- while(newDDE && nsearch<3) {
- float dist_eta0,dist_phi0;
- *distance=calculate_distance_eta_phi(newDDE,eta,phi,dist_eta0,dist_phi0);
- bestDDE=newDDE;
- bestdist=*distance;
-
- if(CaloGeoGeometry::m_debug || CaloGeoGeometry::m_debug_identify==newDDE->identify()) {
- cout<<"CaloGeoGeometryLookup::getDDE, search iteration "<<nsearch<<endl;
- cout<<"cell id="<<newDDE->identify()<<" steps "<<*steps<<" steps, eta="<<eta<<" phi="<<phi<<" dist="<<*distance<<" cell_grid_eta="<<cell_grid_eta()<<" cell_grid_phi="<<cell_grid_phi()<<" m_deta_double="<<m_deta_double<<" m_dphi_double="<<m_dphi_double<<" 1st_eta_corr="<<best_eta_corr<<" 1st_phi_corr="<<best_phi_corr<<endl;
- cout<<" input sampling="<<newDDE->getSampling()<<" eta="<<newDDE->eta()<<" eta_raw="<<newDDE->eta_raw()<<" deta="<<newDDE->deta()<<" ("<<(newDDE->eta_raw()-newDDE->eta())/newDDE->deta()<<") phi="<<newDDE->phi()<<" phi_raw="<<newDDE->phi_raw()<<" dphi="<<newDDE->dphi()<<" ("<<(newDDE->phi_raw()-newDDE->phi())/newDDE->dphi()<<")"<<endl;
- cout<<" ieta="<<ieta<<" iphi="<<iphi<<endl;
- cout<<" dist_eta0="<<dist_eta0<<" ,dist_phi0="<<dist_phi0<<endl;
- }
-
- if(*distance<0) return newDDE;
- //correct ieta and iphi by the observed difference to the hit cell
- ieta+=TMath::Nint(dist_eta0);
- iphi+=TMath::Nint(dist_phi0);
- index_range_adjust(ieta,iphi);
- const CaloGeoDetDescrElement* oldDDE=newDDE;
- newDDE=m_cell_grid[ieta][iphi];
- ++(*steps);
- ++nsearch;
- if(oldDDE==newDDE) break;
- }
- if(CaloGeoGeometry::m_debug && !newDDE) {
- cout<<"CaloGeoGeometryLookup::getDDE, direct search ieta="<<ieta<<" iphi="<<iphi<<" is empty"<<endl;
- }
- float minieta=ieta+TMath::Nint(TMath::Floor(m_mineta_correction/cell_grid_eta()));
- float maxieta=ieta+TMath::Nint(TMath::Ceil(m_maxeta_correction/cell_grid_eta()));
- float miniphi=iphi+TMath::Nint(TMath::Floor(m_minphi_correction/cell_grid_phi()));
- float maxiphi=iphi+TMath::Nint(TMath::Ceil(m_maxphi_correction/cell_grid_phi()));
- if(minieta<0) minieta=0;
- if(maxieta>=m_cell_grid_eta) maxieta=m_cell_grid_eta-1;
- for(int iieta=minieta;iieta<=maxieta;++iieta) {
- for(int iiphi=miniphi;iiphi<=maxiphi;++iiphi) {
- ieta=iieta;
- iphi=iiphi;
- index_range_adjust(ieta,iphi);
- newDDE=m_cell_grid[ieta][iphi];
- ++(*steps);
- if(newDDE) {
- float dist_eta0,dist_phi0;
- *distance=calculate_distance_eta_phi(newDDE,eta,phi,dist_eta0,dist_phi0);
-
- if(CaloGeoGeometry::m_debug || CaloGeoGeometry::m_debug_identify==newDDE->identify()) {
- cout<<"CaloGeoGeometryLookup::getDDE, windows search! minieta="<<m_mineta_correction/cell_grid_eta()<<" maxieta="<<m_maxeta_correction/cell_grid_eta()<<" miniphi="<<m_minphi_correction/cell_grid_phi()<<" maxiphi="<<m_maxphi_correction/cell_grid_phi()<<endl;
- cout<<"cell id="<<newDDE->identify()<<" steps "<<*steps<<" steps, eta="<<eta<<" phi="<<phi<<" dist="<<*distance<<endl;
- cout<<" input sampling="<<newDDE->getSampling()<<" eta="<<newDDE->eta()<<" eta_raw="<<newDDE->eta_raw()<<" deta="<<newDDE->deta()<<" ("<<(newDDE->eta_raw()-newDDE->eta())/newDDE->deta()<<") phi="<<newDDE->phi()<<" phi_raw="<<newDDE->phi_raw()<<" dphi="<<newDDE->dphi()<<" ("<<(newDDE->phi_raw()-newDDE->phi())/newDDE->dphi()<<")"<<endl;
- cout<<" ieta="<<ieta<<" iphi="<<iphi<<endl;
- cout<<" dist_eta0="<<dist_eta0<<" ,dist_phi0="<<dist_phi0<<endl;
- }
-
- if(*distance<0) return newDDE;
- if(*distance<bestdist) {
- bestDDE=newDDE;
- bestdist=*distance;
- }
- } else if(CaloGeoGeometry::m_debug) {
- cout<<"CaloGeoGeometryLookup::getDDE, windows search ieta="<<ieta<<" iphi="<<iphi<<" is empty"<<endl;
- }
- }
- }
- *distance=bestdist;
- return bestDDE;
-}
-
-
-
-/*
-void CaloGeometryLookup::CalculateTransformation()
-{
- gROOT->cd();
-
- TTree* tmap = new TTree( "mapping" , "mapping" );
-
- float eta,phi,Deta_raw,Dphi_raw;
- tmap->Branch("eta", &eta,"eta/F");
- tmap->Branch("phi", &phi,"phi/F");
- tmap->Branch("Deta_raw", &Deta_raw,"Deta_raw/F");
- tmap->Branch("Dphi_raw", &Dphi_raw,"Dphi_raw/F");
-
- if(m_cells.size()==0) return;
-
- int sampling=0;
- for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
- CaloGeoDetDescrElement* refcell=ic->second;
- sampling=refcell->getSampling();
- if(sampling<21) {
- eta=refcell->eta();
- phi=refcell->phi();
- Deta_raw=refcell->eta_raw()-eta;
- Dphi_raw=refcell->phi_raw()-phi;
- } else {
- eta=refcell->x();
- phi=refcell->y();
- Deta_raw=refcell->x_raw()-eta;
- Dphi_raw=refcell->y_raw()-phi;
- }
- tmap->Fill();
- tmap->Fill(); //Fill twice to have all events and training and test tree
- }
- tmap->Print();
-
- TString outfileName( Form("Mapping%d_%d.root",sampling,m_index) );
- TFile* outputFile = new TFile( outfileName, "RECREATE" );
- //TFile* outputFile = new TMemFile( outfileName, "RECREATE" );
-
- TMVA::Factory *factory = new TMVA::Factory( Form("Mapping%d_%d.root",sampling,m_index) , outputFile, "!V:!Silent:Color:DrawProgressBar" );
-
- factory->AddVariable( "eta", "calo eta", "1", 'F' );
- factory->AddVariable( "phi", "calo phi", "1", 'F' );
- factory->AddTarget( "Deta_raw" );
- factory->AddTarget( "Dphi_raw" );
-
- factory->AddRegressionTree( tmap );
-
- //factory->PrepareTrainingAndTestTree( "" , Form("nTrain_Regression=%d:nTest_Regression=%d:SplitMode=Random:NormMode=NumEvents:!V",(int)m_cells.size(),(int)m_cells.size()) );
- factory->PrepareTrainingAndTestTree( "" , "nTrain_Regression=0:nTest_Regression=0:SplitMode=Alternate:NormMode=NumEvents:!V" );
-
- factory->BookMethod( TMVA::Types::kMLP, "MLP", "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+5:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator" );
-
- cout<<"=== Start trainging ==="<<endl;
- // Train MVAs using the set of training events
- factory->TrainAllMethods();
-
- cout<<"=== Start testing ==="<<endl;
- // ---- Evaluate all MVAs using the set of test events
- factory->TestAllMethods();
-
- cout<<"=== Start evaluation ==="<<endl;
- // ----- Evaluate and compare performance of all configured MVAs
- factory->EvaluateAllMethods();
-
- outputFile->Close();
-
- delete factory;
-
- delete tmap;
-}
-*/
-
-CaloGeoGeometry::CaloGeoGeometry() : m_cells_in_sampling(MAX_SAMPLING),m_cells_in_sampling_for_phi0(MAX_SAMPLING),m_cells_in_regions(MAX_SAMPLING),m_isCaloBarrel(MAX_SAMPLING),m_dographs(false),m_FCal_ChannelMap(0)
-{
- //TMVA::Tools::Instance();
- for(int i=0;i<2;++i) {
- m_min_eta_sample[i].resize(MAX_SAMPLING); //[side][calosample]
- m_max_eta_sample[i].resize(MAX_SAMPLING); //[side][calosample]
- m_rmid_map[i].resize(MAX_SAMPLING); //[side][calosample]
- m_zmid_map[i].resize(MAX_SAMPLING); //[side][calosample]
- m_rent_map[i].resize(MAX_SAMPLING); //[side][calosample]
- m_zent_map[i].resize(MAX_SAMPLING); //[side][calosample]
- m_rext_map[i].resize(MAX_SAMPLING); //[side][calosample]
- m_zext_map[i].resize(MAX_SAMPLING); //[side][calosample]
- }
- m_graph_layers.resize(MAX_SAMPLING);
- for(int i=CaloCell_ID_FCS::FirstSample;i<CaloCell_ID_FCS::MaxSample;++i) {
- m_graph_layers[i]=0;
- CaloSampling::CaloSample s=static_cast<CaloSampling::CaloSample>(i);
- m_isCaloBarrel[i]=(CaloSampling::barrelPattern() & CaloSampling::getSamplingPattern(s))!=0;
- }
- m_isCaloBarrel[CaloCell_ID_FCS::TileGap3]=false;
-}
-
-CaloGeoGeometry::~CaloGeoGeometry()
-{
-}
-
-void CaloGeoGeometry::addcell(const CaloGeoDetDescrElement* cell)
-{
- int sampling=cell->getSampling();
- Identifier identify=cell->identify();
-
- //m_cells[identify]=cell;
- //m_cells_in_sampling[sampling][identify]=cell;
-
- m_cells[identify]= new CaloGeoDetDescrElement(*cell);
- m_cells_in_sampling[sampling][identify]= new CaloGeoDetDescrElement(*cell);
-
- CaloGeoGeometryLookup* lookup=0;
- for(unsigned int i=0;i<m_cells_in_regions[sampling].size();++i) {
- if(m_cells_in_regions[sampling][i]->IsCompatible(cell)) {
- lookup=m_cells_in_regions[sampling][i];
- //cout<<sampling<<": found compatible map"<<endl;
- break;
- }
- }
- if(!lookup) {
- lookup=new CaloGeoGeometryLookup(m_cells_in_regions[sampling].size());
- m_cells_in_regions[sampling].push_back(lookup);
- }
- lookup->add(cell);
-}
-
-void CaloGeoGeometry::PrintMapInfo(int i, int j)
-{
- cout<<" map "<<j<<": "<<m_cells_in_regions[i][j]->size()<<" cells";
- if(i<21) {
- cout<<", "<<m_cells_in_regions[i][j]->cell_grid_eta()<<"*"<<m_cells_in_regions[i][j]->cell_grid_phi();
- cout<<", deta="<<m_cells_in_regions[i][j]->deta().mean()<<"+-"<<m_cells_in_regions[i][j]->deta().rms();
- cout<<", dphi="<<m_cells_in_regions[i][j]->dphi().mean()<<"+-"<<m_cells_in_regions[i][j]->dphi().rms();
- cout<<", mineta="<<m_cells_in_regions[i][j]->mineta()<<", maxeta="<<m_cells_in_regions[i][j]->maxeta();
- cout<<", minphi="<<m_cells_in_regions[i][j]->minphi()<<", maxphi="<<m_cells_in_regions[i][j]->maxphi();
- cout<<endl<<" ";
- cout<<", mineta_raw="<<m_cells_in_regions[i][j]->mineta_raw()<<", maxeta_raw="<<m_cells_in_regions[i][j]->maxeta_raw();
- cout<<", minphi_raw="<<m_cells_in_regions[i][j]->minphi_raw()<<", maxphi_raw="<<m_cells_in_regions[i][j]->maxphi_raw();
- cout<<endl;
- } else {
- cout<<", "<<m_cells_in_regions[i][j]->cell_grid_eta()<<"*"<<m_cells_in_regions[i][j]->cell_grid_phi();
- cout<<", dx="<<m_cells_in_regions[i][j]->dx().mean()<<"+-"<<m_cells_in_regions[i][j]->dx().rms();
- cout<<", dy="<<m_cells_in_regions[i][j]->dy().mean()<<"+-"<<m_cells_in_regions[i][j]->dy().rms();
- cout<<", mindx="<<m_cells_in_regions[i][j]->mindx();
- cout<<", mindy="<<m_cells_in_regions[i][j]->mindy();
- cout<<", minx="<<m_cells_in_regions[i][j]->minx()<<", maxx="<<m_cells_in_regions[i][j]->maxx();
- cout<<", miny="<<m_cells_in_regions[i][j]->miny()<<", maxy="<<m_cells_in_regions[i][j]->maxy();
- cout<<endl<<" ";
- cout<<", minx_raw="<<m_cells_in_regions[i][j]->minx_raw()<<", maxx_raw="<<m_cells_in_regions[i][j]->maxx_raw();
- cout<<", miny_raw="<<m_cells_in_regions[i][j]->miny_raw()<<", maxy_raw="<<m_cells_in_regions[i][j]->maxy_raw();
- cout<<endl;
- }
-}
-
-void CaloGeoGeometry::post_process(int sampling)
-{
- //cout<<"post processing sampling "<<sampling<<endl;
- bool found_overlap=false;
- for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
- /*
- cout<<"Sample "<<sampling<<": checking map "<<j<<"/"<<m_cells_in_regions[sampling].size();
- for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
- cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
- }
- cout<<endl;
- */
- for(unsigned int i=j+1;i<m_cells_in_regions[sampling].size();++i) {
- if(m_cells_in_regions[sampling][i]->has_overlap(m_cells_in_regions[sampling][j])) {
- if(!found_overlap) {
- cout<<"Sample "<<sampling<<", starting maps : "<<m_cells_in_regions[sampling].size();
- for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
- cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
- }
- cout<<endl;
- }
- found_overlap=true;
- /*
- cout<<"Sample "<<sampling<<": Found overlap between map "<<j<<" and "<<i<<", "
- <<m_cells_in_regions[sampling].size()<<" total maps"<<endl;
- cout<<" current configuration map "<<j<<"/"<<m_cells_in_regions[sampling].size();
- for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
- cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
- }
- cout<<endl;
-
- PrintMapInfo(sampling,j);
- PrintMapInfo(sampling,i);
- */
-
- CaloGeoGeometryLookup* toremove=m_cells_in_regions[sampling][i];
- toremove->merge_into_ref(m_cells_in_regions[sampling][j]);
-
- /*
- cout<<" New ";
- PrintMapInfo(sampling,j);
- */
-
- for(unsigned int k=i;k<m_cells_in_regions[sampling].size()-1;++k) {
- m_cells_in_regions[sampling][k]=m_cells_in_regions[sampling][k+1];
- m_cells_in_regions[sampling][k]->set_index(k);
- }
- m_cells_in_regions[sampling].resize(m_cells_in_regions[sampling].size()-1);
-
- /*
- cout<<" new configuration map "<<j<<"/"<<m_cells_in_regions[sampling].size();
- for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
- cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
- }
- cout<<endl;
- */
-
- --i;
- }
- }
- }
-
- if(found_overlap) {
- cout<<"Sample "<<sampling<<", final maps : "<<m_cells_in_regions[sampling].size();
- for(unsigned int k=0;k<m_cells_in_regions[sampling].size();++k) {
- cout<<", "<<k<<":"<<m_cells_in_regions[sampling][k]->size();
- }
- cout<<endl;
- }
-
- if(found_overlap) {
- cout<<"Run another round of ";
- post_process(sampling);
- }
-}
-
-void CaloGeoGeometry::InitRZmaps()
-{
- std::cout<<"CHECKME InitRZmaps"<<std::endl;
-
- int nok=0;
-
- FSmap< double , double > rz_map_eta [2][MAX_SAMPLING];
- FSmap< double , double > rz_map_rmid[2][MAX_SAMPLING];
- FSmap< double , double > rz_map_zmid[2][MAX_SAMPLING];
- FSmap< double , double > rz_map_rent[2][MAX_SAMPLING];
- FSmap< double , double > rz_map_zent[2][MAX_SAMPLING];
- FSmap< double , double > rz_map_rext[2][MAX_SAMPLING];
- FSmap< double , double > rz_map_zext[2][MAX_SAMPLING];
- FSmap< double , int > rz_map_n [2][MAX_SAMPLING];
-
-
- for(int side=0;side<=1;++side) for(int sample=0;sample<MAX_SAMPLING;++sample)
- {
- m_min_eta_sample[side][sample]=+1000;
- m_max_eta_sample[side][sample]=-1000;
- }
-
-
- for(t_cellmap::iterator calo_iter=m_cells.begin();calo_iter!=m_cells.end();++calo_iter)
- {
- const CaloGeoDetDescrElement* theDDE=(*calo_iter).second;
-
- if(theDDE)
- {
- ++nok;
- int sample=theDDE->getSampling();
-
- int side=0;
- int sign_side=-1;
- double eta_raw=theDDE->eta_raw();
- if(eta_raw>0) {
- side=1;
- sign_side=+1;
- }
-
- if(!m_cells_in_sampling_for_phi0[sample][eta_raw]) {
- m_cells_in_sampling_for_phi0[sample][eta_raw]=theDDE;
- } else {
- if(TMath::Abs(theDDE->phi()) < TMath::Abs(m_cells_in_sampling_for_phi0[sample][eta_raw]->phi())) {
- m_cells_in_sampling_for_phi0[sample][eta_raw]=theDDE;
- }
- }
-
- double min_eta=theDDE->eta()-theDDE->deta()/2;
- double max_eta=theDDE->eta()+theDDE->deta()/2;
- if(min_eta<m_min_eta_sample[side][sample]) m_min_eta_sample[side][sample]=min_eta;
- if(max_eta>m_max_eta_sample[side][sample]) m_max_eta_sample[side][sample]=max_eta;
-
- if(rz_map_eta[side][sample].find(eta_raw)==rz_map_eta[side][sample].end()) {
- rz_map_eta [side][sample][eta_raw]=0;
- rz_map_rmid[side][sample][eta_raw]=0;
- rz_map_zmid[side][sample][eta_raw]=0;
- rz_map_rent[side][sample][eta_raw]=0;
- rz_map_zent[side][sample][eta_raw]=0;
- rz_map_rext[side][sample][eta_raw]=0;
- rz_map_zext[side][sample][eta_raw]=0;
- rz_map_n [side][sample][eta_raw]=0;
- }
- rz_map_eta [side][sample][eta_raw]+=theDDE->eta();
- rz_map_rmid[side][sample][eta_raw]+=theDDE->r();
- rz_map_zmid[side][sample][eta_raw]+=theDDE->z();
- double drh=theDDE->dr()/2;
- double dzh=theDDE->dz();
- if(sample>=CaloSampling::PreSamplerB && sample<=CaloSampling::EMB3) {
- drh=theDDE->dr();
- }
- rz_map_rent[side][sample][eta_raw]+=theDDE->r()-drh;
- rz_map_zent[side][sample][eta_raw]+=theDDE->z()-dzh*sign_side;
- rz_map_rext[side][sample][eta_raw]+=theDDE->r()+drh;
- rz_map_zext[side][sample][eta_raw]+=theDDE->z()+dzh*sign_side;
- rz_map_n [side][sample][eta_raw]++;
-
- }
- }
-
- for(int side=0;side<=1;++side)
- {
- for(int sample=0;sample<MAX_SAMPLING;++sample)
- {
-
- if(rz_map_n[side][sample].size()>0)
- {
- for(FSmap< double , int >::iterator iter=rz_map_n[side][sample].begin();iter!=rz_map_n[side][sample].end();++iter)
- {
- double eta_raw=iter->first;
- if(iter->second<1)
- {
- //ATH_MSG_WARNING("rz-map for side="<<side<<" sample="<<sample<<" eta_raw="<<eta_raw<<" : #cells="<<iter->second<<" !!!");
- }
- else
- {
- double eta =rz_map_eta[side][sample][eta_raw]/iter->second;
- double rmid=rz_map_rmid[side][sample][eta_raw]/iter->second;
- double zmid=rz_map_zmid[side][sample][eta_raw]/iter->second;
- double rent=rz_map_rent[side][sample][eta_raw]/iter->second;
- double zent=rz_map_zent[side][sample][eta_raw]/iter->second;
- double rext=rz_map_rext[side][sample][eta_raw]/iter->second;
- double zext=rz_map_zext[side][sample][eta_raw]/iter->second;
-
- m_rmid_map[side][sample][eta]=rmid;
- m_zmid_map[side][sample][eta]=zmid;
-
-
-
- m_rent_map[side][sample][eta]=rent;
-
-
- m_zent_map[side][sample][eta]=zent;
- m_rext_map[side][sample][eta]=rext;
- m_zext_map[side][sample][eta]=zext;
- }
- }
- //ATH_MSG_DEBUG("rz-map for side="<<side<<" sample="<<sample<<" #etas="<<m_rmid_map[side][sample].size());
- }
- else
- {
- std::cout<<"rz-map for side="<<side<<" sample="<<sample<<" is empty!!!"<<std::endl;
- }
- } //sample
- } //side
-
- if(DoGraphs()) {
- int calocol[24]={1,2,3,4, // LAr barrel
- 1,2,3,4, // LAr EM endcap
- 1,2,3,4, // Hadronic end cap cal.
- 1,2,3, // Tile barrel
- 6,28,42, // Tile gap (ITC & scint)
- 1,2,3, // Tile extended barrel
- 1,2,3 // Forward EM endcap
- };
-
- for(int sample=0;sample<MAX_SAMPLING;++sample) {
- m_graph_layers[sample]=new TGraphErrors(rz_map_n[0][sample].size()+rz_map_n[1][sample].size());
- m_graph_layers[sample]->SetMarkerColor(calocol[sample]);
- m_graph_layers[sample]->SetLineColor(calocol[sample]);
- int np=0;
- for(int side=0;side<=1;++side) {
- for(FSmap< double , int >::iterator iter=rz_map_n[side][sample].begin();iter!=rz_map_n[side][sample].end();++iter) {
- double eta_raw=iter->first;
- int sign_side=-1;
- if(eta_raw>0) sign_side=+1;
- //double eta =rz_map_eta[side][sample][eta_raw]/iter->second;
- double rmid=rz_map_rmid[side][sample][eta_raw]/iter->second;
- double zmid=rz_map_zmid[side][sample][eta_raw]/iter->second;
- //double rent=rz_map_rent[side][sample][eta_raw]/iter->second;
- //double zent=rz_map_zent[side][sample][eta_raw]/iter->second;
- double rext=rz_map_rext[side][sample][eta_raw]/iter->second;
- double zext=rz_map_zext[side][sample][eta_raw]/iter->second;
- m_graph_layers[sample]->SetPoint(np,zmid,rmid);
- /*
- if(isCaloBarrel(sample)) {
- m_graph_layers[sample]->SetPointError(np,0,rext-rmid);
- } else {
- m_graph_layers[sample]->SetPointError(np,(zext-zent)*sign_side,0);
- }
- */
- m_graph_layers[sample]->SetPointError(np,(zext-zmid)*sign_side,rext-rmid);
- ++np;
- }
- }
- }
- }
-}
-
-TCanvas* CaloGeoGeometry::DrawGeoForPhi0()
-{
- TCanvas* c=new TCanvas("CaloGeoForPhi0","Calo geometry for #phi~0");
- TH2D* hcalolayout=new TH2D("hcalolayoutPhi0","Reconstruction geometry: calorimeter layout;z [mm];r [mm]",50,-7000,7000,50,0,4000);
- hcalolayout->Draw();
- hcalolayout->SetStats(0);
- hcalolayout->GetYaxis()->SetTitleOffset(1.4);
-
- int calocol[MAX_SAMPLING]={1,2,3,4, // LAr barrel
- 1,2,3,4, // LAr EM endcap
- 1,2,3,4, // Hadronic end cap cal.
- 1,2,3, // Tile barrel
- -42,-28,-6, // Tile gap (ITC & scint)
- 1,2,3, // Tile extended barrel
- 1,2,3 // Forward EM endcap
- };
-
- TLegend* leg=new TLegend(0.30,0.13,0.70,0.37);
- leg->SetFillStyle(0);
- leg->SetFillColor(10);
- leg->SetBorderSize(1);
- leg->SetNColumns(2);
-
- for(int sample=0;sample<MAX_SAMPLING;++sample) {
-// for(int sample=21;sample<22;++sample) {
- cout<<"Start sample "<<sample<<" ("<<SamplingName(sample)<<")"<<endl;
- int ngr=0;
- for(t_eta_cellmap::iterator calo_iter=m_cells_in_sampling_for_phi0[sample].begin();calo_iter!=m_cells_in_sampling_for_phi0[sample].end();++calo_iter) {
- const CaloGeoDetDescrElement* theDDE=(*calo_iter).second;
- if(theDDE) {
- TVector3 cv;
- TGraph* gr=new TGraph(5);
- gr->SetLineColor(TMath::Abs(calocol[sample]));
- gr->SetFillColor(TMath::Abs(calocol[sample]));
- if(calocol[sample]<0) {
- gr->SetFillStyle(1001);
- } else {
- gr->SetFillStyle(0);
- }
- gr->SetLineWidth(2);
- double r=theDDE->r();
- double dr=theDDE->dr();
- double x=theDDE->x();
- double dx=theDDE->dx();
- double y=theDDE->y();
- double dy=theDDE->dy();
- double z=theDDE->z();
- double dz=theDDE->dz()*2;
- double eta=theDDE->eta();
- double deta=theDDE->deta();
- if(CaloSampling::PreSamplerB<=sample && sample<=CaloSampling::EMB3) {
- dr*=2;
- }
-
- if(isCaloBarrel(sample)) {
- cv.SetPtEtaPhi(r-dr/2,eta-deta/2,0);
- gr->SetPoint(0,cv.Z(),cv.Pt());
- gr->SetPoint(4,cv.Z(),cv.Pt());
- cv.SetPtEtaPhi(r-dr/2,eta+deta/2,0);
- gr->SetPoint(1,cv.Z(),cv.Pt());
- cv.SetPtEtaPhi(r+dr/2,eta+deta/2,0);
- gr->SetPoint(2,cv.Z(),cv.Pt());
- cv.SetPtEtaPhi(r+dr/2,eta-deta/2,0);
- gr->SetPoint(3,cv.Z(),cv.Pt());
- } else {
- if(sample<CaloSampling::FCAL0) {
- cv.SetPtEtaPhi(1,eta-deta/2,0);cv*=(z-dz/2)/cv.Z();
- gr->SetPoint(0,cv.Z(),cv.Pt());
- gr->SetPoint(4,cv.Z(),cv.Pt());
- cv.SetPtEtaPhi(1,eta+deta/2,0);cv*=(z-dz/2)/cv.Z();
- gr->SetPoint(1,cv.Z(),cv.Pt());
- cv.SetPtEtaPhi(1,eta+deta/2,0);cv*=(z+dz/2)/cv.Z();
- gr->SetPoint(2,cv.Z(),cv.Pt());
- cv.SetPtEtaPhi(1,eta-deta/2,0);cv*=(z+dz/2)/cv.Z();
- gr->SetPoint(3,cv.Z(),cv.Pt());
- } else {
- double minr=r;
- double maxr=r;
- for(double px=x-dx/2;px<=x+dx/2;px+=dx) {
- for(double py=y-dy/2;py<=y+dy/2;py+=dy) {
- double pr=TMath::Sqrt(px*px+py*py);
- minr=TMath::Min(minr,pr);
- maxr=TMath::Max(maxr,pr);
- }
- }
- cv.SetXYZ(minr,0,z-dz/2);
- gr->SetPoint(0,cv.Z(),cv.Pt());
- gr->SetPoint(4,cv.Z(),cv.Pt());
- cv.SetXYZ(maxr,0,z-dz/2);
- gr->SetPoint(1,cv.Z(),cv.Pt());
- cv.SetXYZ(maxr,0,z+dz/2);
- gr->SetPoint(2,cv.Z(),cv.Pt());
- cv.SetXYZ(minr,0,z+dz/2);
- gr->SetPoint(3,cv.Z(),cv.Pt());
- }
- }
- //if(calocol[sample]>0) gr->Draw("Lsame");
- // else gr->Draw("LFsame");
- gr->Draw("LFsame");
- if(ngr==0) {
- std::string sname=Form("Sampling %2d : ",sample);
- sname+=SamplingName(sample);
- leg->AddEntry(gr,sname.c_str(),"LF");
- }
- ++ngr;
- }
- }
- cout<<"Done sample "<<sample<<" ("<<SamplingName(sample)<<")="<<ngr<<endl;
- }
- leg->Draw();
- return c;
-}
-
-const CaloGeoDetDescrElement* CaloGeoGeometry::getDDE(Identifier identify)
-{
- return m_cells[identify];
-}
-const CaloGeoDetDescrElement* CaloGeoGeometry::getDDE(int sampling,Identifier identify)
-{
- return m_cells_in_sampling[sampling][identify];
-}
-
-const CaloGeoDetDescrElement* CaloGeoGeometry::getDDE(int sampling,float eta,float phi,float* distance,int* steps)
-{
- if(sampling<0) return 0;
- if(sampling>=MAX_SAMPLING) return 0;
- if(m_cells_in_regions[sampling].size()==0) return 0;
-
- float dist;
- const CaloGeoDetDescrElement* bestDDE=0;
- if(!distance) distance=&dist;
- *distance=+10000000;
- int intsteps;
- int beststeps;
- if(steps) beststeps=(*steps);
- else beststeps=0;
-
- if(sampling<21) {
- for(int skip_range_check=0;skip_range_check<=1;++skip_range_check) {
- for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
- if(!skip_range_check) {
- if(eta<m_cells_in_regions[sampling][j]->mineta()) continue;
- if(eta>m_cells_in_regions[sampling][j]->maxeta()) continue;
- }
- if(steps) intsteps=(*steps);
- else intsteps=0;
- if(m_debug) {
- cout<<"CaloGeoGeometry::getDDE : check map"<<j<<" skip_range_check="<<skip_range_check<<endl;
- }
- float newdist;
- const CaloGeoDetDescrElement* newDDE=m_cells_in_regions[sampling][j]->getDDE(eta,phi,&newdist,&intsteps);
- if(m_debug) {
- cout<<"CaloGeoGeometry::getDDE : map"<<j<<" dist="<<newdist<<" best dist="<<*distance<<" steps="<<intsteps<<endl;
- }
- if(newdist<*distance) {
- bestDDE=newDDE;
- *distance=newdist;
- if(steps) beststeps=intsteps;
- if(newdist<-0.1) break; //stop, we are well within the hit cell
- }
- }
- if(bestDDE) break;
- }
- } else {
- return 0;
- //for(int skip_range_check=0;skip_range_check<=1;++skip_range_check) {
- //for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
- //if(!skip_range_check) {
- //if(eta<m_cells_in_regions[sampling][j]->minx()) continue;
- //if(eta>m_cells_in_regions[sampling][j]->maxx()) continue;
- //if(phi<m_cells_in_regions[sampling][j]->miny()) continue;
- //if(phi>m_cells_in_regions[sampling][j]->maxy()) continue;
- //}
- //if(steps) intsteps=*steps;
- //else intsteps=0;
- //if(m_debug) {
- //cout<<"CaloGeoGeometry::getDDE : check map"<<j<<" skip_range_check="<<skip_range_check<<endl;
- //}
- //float newdist;
- //const CaloGeoDetDescrElement* newDDE=m_cells_in_regions[sampling][j]->getDDE(eta,phi,&newdist,&intsteps);
- //if(m_debug) {
- //cout<<"CaloGeoGeometry::getDDE : map"<<j<<" dist="<<newdist<<" best dist="<<*distance<<" steps="<<intsteps<<endl;
- //}
- //if(newdist<*distance) {
- //bestDDE=newDDE;
- //*distance=newdist;
- //if(steps) beststeps=intsteps;
- //if(newdist<-0.1) break; //stop, we are well within the hit cell
- //}
- //}
- //if(bestDDE) break;
- //}
- }
- if(steps) *steps=beststeps;
- return bestDDE;
-}
-
-const CaloGeoDetDescrElement* CaloGeoGeometry::getFCalDDE(int sampling,float x,float y,float z){
- int isam = sampling - 20;
- int iphi,ieta;
- Long64_t mask1[]{0x34,0x34,0x35};
- Long64_t mask2[]{0x36,0x36,0x37};
- m_FCal_ChannelMap.getTileID(isam, x, y, ieta, iphi);
- //cout << ieta << ""
- Long64_t id = (ieta << 5) + 2*iphi;
- if(isam==2)id+= (8<<8);
-
- if(z>0) id+=(mask2[isam-1] << 12);
- else id+=(mask1[isam-1] << 12);
- //if(z<0) cout << "CaloGeoGeometry::getFCalDDE: Identifier: " << (id << 44) << " " << m_cells[id << 44]->identify() << endl;
- id = id << 44;
- Identifier identify((unsigned long long)id);
- return m_cells[identify];
- //return m_cells_in_sampling[sampling][id << 12];
-}
-
-
-bool CaloGeoGeometry::PostProcessGeometry()
-{
- for(int i=0;i<MAX_SAMPLING;++i) {
- post_process(i);
- for(unsigned int j=0;j<m_cells_in_regions[i].size();++j) {
- m_cells_in_regions[i][j]->post_process();
- }
- //if(i>=21) break;
- }
-
- InitRZmaps();
-
- /*
- cout<<"all : "<<m_cells.size()<<endl;
- for(int sampling=0;sampling<MAX_SAMPLING;++sampling) {
- cout<<"cells sampling "<<sampling<<" : "<<m_cells_in_sampling[sampling].size()<<" cells";
- cout<<", "<<m_cells_in_regions[sampling].size()<<" lookup map(s)"<<endl;
- for(unsigned int j=0;j<m_cells_in_regions[sampling].size();++j) {
- PrintMapInfo(sampling,j);
- //break;
- }
- //if(i>0) break;
- }
- */
-
- return true;
-}
-
-void CaloGeoGeometry::Validate()
-{
- int ntest=0;
- cout<<"start CaloGeoGeometry::Validate()"<<endl;
- for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
- const CaloGeoDetDescrElement* cell=ic->second;
- int sampling=cell->getSampling();
- if(sampling>=21) continue;
-
- if(m_debug_identify==cell->identify()) {
- cout<<"CaloGeoGeometry::Validate(), cell "<<ntest<<" id="<<cell->identify()<<endl;
- m_debug=true;
- }
-
- const int nrnd=100;
- for(int irnd=0;irnd<nrnd;++irnd) {
- int steps=0;
- float eta=cell->eta()+1.5*(gRandom->Rndm()-0.5)*cell->deta();
- float phi=cell->phi()+1.5*(gRandom->Rndm()-0.5)*cell->dphi();
- float distance;
- const CaloGeoDetDescrElement* foundcell=getDDE(sampling,eta,phi,&distance,&steps);
- if(m_debug && foundcell) {
- cout<<"CaloGeoGeometry::Validate(), irnd="<<irnd<<", foundcell id="<<foundcell->identify()<<", "<<steps<<" steps"<<endl;
- }
- if(cell==foundcell) {
- if(steps>3 && distance<-0.01) {
- cout<<"cell id="<<cell->identify()<<", sampling="<<sampling<<" found in "<<steps<<" steps, dist="<<distance<<" eta="<<eta<<" phi="<<phi<<endl;
- cout<<" eta="<<cell->eta()<<" eta_raw="<<cell->eta_raw()<<" deta="<<cell->deta()<<" ("<<(cell->eta_raw()-cell->eta())/cell->deta()<<") phi="<<cell->phi()<<" phi_raw="<<cell->phi_raw()<<" dphi="<<cell->dphi()<<" ("<<(cell->phi_raw()-cell->phi())/cell->dphi()<<")"<<endl;
- //if(steps>3 && distance>0.01) return;
- }
- } else {
- if( TMath::Abs( (eta-cell->eta())/cell->deta() )<0.45 && TMath::Abs( (phi-cell->phi())/cell->dphi() )<0.45 ) {
- cout<<"cell id="<<cell->identify()<<" not found! Found instead id=";
- if (foundcell) cout << foundcell->identify();
- cout <<" in "<<steps<<" steps, dist="<<distance<<" eta="<<eta<<" phi="<<phi<<endl;
- cout<<" input sampling="<<sampling<<" eta="<<cell->eta()<<" eta_raw="<<cell->eta_raw()<<" deta="<<cell->deta()<<" ("<<(cell->eta_raw()-cell->eta())/cell->deta()<<") phi="<<cell->phi()<<" phi_raw="<<cell->phi_raw()<<" dphi="<<cell->dphi()<<" ("<<(cell->phi_raw()-cell->phi())/cell->dphi()<<")"<<endl;
- cout<<" output sampling="<<(foundcell?(foundcell->getSampling()):-1)<<" eta="<<(foundcell?(foundcell->eta()):0)<<" eta_raw="<<(foundcell?(foundcell->eta_raw()):0)<<" deta="<<(foundcell?(foundcell->deta()):0)<<" ("<<(foundcell?((foundcell->eta_raw()-foundcell->eta())/foundcell->deta()):0)<<") phi="<<(foundcell?(foundcell->phi()):0)<<" phi_raw="<<(foundcell?(foundcell->phi_raw()):0)<<" dphi="<<(foundcell?(foundcell->dphi()):0)<<" ("<<(foundcell?((foundcell->phi_raw()-foundcell->phi())/cell->dphi()):0)<<")"<<endl;
- return;
- }
- if(!foundcell) {
- cout<<"nothing found close to cell id="<<cell->identify()<<" in "<<steps<<" steps, dist="<<distance<<" eta="<<eta<<" phi="<<phi<<endl;
- cout<<" input sampling="<<sampling<<" eta="<<cell->eta()<<" eta_raw="<<cell->eta_raw()<<" deta="<<cell->deta()<<" ("<<(cell->eta_raw()-cell->eta())/cell->deta()<<") phi="<<cell->phi()<<" phi_raw="<<cell->phi_raw()<<" dphi="<<cell->dphi()<<" ("<<(cell->phi_raw()-cell->phi())/cell->dphi()<<")"<<endl;
- return;
- }
- }
- //if(ntest>60000) break;
- }
- m_debug=false;
- if(ntest%25000==0) cout<<"Validate cell "<<ntest<<" with "<<nrnd<<" random hits"<<endl;
- ++ntest;
- }
- cout<<"end CaloGeoGeometry::Validate()"<<endl;
-}
-
-double CaloGeoGeometry::deta(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- double mineta=m_min_eta_sample[side][sample];
- double maxeta=m_max_eta_sample[side][sample];
-
- if(eta<mineta)
- {
- return fabs(eta-mineta);
- }
- else if(eta>maxeta)
- {
- return fabs(eta-maxeta);
- }
- else
- {
- double d1=fabs(eta-mineta);
- double d2=fabs(eta-maxeta);
- if(d1<d2) return -d1;
- else return -d2;
- }
-}
-
-
-void CaloGeoGeometry::minmaxeta(int sample,double eta,double& mineta,double& maxeta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- mineta=m_min_eta_sample[side][sample];
- maxeta=m_max_eta_sample[side][sample];
-}
-
-double CaloGeoGeometry::rmid(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- return m_rmid_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::zmid(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- return m_zmid_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::rzmid(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- if(isCaloBarrel(sample)) return m_rmid_map[side][sample].find_closest(eta)->second;
- else return m_zmid_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::rent(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- return m_rent_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::zent(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- return m_zent_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::rzent(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- if(isCaloBarrel(sample)) return m_rent_map[side][sample].find_closest(eta)->second;
- else return m_zent_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::rext(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- return m_rext_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::zext(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- return m_zext_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::rzext(int sample,double eta) const
-{
- int side=0;
- if(eta>0) side=1;
-
- if(isCaloBarrel(sample)) return m_rext_map[side][sample].find_closest(eta)->second;
- else return m_zext_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::rpos(int sample,double eta,int subpos) const
-{
-
- int side=0;
- if(eta>0) side=1;
-
- if(subpos==CaloSubPos::SUBPOS_ENT) return m_rent_map[side][sample].find_closest(eta)->second;
- if(subpos==CaloSubPos::SUBPOS_EXT) return m_rext_map[side][sample].find_closest(eta)->second;
-
- return m_rmid_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::zpos(int sample,double eta,int subpos) const
-{
- int side=0;
- if(eta>0) side=1;
-
- if(subpos==CaloSubPos::SUBPOS_ENT) return m_zent_map[side][sample].find_closest(eta)->second;
- if(subpos==CaloSubPos::SUBPOS_EXT) return m_zext_map[side][sample].find_closest(eta)->second;
- return m_zmid_map[side][sample].find_closest(eta)->second;
-}
-
-double CaloGeoGeometry::rzpos(int sample,double eta,int subpos) const
-{
- int side=0;
- if(eta>0) side=1;
-
- if(isCaloBarrel(sample)) {
- if(subpos==CaloSubPos::SUBPOS_ENT) return m_rent_map[side][sample].find_closest(eta)->second;
- if(subpos==CaloSubPos::SUBPOS_EXT) return m_rext_map[side][sample].find_closest(eta)->second;
- return m_rmid_map[side][sample].find_closest(eta)->second;
- } else {
- if(subpos==CaloSubPos::SUBPOS_ENT) return m_zent_map[side][sample].find_closest(eta)->second;
- if(subpos==CaloSubPos::SUBPOS_EXT) return m_zext_map[side][sample].find_closest(eta)->second;
- return m_zmid_map[side][sample].find_closest(eta)->second;
- }
-}
-
-std::string CaloGeoGeometry::SamplingName(int sample)
-{
- return CaloSampling::getSamplingName(sample);
-}
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx
index 6caff6fc20be..ca49c0552219 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometry.cxx
@@ -11,6 +11,7 @@
#include <TGraphErrors.h>
#include <TVector3.h>
#include <TLegend.h>
+#include <fstream>
#include "CaloDetDescr/CaloDetDescrElement.h"
//#include "ISF_FastCaloSimParametrization/CaloDetDescrElement.h"
@@ -26,420 +27,6 @@ const int CaloGeometry::MAX_SAMPLING = CaloCell_ID_FCS::MaxSample; //number of c
Identifier CaloGeometry::m_debug_identify;
bool CaloGeometry::m_debug=false;
-CaloGeometryLookup::CaloGeometryLookup(int ind):m_xy_grid_adjustment_factor(0.75),m_index(ind)
-{
- m_mineta=+10000;
- m_maxeta=-10000;
- m_minphi=+10000;
- m_maxphi=-10000;
-
- m_mineta_raw=+10000;
- m_maxeta_raw=-10000;
- m_minphi_raw=+10000;
- m_maxphi_raw=-10000;
-
- m_mindeta=10000;
- m_mindphi=10000;
-
- m_mineta_correction=+10000;
- m_maxeta_correction=-10000;
- m_minphi_correction=+10000;
- m_maxphi_correction=-10000;
-
- m_cell_grid_eta=0.;
- m_cell_grid_phi=0.;
- m_deta_double =0.;
- m_dphi_double =0.;
-}
-
-CaloGeometryLookup::~CaloGeometryLookup()
-{
-}
-
-bool CaloGeometryLookup::has_overlap(CaloGeometryLookup* ref)
-{
- if(m_cells.size()==0) return false;
- for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
- const CaloDetDescrElement* cell=ic->second;
- if(ref->IsCompatible(cell)) return true;
- }
- return false;
-}
-
-void CaloGeometryLookup::merge_into_ref(CaloGeometryLookup* ref)
-{
- for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
- const CaloDetDescrElement* cell=ic->second;
- ref->add(cell);
- }
-}
-
-
-bool CaloGeometryLookup::IsCompatible(const CaloDetDescrElement* cell)
-{
- if(m_cells.size()==0) return true;
- t_cellmap::iterator ic=m_cells.begin();
- const CaloDetDescrElement* refcell=ic->second;
- int sampling=refcell->getSampling();
- if(cell->getSampling()!=sampling) return false;
- if(cell->eta_raw()*refcell->eta_raw()<0) return false;
- if(sampling<21) { // keep only cells reasonable close in eta to each other
- if(TMath::Abs(cell->deta()-refcell->deta())>0.001) return false;
- if(TMath::Abs(cell->dphi()-refcell->dphi())>0.001) return false;
-
- if( cell->eta()<mineta()-2.1*cell->deta() ) return false;
- if( cell->eta()>maxeta()+2.1*cell->deta() ) return false;
-
- double delta_eta=(cell->eta_raw()-refcell->eta_raw())/cell->deta();
- double delta_phi=(cell->phi_raw()-refcell->phi_raw())/cell->dphi();
- if(TMath::Abs(delta_eta-TMath::Nint(delta_eta)) > 0.01) return false;
- if(TMath::Abs(delta_phi-TMath::Nint(delta_phi)) > 0.01) return false;
-
- /*
- cout<<"Compatible: s="<<cell->getSampling()<<"~"<<refcell->getSampling()<<"; ";
- cout<<"eta="<<cell->eta_raw()<<","<<refcell->eta_raw()<<"; ";
- cout<<"phi="<<cell->phi_raw()<<","<<refcell->phi_raw()<<"; ";
- cout<<"deta="<<cell->deta()<<"~"<<refcell->deta()<<"; ";
- cout<<"dphi="<<cell->dphi()<<"~"<<refcell->dphi()<<";";
- cout<<"mineta="<<mineta()<<", maxeta="<<maxeta()<<"; ";
- cout<<endl;
- */
- } else {
- //FCal is not sufficiently regular to use submaps with regular mapping
- return true;
- /*
- if(TMath::Abs(cell->dx()-refcell->dx())>0.001) return false;
- if(TMath::Abs(cell->dy()-refcell->dy())>0.001) return false;
- if( cell->x()<minx()-20*cell->dx() ) return false;
- if( cell->x()>maxx()+20*cell->dx() ) return false;
- if( cell->y()<miny()-20*cell->dy() ) return false;
- if( cell->y()>maxy()+20*cell->dy() ) return false;
-
- double delta_x=(cell->x_raw()-refcell->x_raw())/cell->dx();
- double delta_y=(cell->y_raw()-refcell->y_raw())/cell->dy();
- if(TMath::Abs(delta_x-TMath::Nint(delta_x)) > 0.01) return false;
- if(TMath::Abs(delta_y-TMath::Nint(delta_y)) > 0.01) return false;
- */
- }
-
- return true;
-}
-
-void CaloGeometryLookup::add(const CaloDetDescrElement* cell)
-{
- if(cell->getSampling()<21) {
- m_deta.add(cell->deta());
- m_dphi.add(cell->dphi());
- m_mindeta=TMath::Min(cell->deta(),m_mindeta);
- m_mindphi=TMath::Min(cell->dphi(),m_mindphi);
-
- m_eta_correction.add(cell->eta_raw()-cell->eta());
- m_phi_correction.add(cell->phi_raw()-cell->phi());
- m_mineta_correction=TMath::Min(cell->eta_raw()-cell->eta() , m_mineta_correction);
- m_maxeta_correction=TMath::Max(cell->eta_raw()-cell->eta() , m_maxeta_correction);
- m_minphi_correction=TMath::Min(cell->phi_raw()-cell->phi() , m_minphi_correction);
- m_maxphi_correction=TMath::Max(cell->phi_raw()-cell->phi() , m_maxphi_correction);
-
- m_mineta=TMath::Min(cell->eta()-cell->deta()/2,m_mineta);
- m_maxeta=TMath::Max(cell->eta()+cell->deta()/2,m_maxeta);
- m_minphi=TMath::Min(cell->phi()-cell->dphi()/2,m_minphi);
- m_maxphi=TMath::Max(cell->phi()+cell->dphi()/2,m_maxphi);
-
- m_mineta_raw=TMath::Min(cell->eta_raw()-cell->deta()/2,m_mineta_raw);
- m_maxeta_raw=TMath::Max(cell->eta_raw()+cell->deta()/2,m_maxeta_raw);
- m_minphi_raw=TMath::Min(cell->phi_raw()-cell->dphi()/2,m_minphi_raw);
- m_maxphi_raw=TMath::Max(cell->phi_raw()+cell->dphi()/2,m_maxphi_raw);
- } else {
- m_deta.add(cell->dx());
- m_dphi.add(cell->dy());
- m_mindeta=TMath::Min(cell->dx(),m_mindeta);
- m_mindphi=TMath::Min(cell->dy(),m_mindphi);
-
- m_eta_correction.add(cell->x_raw()-cell->x());
- m_phi_correction.add(cell->y_raw()-cell->y());
- m_mineta_correction=TMath::Min(cell->x_raw()-cell->x() , m_mineta_correction);
- m_maxeta_correction=TMath::Max(cell->x_raw()-cell->x() , m_maxeta_correction);
- m_minphi_correction=TMath::Min(cell->y_raw()-cell->y() , m_minphi_correction);
- m_maxphi_correction=TMath::Max(cell->y_raw()-cell->y() , m_maxphi_correction);
-
- m_mineta=TMath::Min(cell->x()-cell->dx()/2,m_mineta);
- m_maxeta=TMath::Max(cell->x()+cell->dx()/2,m_maxeta);
- m_minphi=TMath::Min(cell->y()-cell->dy()/2,m_minphi);
- m_maxphi=TMath::Max(cell->y()+cell->dy()/2,m_maxphi);
-
- m_mineta_raw=TMath::Min(cell->x_raw()-cell->dx()/2,m_mineta_raw);
- m_maxeta_raw=TMath::Max(cell->x_raw()+cell->dx()/2,m_maxeta_raw);
- m_minphi_raw=TMath::Min(cell->y_raw()-cell->dy()/2,m_minphi_raw);
- m_maxphi_raw=TMath::Max(cell->y_raw()+cell->dy()/2,m_maxphi_raw);
- }
- m_cells[cell->identify()]=cell;
-}
-
-bool CaloGeometryLookup::index_range_adjust(int& ieta,int& iphi)
-{
- while(iphi<0) {iphi+=m_cell_grid_phi;};
- while(iphi>=m_cell_grid_phi) {iphi-=m_cell_grid_phi;};
- if(ieta<0) {
- ieta=0;
- return false;
- }
- if(ieta>=m_cell_grid_eta) {
- ieta=m_cell_grid_eta-1;
- return false;
- }
- return true;
-}
-
-void CaloGeometryLookup::post_process()
-{
- if(size()==0) return;
- t_cellmap::iterator ic=m_cells.begin();
- const CaloDetDescrElement* refcell=ic->second;
- int sampling=refcell->getSampling();
- if(sampling<21) {
- double rneta=neta_double()-neta();
- double rnphi=nphi_double()-nphi();
- if(TMath::Abs(rneta)>0.05 || TMath::Abs(rnphi)>0.05) {
- cout<<"ERROR: can't map cells into a regular grid for sampling "<<sampling<<", map index "<<index()<<endl;
- return;
- }
- m_cell_grid_eta=neta();
- m_cell_grid_phi=nphi();
- m_deta_double=deta();
- m_dphi_double=dphi();
- } else {
- //double rnx=nx_double()-nx();
- //double rny=ny_double()-ny();
- //if(TMath::Abs(rnx)>0.05 || TMath::Abs(rny)>0.05) {
- // cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": mapping cells into a regular grid, although cells don't fit well"<<endl;
- //}
-
- m_cell_grid_eta=TMath::Nint(TMath::Ceil(nx_double()/m_xy_grid_adjustment_factor));
- m_cell_grid_phi=TMath::Nint(TMath::Ceil(ny_double()/m_xy_grid_adjustment_factor));
- m_deta_double=mindx()*m_xy_grid_adjustment_factor;
- m_dphi_double=mindy()*m_xy_grid_adjustment_factor;
- }
-
- m_cell_grid.resize(m_cell_grid_eta);
- for(int ieta=0;ieta<m_cell_grid_eta;++ieta) {
- m_cell_grid[ieta].resize(m_cell_grid_phi,0);
- }
-
- for(ic=m_cells.begin();ic!=m_cells.end();++ic) {
- refcell=ic->second;
- int ieta,iphi;
- if(sampling<21) {
- ieta=raw_eta_position_to_index(refcell->eta_raw());
- iphi=raw_phi_position_to_index(refcell->phi_raw());
- } else {
- ieta=raw_eta_position_to_index(refcell->x_raw());
- iphi=raw_phi_position_to_index(refcell->y_raw());
- }
- if(index_range_adjust(ieta,iphi)) {
- if(m_cell_grid[ieta][iphi]) {
- cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": Already cell found at pos ("<<ieta<<","<<iphi<<"): id="
- <<m_cell_grid[ieta][iphi]->identify()<<"->"<<refcell->identify()<<endl;
- cout<<" x="<<m_cells[m_cell_grid[ieta][iphi]->identify()]->x_raw()<<" -> "<<refcell->x_raw();
- cout<<" , y="<<m_cells[m_cell_grid[ieta][iphi]->identify()]->y_raw()<<" -> "<<refcell->y_raw();
- cout<<" mindx="<<m_deta_double<<" mindy="<<m_dphi_double<<endl;
- cout<<endl;
- }
- m_cell_grid[ieta][iphi]=refcell;
- } else {
- cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": Cell at pos ("<<ieta<<","<<iphi<<"): id="
- <<refcell->identify()<<" outside eta range"<<endl;
- }
- }
-
- int ncells=0;
- int nempty=0;
- for(int ieta=0;ieta<m_cell_grid_eta;++ieta) {
- for(int iphi=0;iphi<m_cell_grid_phi;++iphi) {
- if(!m_cell_grid[ieta][iphi]) {
- ++nempty;
- //cout<<"Sampling "<<sampling<<"_"<<index()<<": No cell at pos ("<<ieta<<","<<iphi<<")"<<endl;
- } else {
- ++ncells;
- }
- }
- }
- // cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": "<<ncells<<"/"<<size()<<" cells filled, "<<nempty<<" empty grid positions deta="<<m_deta_double<<" dphi="<<m_dphi_double<<endl;
-}
-
-float CaloGeometryLookup::calculate_distance_eta_phi(const CaloDetDescrElement* DDE,float eta,float phi,float& dist_eta0,float& dist_phi0)
-{
- dist_eta0=(eta - DDE->eta())/m_deta_double;
- dist_phi0=(TVector2::Phi_mpi_pi(phi - DDE->phi()))/m_dphi_double;
- float abs_dist_eta0=TMath::Abs(dist_eta0);
- float abs_dist_phi0=TMath::Abs(dist_phi0);
- return TMath::Max(abs_dist_eta0,abs_dist_phi0)-0.5;
-}
-
-const CaloDetDescrElement* CaloGeometryLookup::getDDE(float eta,float phi,float* distance,int* steps)
-{
- float dist;
- const CaloDetDescrElement* bestDDE=0;
- if(!distance) distance=&dist;
- *distance=+10000000;
- int intsteps=0;
- if(!steps) steps=&intsteps;
-
- float best_eta_corr=m_eta_correction;
- float best_phi_corr=m_phi_correction;
-
- float raw_eta=eta+best_eta_corr;
- float raw_phi=phi+best_phi_corr;
-
- int ieta=raw_eta_position_to_index(raw_eta);
- int iphi=raw_phi_position_to_index(raw_phi);
- index_range_adjust(ieta,iphi);
-
- const CaloDetDescrElement* newDDE=m_cell_grid[ieta][iphi];
- float bestdist=+10000000;
- ++(*steps);
- int nsearch=0;
- while(newDDE && nsearch<3) {
- float dist_eta0,dist_phi0;
- *distance=calculate_distance_eta_phi(newDDE,eta,phi,dist_eta0,dist_phi0);
- bestDDE=newDDE;
- bestdist=*distance;
-
- if(CaloGeometry::m_debug || CaloGeometry::m_debug_identify==newDDE->identify()) {
- cout<<"CaloGeometryLookup::getDDE, search iteration "<<nsearch<<endl;
- cout<<"cell id="<<newDDE->identify()<<" steps "<<*steps<<" steps, eta="<<eta<<" phi="<<phi<<" dist="<<*distance<<" cell_grid_eta="<<cell_grid_eta()<<" cell_grid_phi="<<cell_grid_phi()<<" m_deta_double="<<m_deta_double<<" m_dphi_double="<<m_dphi_double<<" 1st_eta_corr="<<best_eta_corr<<" 1st_phi_corr="<<best_phi_corr<<endl;
- cout<<" input sampling="<<newDDE->getSampling()<<" eta="<<newDDE->eta()<<" eta_raw="<<newDDE->eta_raw()<<" deta="<<newDDE->deta()<<" ("<<(newDDE->eta_raw()-newDDE->eta())/newDDE->deta()<<") phi="<<newDDE->phi()<<" phi_raw="<<newDDE->phi_raw()<<" dphi="<<newDDE->dphi()<<" ("<<(newDDE->phi_raw()-newDDE->phi())/newDDE->dphi()<<")"<<endl;
- cout<<" ieta="<<ieta<<" iphi="<<iphi<<endl;
- cout<<" dist_eta0="<<dist_eta0<<" ,dist_phi0="<<dist_phi0<<endl;
- }
-
- if(*distance<0) return newDDE;
- //correct ieta and iphi by the observed difference to the hit cell
- ieta+=TMath::Nint(dist_eta0);
- iphi+=TMath::Nint(dist_phi0);
- index_range_adjust(ieta,iphi);
- const CaloDetDescrElement* oldDDE=newDDE;
- newDDE=m_cell_grid[ieta][iphi];
- ++(*steps);
- ++nsearch;
- if(oldDDE==newDDE) break;
- }
- if(CaloGeometry::m_debug && !newDDE) {
- cout<<"CaloGeometryLookup::getDDE, direct search ieta="<<ieta<<" iphi="<<iphi<<" is empty"<<endl;
- }
- float minieta=ieta+TMath::Nint(TMath::Floor(m_mineta_correction/cell_grid_eta()));
- float maxieta=ieta+TMath::Nint(TMath::Ceil(m_maxeta_correction/cell_grid_eta()));
- float miniphi=iphi+TMath::Nint(TMath::Floor(m_minphi_correction/cell_grid_phi()));
- float maxiphi=iphi+TMath::Nint(TMath::Ceil(m_maxphi_correction/cell_grid_phi()));
- if(minieta<0) minieta=0;
- if(maxieta>=m_cell_grid_eta) maxieta=m_cell_grid_eta-1;
- for(int iieta=minieta;iieta<=maxieta;++iieta) {
- for(int iiphi=miniphi;iiphi<=maxiphi;++iiphi) {
- ieta=iieta;
- iphi=iiphi;
- index_range_adjust(ieta,iphi);
- newDDE=m_cell_grid[ieta][iphi];
- ++(*steps);
- if(newDDE) {
- float dist_eta0,dist_phi0;
- *distance=calculate_distance_eta_phi(newDDE,eta,phi,dist_eta0,dist_phi0);
-
- if(CaloGeometry::m_debug || CaloGeometry::m_debug_identify==newDDE->identify()) {
- cout<<"CaloGeometryLookup::getDDE, windows search! minieta="<<m_mineta_correction/cell_grid_eta()<<" maxieta="<<m_maxeta_correction/cell_grid_eta()<<" miniphi="<<m_minphi_correction/cell_grid_phi()<<" maxiphi="<<m_maxphi_correction/cell_grid_phi()<<endl;
- cout<<"cell id="<<newDDE->identify()<<" steps "<<*steps<<" steps, eta="<<eta<<" phi="<<phi<<" dist="<<*distance<<endl;
- cout<<" input sampling="<<newDDE->getSampling()<<" eta="<<newDDE->eta()<<" eta_raw="<<newDDE->eta_raw()<<" deta="<<newDDE->deta()<<" ("<<(newDDE->eta_raw()-newDDE->eta())/newDDE->deta()<<") phi="<<newDDE->phi()<<" phi_raw="<<newDDE->phi_raw()<<" dphi="<<newDDE->dphi()<<" ("<<(newDDE->phi_raw()-newDDE->phi())/newDDE->dphi()<<")"<<endl;
- cout<<" ieta="<<ieta<<" iphi="<<iphi<<endl;
- cout<<" dist_eta0="<<dist_eta0<<" ,dist_phi0="<<dist_phi0<<endl;
- }
-
- if(*distance<0) return newDDE;
- if(*distance<bestdist) {
- bestDDE=newDDE;
- bestdist=*distance;
- }
- } else if(CaloGeometry::m_debug) {
- cout<<"CaloGeometryLookup::getDDE, windows search ieta="<<ieta<<" iphi="<<iphi<<" is empty"<<endl;
- }
- }
- }
- *distance=bestdist;
- return bestDDE;
-}
-
-
-
-/*
-void CaloGeometryLookup::CalculateTransformation()
-{
- gROOT->cd();
-
- TTree* tmap = new TTree( "mapping" , "mapping" );
-
- float eta,phi,Deta_raw,Dphi_raw;
- tmap->Branch("eta", &eta,"eta/F");
- tmap->Branch("phi", &phi,"phi/F");
- tmap->Branch("Deta_raw", &Deta_raw,"Deta_raw/F");
- tmap->Branch("Dphi_raw", &Dphi_raw,"Dphi_raw/F");
-
- if(m_cells.size()==0) return;
-
- int sampling=0;
- for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
- CaloGeoDetDescrElement* refcell=ic->second;
- sampling=refcell->getSampling();
- if(sampling<21) {
- eta=refcell->eta();
- phi=refcell->phi();
- Deta_raw=refcell->eta_raw()-eta;
- Dphi_raw=refcell->phi_raw()-phi;
- } else {
- eta=refcell->x();
- phi=refcell->y();
- Deta_raw=refcell->x_raw()-eta;
- Dphi_raw=refcell->y_raw()-phi;
- }
- tmap->Fill();
- tmap->Fill(); //Fill twice to have all events and training and test tree
- }
- tmap->Print();
-
- TString outfileName( Form("Mapping%d_%d.root",sampling,m_index) );
- TFile* outputFile = new TFile( outfileName, "RECREATE" );
- //TFile* outputFile = new TMemFile( outfileName, "RECREATE" );
-
- TMVA::Factory *factory = new TMVA::Factory( Form("Mapping%d_%d.root",sampling,m_index) , outputFile, "!V:!Silent:Color:DrawProgressBar" );
-
- factory->AddVariable( "eta", "calo eta", "1", 'F' );
- factory->AddVariable( "phi", "calo phi", "1", 'F' );
- factory->AddTarget( "Deta_raw" );
- factory->AddTarget( "Dphi_raw" );
-
- factory->AddRegressionTree( tmap );
-
- //factory->PrepareTrainingAndTestTree( "" , Form("nTrain_Regression=%d:nTest_Regression=%d:SplitMode=Random:NormMode=NumEvents:!V",(int)m_cells.size(),(int)m_cells.size()) );
- factory->PrepareTrainingAndTestTree( "" , "nTrain_Regression=0:nTest_Regression=0:SplitMode=Alternate:NormMode=NumEvents:!V" );
-
- factory->BookMethod( TMVA::Types::kMLP, "MLP", "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+5:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator" );
-
- cout<<"=== Start trainging ==="<<endl;
- // Train MVAs using the set of training events
- factory->TrainAllMethods();
-
- cout<<"=== Start testing ==="<<endl;
- // ---- Evaluate all MVAs using the set of test events
- factory->TestAllMethods();
-
- cout<<"=== Start evaluation ==="<<endl;
- // ----- Evaluate and compare performance of all configured MVAs
- factory->EvaluateAllMethods();
-
- outputFile->Close();
-
- delete factory;
-
- delete tmap;
-}
-*/
-
CaloGeometry::CaloGeometry() : m_cells_in_sampling(MAX_SAMPLING),m_cells_in_sampling_for_phi0(MAX_SAMPLING),m_cells_in_regions(MAX_SAMPLING),m_isCaloBarrel(MAX_SAMPLING),m_dographs(false),m_FCal_ChannelMap(0)
{
//TMVA::Tools::Instance();
@@ -1214,3 +801,81 @@ std::string CaloGeometry::SamplingName(int sample)
return CaloSampling::getSamplingName(sample);
}
+void CaloGeometry::LoadFCalGeometryFromFiles(TString filename1,TString filename2,TString filename3){
+
+ vector<ifstream*> electrodes(3);
+
+ electrodes[0]=new ifstream(filename1);
+ electrodes[1]=new ifstream(filename2);
+ electrodes[2]=new ifstream(filename3);
+
+
+ int thisTubeId;
+ int thisTubeI;
+ int thisTubeJ;
+ //int thisTubeID;
+ //int thisTubeMod;
+ double thisTubeX;
+ double thisTubeY;
+ TString tubeName;
+
+ //int second_column;
+ string seventh_column;
+ string eight_column;
+ int ninth_column;
+
+
+
+
+
+ int i;
+ for(int imodule=1;imodule<=3;imodule++){
+
+ i=0;
+ //while(i<50){
+ while(1){
+
+ //cout << electrodes[imodule-1]->eof() << endl;
+ (*electrodes[imodule-1]) >> tubeName;
+ if(electrodes[imodule-1]->eof())break;
+ (*electrodes[imodule-1]) >> thisTubeId; // ?????
+ (*electrodes[imodule-1]) >> thisTubeI;
+ (*electrodes[imodule-1]) >> thisTubeJ;
+ (*electrodes[imodule-1]) >> thisTubeX;
+ (*electrodes[imodule-1]) >> thisTubeY;
+ (*electrodes[imodule-1]) >> seventh_column;
+ (*electrodes[imodule-1]) >> eight_column;
+ (*electrodes[imodule-1]) >> ninth_column;
+
+ tubeName.ReplaceAll("'","");
+ string tubeNamestring=tubeName.Data();
+
+ std::istringstream tileStream1(std::string(tubeNamestring,1,1));
+ std::istringstream tileStream2(std::string(tubeNamestring,3,2));
+ std::istringstream tileStream3(std::string(tubeNamestring,6,3));
+ int a1=0,a2=0,a3=0;
+ if (tileStream1) tileStream1 >> a1;
+ if (tileStream2) tileStream2 >> a2;
+ if (tileStream3) tileStream3 >> a3;
+
+ //unsigned int tileName= (a3 << 16) + a2;
+ stringstream s;
+
+
+ m_FCal_ChannelMap.add_tube(tubeNamestring, imodule, thisTubeId, thisTubeI,thisTubeJ, thisTubeX, thisTubeY,seventh_column);
+
+
+
+ //cout << "FCal electrodes: " << tubeName << " " << second_column << " " << thisTubeI << " " << thisTubeJ << " " << thisTubeX << " " << thisTubeY << " " << seventh_column << " " << eight_column << " " << ninth_column << endl;
+ //cout << tileStream1.str() << " " << tileStream2.str() << " " << tileStream3.str() << endl;
+ //cout << a1 << " " << a2 << " " << a3 << " " << tileName << endl;
+ i++;
+ }
+ }
+
+
+ m_FCal_ChannelMap.finish(); // Creates maps
+
+
+}
+
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.cxx
index a309a2b7010a..ffc2ad9fae24 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.cxx
@@ -2,7 +2,7 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
-#include "CaloGeometryFromCaloDDM.h"
+#include "ISF_FastCaloSimParametrization/CaloGeometryFromCaloDDM.h"
#include "CaloDetDescr/CaloDetDescrElement.h"
//#include "ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h"
#include "CaloDetDescr/CaloDetDescrManager.h"
@@ -17,18 +17,18 @@ CaloGeometryFromCaloDDM::~CaloGeometryFromCaloDDM()
{
}
-bool CaloGeometryFromCaloDDM::LoadGeometryFromCaloDDM(const CaloDetDescrManager* /*calo_dd_man*/)
+bool CaloGeometryFromCaloDDM::LoadGeometryFromCaloDDM(const CaloDetDescrManager* calo_dd_man)
{
- //int jentry=0;
- //for(CaloDetDescrManager::calo_element_const_iterator calo_iter=calo_dd_man->element_begin();calo_iter<calo_dd_man->element_end();++calo_iter) {
- //const CaloGeoDetDescrElement* pcell=*calo_iter;
- //addcell(pcell);
-
- //if(jentry%25000==0) {
- //cout<<jentry<<" : "<<pcell->getSampling()<<", "<<pcell->identify()<<endl;
- //}
- //++jentry;
- //}
+ int jentry=0;
+ for(CaloDetDescrManager::calo_element_const_iterator calo_iter=calo_dd_man->element_begin();calo_iter<calo_dd_man->element_end();++calo_iter) {
+ const CaloDetDescrElement* pcell=*calo_iter;
+ addcell(pcell);
+
+ if(jentry%25000==0) {
+ cout<<jentry<<" : "<<pcell->getSampling()<<", "<<pcell->identify()<<endl;
+ }
+ ++jentry;
+ }
return PostProcessGeometry();
}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.h
deleted file mode 100644
index 169eff26fae2..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromCaloDDM.h
+++ /dev/null
@@ -1,21 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef CaloGeometryFromCaloDDM_h
-#define CaloGeometryFromCaloDDM_h
-
-#include "ISF_FastCaloSimParametrization/CaloGeometry.h"
-
-class CaloDetDescrManager;
-
-class CaloGeometryFromCaloDDM:public CaloGeometry {
-public :
- CaloGeometryFromCaloDDM();
- virtual ~CaloGeometryFromCaloDDM();
-
- virtual bool LoadGeometryFromCaloDDM(const CaloDetDescrManager* calo_dd_man);
-};
-
-#endif
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromFile.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromFile.cxx
deleted file mode 100644
index 76c98c6f3f07..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryFromFile.cxx
+++ /dev/null
@@ -1,373 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#include "ISF_FastCaloSimParametrization/CaloGeometryFromFile.h"
-#include <TTree.h>
-#include <TFile.h>
-#include "ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h"
-//#include "CaloDetDescr/CaloGeoDetDescrElement.h"
-#include <fstream>
-#include <sstream>
-#include <TGraph.h>
-#include "TH2D.h"
-
-#include <fstream>
-#include <map>
-#include <string>
-#include <sstream>
-
-using namespace std;
-
-map<unsigned long long, unsigned long long> g_cellId_vs_cellHashId_map;
-
-CaloGeometryFromFile::CaloGeometryFromFile() : CaloGeoGeometry()
-{
- ifstream textfile("/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/cellId_vs_cellHashId_map.txt");
- unsigned long long id, hash_id;
- cout << "Loading cellId_vs_cellHashId_map" << endl;
- int i=0;
- string line;
- stringstream s;
- while(1){
- //getline(textfile,line);
- s.str(line);
- textfile >> id;
- if(textfile.eof())break;
- textfile >> hash_id;
- g_cellId_vs_cellHashId_map[id]=hash_id;
- if(i%10000==0)cout << "Line: " << i << " line " << line << " id " << hex << id << " hash_id " << dec << hash_id << endl;
- i++;
- }
- cout << "Done." << endl;
-
-}
-
-CaloGeometryFromFile::~CaloGeometryFromFile()
-{
-}
-
-bool CaloGeometryFromFile::LoadGeometryFromFile(TString filename,TString treename)
-{
- TTree *tree;
- TFile *f = TFile::Open(filename);
- if(!f) return false;
- f->GetObject(treename,tree);
- if(!tree) return false;
-
- TTree* fChain = tree;
-
- CaloGeoDetDescrElement cell;
-
- // List of branches
- TBranch *b_identifier; //!
- TBranch *b_calosample; //!
- TBranch *b_eta; //!
- TBranch *b_phi; //!
- TBranch *b_r; //!
- TBranch *b_eta_raw; //!
- TBranch *b_phi_raw; //!
- TBranch *b_r_raw; //!
- TBranch *b_x; //!
- TBranch *b_y; //!
- TBranch *b_z; //!
- TBranch *b_x_raw; //!
- TBranch *b_y_raw; //!
- TBranch *b_z_raw; //!
- TBranch *b_deta; //!
- TBranch *b_dphi; //!
- TBranch *b_dr; //!
- TBranch *b_dx; //!
- TBranch *b_dy; //!
- TBranch *b_dz; //!
-
- fChain->SetMakeClass(1);
- fChain->SetBranchAddress("identifier", &cell.m_identify, &b_identifier);
- fChain->SetBranchAddress("calosample", &cell.m_calosample, &b_calosample);
- fChain->SetBranchAddress("eta", &cell.m_eta, &b_eta);
- fChain->SetBranchAddress("phi", &cell.m_phi, &b_phi);
- fChain->SetBranchAddress("r", &cell.m_r, &b_r);
- fChain->SetBranchAddress("eta_raw", &cell.m_eta_raw, &b_eta_raw);
- fChain->SetBranchAddress("phi_raw", &cell.m_phi_raw, &b_phi_raw);
- fChain->SetBranchAddress("r_raw", &cell.m_r_raw, &b_r_raw);
- fChain->SetBranchAddress("x", &cell.m_x, &b_x);
- fChain->SetBranchAddress("y", &cell.m_y, &b_y);
- fChain->SetBranchAddress("z", &cell.m_z, &b_z);
- fChain->SetBranchAddress("x_raw", &cell.m_x_raw, &b_x_raw);
- fChain->SetBranchAddress("y_raw", &cell.m_y_raw, &b_y_raw);
- fChain->SetBranchAddress("z_raw", &cell.m_z_raw, &b_z_raw);
- fChain->SetBranchAddress("deta", &cell.m_deta, &b_deta);
- fChain->SetBranchAddress("dphi", &cell.m_dphi, &b_dphi);
- fChain->SetBranchAddress("dr", &cell.m_dr, &b_dr);
- fChain->SetBranchAddress("dx", &cell.m_dx, &b_dx);
- fChain->SetBranchAddress("dy", &cell.m_dy, &b_dy);
- fChain->SetBranchAddress("dz", &cell.m_dz, &b_dz);
-
- Long64_t nentries = fChain->GetEntriesFast();
- for (Long64_t jentry=0; jentry<nentries;jentry++) {
- Long64_t ientry = fChain->LoadTree(jentry);
- if (ientry < 0) break;
- fChain->GetEntry(jentry);
-
- if (g_cellId_vs_cellHashId_map.find(cell.m_identify)!=g_cellId_vs_cellHashId_map.end()) cell.m_hash_id=g_cellId_vs_cellHashId_map[cell.m_identify];
- else cout << endl << "ERROR: Cell id not found in the cellId_vs_cellHashId_map!!!" << endl << endl;
-
- const CaloGeoDetDescrElement* pcell=new CaloGeoDetDescrElement(cell);
- int sampling = pcell->getSampling();
-
-
- if(sampling >20) continue;
- this->addcell(pcell);
-
- if(jentry%25000==0) {
- //if(jentry==nentries-1) {
- cout << "Checking loading cells from file" << endl;
- cout<<jentry<<" : "<<pcell->getSampling()<<", "<<pcell->identify()<<endl;
-
-
-
-
- //if(jentry>5) break;
- }
- }
- //cout<<"all : "<<m_cells.size()<<endl;
- //unsigned long long max(0);
- //unsigned long long min_id=m_cells_in_sampling[0].begin()->first;
- //for(int i=0;i<21;++i) {
- ////cout<<" cells sampling "<<i<<" : "<<m_cells_in_sampling[i].size()<<" cells";
- ////cout<<", "<<m_cells_in_regions[i].size()<<" lookup map(s)"<<endl;
- //for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
- ////cout << it->second->getSampling() << " " << it->first << endl;
- //if(min_id/10 >= it->first){
- ////cout << "Warning: Identifiers are not in increasing order!!!!" << endl;
- ////cout << min_id << " " << it->first << endl;
-
- //}
- //if(min_id > it->first)min_id = it->first;
- //}
- //}
- //cout << "Min id for samplings < 21: " << min_id << endl;
- delete f;
- //return true;
- bool ok = PostProcessGeometry();
-
- cout << "Result of PostProcessGeometry(): " << ok << endl;
- const CaloGeoDetDescrElement* mcell=0;
- unsigned long long cellid64(3179554531063103488);
- Identifier cellid(cellid64);
- mcell=this->getDDE(cellid); //This is working also for the FCal
-
- std::cout << "\n \n";
- std::cout << "Testing whether CaloGeoGeometry is loaded properly" << std::endl;
- if(!mcell)std::cout << "Cell is not found" << std::endl;
- std::cout << "Identifier " << mcell->identify() <<" sampling " << mcell->getSampling() << " eta: " << mcell->eta() << " phi: " << mcell->phi() << std::endl<< std::endl;
-
- const CaloGeoDetDescrElement* mcell2 = this->getDDE(mcell->getSampling(),mcell->eta(),mcell->phi());
- std::cout << "Identifier " << mcell2->identify() <<" sampling " << mcell2->getSampling() << " eta: " << mcell2->eta() << " phi: " << mcell2->phi() << std::endl<< std::endl;
-
- return ok;
-}
-
-
-bool CaloGeometryFromFile::LoadFCalGeometryFromFiles(TString filename1,TString filename2,TString filename3){
-
- vector<ifstream*> electrodes(3);
-
- electrodes[0]=new ifstream(filename1);
- electrodes[1]=new ifstream(filename2);
- electrodes[2]=new ifstream(filename3);
-
-
- int thisTubeId;
- int thisTubeI;
- int thisTubeJ;
- //int thisTubeID;
- //int thisTubeMod;
- double thisTubeX;
- double thisTubeY;
- TString tubeName;
-
- //int second_column;
- string seventh_column;
- string eight_column;
- int ninth_column;
-
-
-
-
-
- int i;
- for(int imodule=1;imodule<=3;imodule++){
-
- i=0;
- //while(i<50){
- while(1){
-
- //cout << electrodes[imodule-1]->eof() << endl;
- (*electrodes[imodule-1]) >> tubeName;
- if(electrodes[imodule-1]->eof())break;
- (*electrodes[imodule-1]) >> thisTubeId; // ?????
- (*electrodes[imodule-1]) >> thisTubeI;
- (*electrodes[imodule-1]) >> thisTubeJ;
- (*electrodes[imodule-1]) >> thisTubeX;
- (*electrodes[imodule-1]) >> thisTubeY;
- (*electrodes[imodule-1]) >> seventh_column;
- (*electrodes[imodule-1]) >> eight_column;
- (*electrodes[imodule-1]) >> ninth_column;
-
- tubeName.ReplaceAll("'","");
- string tubeNamestring=tubeName.Data();
-
- std::istringstream tileStream1(std::string(tubeNamestring,1,1));
- std::istringstream tileStream2(std::string(tubeNamestring,3,2));
- std::istringstream tileStream3(std::string(tubeNamestring,6,3));
- int a1=0,a2=0,a3=0;
- if (tileStream1) tileStream1 >> a1;
- if (tileStream2) tileStream2 >> a2;
- if (tileStream3) tileStream3 >> a3;
-
- //unsigned int tileName= (a3 << 16) + a2;
- stringstream s;
-
-
- m_FCal_ChannelMap.add_tube(tubeNamestring, imodule, thisTubeId, thisTubeI,thisTubeJ, thisTubeX, thisTubeY,seventh_column);
-
-
-
- //cout << "FCal electrodes: " << tubeName << " " << second_column << " " << thisTubeI << " " << thisTubeJ << " " << thisTubeX << " " << thisTubeY << " " << seventh_column << " " << eight_column << " " << ninth_column << endl;
- //cout << tileStream1.str() << " " << tileStream2.str() << " " << tileStream3.str() << endl;
- //cout << a1 << " " << a2 << " " << a3 << " " << tileName << endl;
- i++;
- }
- }
-
-
- m_FCal_ChannelMap.finish(); // Creates maps
-
- unsigned long long phi_index,eta_index;
- float x,y,dx,dy;
- long id;
- //auto it =m_cells_in_sampling[20].rbegin();
- //it--;
- //unsigned long long identify=it->first;
- //for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
- //
- //}
- long mask1[]{0x34,0x34,0x35};
- long mask2[]{0x36,0x36,0x37};
-
- for(int imap=1;imap<=3;imap++){
- for(auto it=m_FCal_ChannelMap.begin(imap);it!=m_FCal_ChannelMap.end(imap);it++){
- phi_index = it->first & 0xffff;
- eta_index = it->first >> 16;
- CaloGeoDetDescrElement *DDE =new CaloGeoDetDescrElement; // side C
- CaloGeoDetDescrElement *DDE2 =new CaloGeoDetDescrElement; // side A
- x=it->second.x();
- y=it->second.y();
- m_FCal_ChannelMap.tileSize(imap, eta_index, phi_index,dx,dy);
- //if(imap==2) eta_index+=100;
- //if(imap==3) eta_index+=200;
- //cout << imap << " " << eta_index << " " << phi_index << " " << it->first << " " << (eta_index << 16) + phi_index << endl;
-
- DDE->m_calosample=imap+20;
- DDE->m_x=x;
- DDE->m_y=y;
- DDE->m_x_raw=x;
- DDE->m_y_raw=y;
- DDE->m_dx=dx;
- DDE->m_dy=dy;
-
- DDE2->m_calosample=imap+20;
- DDE2->m_x=x;
- DDE2->m_y=y;
- DDE2->m_x_raw=x;
- DDE2->m_y_raw=y;
- DDE2->m_dx=dx;
- DDE2->m_dy=dy;
-
- id=(mask1[imap-1]<<12) + (eta_index << 5) +2*phi_index;
- if(imap==2) id+= (8<<8);
- DDE->m_identify=(id << 44);
- //DDE->m_identify=(id << 12);
- //cout << DDE->identify() << endl;
- if (g_cellId_vs_cellHashId_map.find(DDE->m_identify)!=g_cellId_vs_cellHashId_map.end()) DDE->m_hash_id=g_cellId_vs_cellHashId_map[DDE->m_identify];
- else cout << endl << "ERROR: Cell id not found in the cellId_vs_cellHashId_map!!!" << endl << endl;
-
- addcell(DDE);
- id=(mask2[imap-1]<<12) + (eta_index << 5) +2*phi_index;
- if(imap==2) id+= (8<<8);
- DDE2->m_identify=(id << 44);
- //DDE->m_identify=(id << 12);
- if (g_cellId_vs_cellHashId_map.find(DDE2->m_identify)!=g_cellId_vs_cellHashId_map.end()) DDE2->m_hash_id=g_cellId_vs_cellHashId_map[DDE2->m_identify];
- else cout << endl << "ERROR: Cell id not found in the cellId_vs_cellHashId_map!!!" << endl << endl;
- addcell(DDE2);
- delete DDE;
- delete DDE2;
- }
-
- }
- //auto it =m_cells_in_sampling[0].begin();
- //it--;
- //unsigned long long identify=it->first;
- ////cout << "min identifier from sampling < 21: " << identify << endl;
-
- //for(int i=0;i<21;++i) {
- ////for(int i=21;i<MAX_SAMPLING;++i) {
- //cout<<" cells sampling "<<i<<" : "<<m_cells_in_sampling[i].size()<<" cells";
- //cout<<", "<<m_cells_in_regions[i].size()<<" lookup map(s)"<<endl;
- //int k=0;
- //for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
- //if(k<10){
- //cout << " cells sampling "<< it->second->getSampling() << " : " << it->first << endl;
- //const CaloGeoDetDescrElement *DDE = this->getDDE(it->second->getSampling(),it->second->eta(),it->second->phi());
- //cout << " cells sampling "<< DDE->getSampling() << " : " << DDE->identify() << endl;
- //DDE = this->getDDE(it->first);
- //cout << " cells sampling "<< DDE->getSampling() << " : " << DDE->identify() << endl << endl;
- //}
- //k++;
- //}
-
- //}
-
-
- for(int imodule=1;imodule<=3;imodule++) delete electrodes[imodule-1];
- electrodes.clear();
-
- return true;
-}
-
-void CaloGeometryFromFile::DrawFCalGraph(int isam,int color){
-
- stringstream ss;
- ss << "FCal" << isam - 20 << endl;
- TString name = ss.str().c_str();
-
- const unsigned int size=m_cells_in_sampling[isam].size();
- double* x = new double[size];
- double* y = new double[size];
- //const CaloGeoDetDescrElement* cell;
- unsigned int i=0;
- for(auto it=m_cells_in_sampling[isam].begin();it!=m_cells_in_sampling[isam].end();it++){
- x[i]=it->second->x();
- y[i]=it->second->y();
- i++;
- }
- // cout << size << endl;
- //TH2D* h = new TH2D("","",10,-0.5,0.5,10,-0.5,0.5);
- //h->SetStats(0);
- //h->Draw();
- TGraph* graph = new TGraph(size,x,y);
- graph->SetLineColor(color);
- graph->SetTitle(name);
- graph->SetMarkerStyle(21);
- graph->SetMarkerColor(color);
- graph->SetMarkerSize(0.5);
- graph->GetXaxis()->SetTitle("x");
- graph->GetYaxis()->SetTitle("y");
-
- graph->Draw("AP");
-
- if (x) { delete x; x = 0; }
- if (y) { delete y; y = 0; }
-
-}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FCAL_ChannelMap.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FCAL_ChannelMap.cxx
deleted file mode 100755
index 4c2cbebd5f5d..000000000000
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FCAL_ChannelMap.cxx
+++ /dev/null
@@ -1,488 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-// ***************************************************************************
-// Liquid Argon FCAL detector description package
-// -----------------------------------------
-//
-//
-// 10-Sep-2000 S.Simion Handling of the FCAL read-out identifiers
-// Jan-2001 R.Sobie Modify for persistency
-// Feb-2002 R.Sobie Use same FCAL geometry files as simulation
-//****************************************************************************
-
-#include "ISF_FastCaloSimParametrization/FCAL_ChannelMap.h"
-//#include "CLHEP/Units/SystemOfUnits.h"
-//#include "boost/io/ios_state.hpp"
-#include <sstream>
-#include <iostream>
-#include <iomanip>
-#include <stdio.h>
-
-/* === Geometrical parameters === */
-//const double cm = 0.01;
-const double cm = 10.;
-//const double FCAL_ChannelMap::m_tubeSpacing[] = {0.75*CLHEP::cm, 0.8179*CLHEP::cm, 0.90*CLHEP::cm};
-const double FCAL_ChannelMap::m_tubeSpacing[] = {0.75*cm, 0.8179*cm, 0.90*cm};
-
-FCAL_ChannelMap::FCAL_ChannelMap( int flag)
-{
-
- /* === Initialize geometrical dimensions */
- for(int i=0; i<3; i++){
- m_tubeDx[i] = m_tubeSpacing[i] / 2.;
- m_tubeDy[i] = m_tubeSpacing[i] * sqrt(3.)/2.;
- }
-
- // FCAL1 small cells are 2x2 tubes
- m_tileDx[0] = 2. * m_tubeSpacing[0];
- m_tileDy[0] = 2. * m_tubeSpacing[0] * sqrt(3.)/2.;
-
- // FCAL2 small cells are 2x3 tubes
- m_tileDx[1] = 2. * m_tubeSpacing[1];
- m_tileDy[1] = 3. * m_tubeSpacing[1] * sqrt(3.)/2.;
-
- // FCAL3 cells are 6x6 tubes
- m_tileDx[2] = 6. * m_tubeSpacing[2];
- m_tileDy[2] = 6. * m_tubeSpacing[2] * sqrt(3.)/2.;
-
-
- m_invert_x = flag & 1;
- m_invert_xy = flag & 2;
-
-}
-
-
-void FCAL_ChannelMap::finish() {
- create_tileMap(1);
- create_tileMap(2);
- create_tileMap(3);
-}
-
-// *********************************************************************
-// Read tube mapping tables
-//
-// Jan 23,2002 R. Sobie
-// ********************************************************************
-
-
-//original
-void FCAL_ChannelMap::add_tube(const std::string & tileName, int mod, int /*id*/, int i, int j, double x, double y) {
-
- // Get three integers from the tileName:
- std::istringstream tileStream1(std::string(tileName,1,1));
- std::istringstream tileStream2(std::string(tileName,3,2));
- std::istringstream tileStream3(std::string(tileName,6,3));
- int a1=0,a2=0,a3=0;
- if (tileStream1) tileStream1 >> a1;
- if (tileStream2) tileStream2 >> a2;
- if (tileStream3) tileStream3 >> a3;
-
- tileName_t tilename = (a3 << 16) + a2;
-
- TubePosition tb(tilename, x*cm, y*cm,"");
- // Add offsets, becaues iy and ix can be negative HMA
-
- i = i+200;
- j = j+200;
- // m_tubeMap[mod-1][(j << 16) + i] = tb;
- unsigned int ThisId = (j<<16) + i;
- tubemap_const_iterator p = m_tubeMap[mod-1].insert(m_tubeMap[mod-1].end(),std::make_pair(ThisId,tb));
- m_tubeIndex[mod-1].push_back(p);
-}
-
-
-//Gabe: new to include HV and LARFCALELECRTODES ID
-void FCAL_ChannelMap::add_tube(const std::string & tileName, int mod, int /*id*/, int i, int j, double x, double y, std::string hvFT) {
-
- // Get three integers from the tileName:
- std::istringstream tileStream1(std::string(tileName,1,1));
- std::istringstream tileStream2(std::string(tileName,3,2));
- std::istringstream tileStream3(std::string(tileName,6,3));
- int a1=0,a2=0,a3=0;
- if (tileStream1) tileStream1 >> a1;
- if (tileStream2) tileStream2 >> a2;
- if (tileStream3) tileStream3 >> a3;
-
- tileName_t tilename = (a3 << 16) + a2;
-
- TubePosition tb(tilename, x*cm,y*cm, hvFT);
- // Add offsets, becaues iy and ix can be negative HMA
-
- i = i+200;
- j = j+200;
- // m_tubeMap[mod-1][(j << 16) + i] = tb;
- unsigned int ThisId = (j<<16) + i;
- tubemap_const_iterator p = m_tubeMap[mod-1].insert(m_tubeMap[mod-1].end(),std::make_pair(ThisId,tb));
- m_tubeIndex[mod-1].push_back(p);
-}
-
-FCAL_ChannelMap::tubemap_const_iterator FCAL_ChannelMap::getTubeByCopyNumber( int isam, int copyNo) const {
- return m_tubeIndex[isam-1][copyNo];
-}
-
-
-// *********************************************************************
-// Create tile mapping tables
-//
-// Jan 23,2002 R. Sobie
-// ********************************************************************
-void
-FCAL_ChannelMap::create_tileMap(int isam)
-{
- tileMap_const_iterator tile;
- tubemap_const_iterator first = m_tubeMap[isam-1].begin();
- tubemap_const_iterator last = m_tubeMap[isam-1].end();
-
- // Loop over tubes -> find unique tiles and fill the descriptors
- while (first != last){
-
- tileName_t tileName = (first->second).get_tileName();
- tile = m_tileMap[isam-1].find(tileName);
-
- if (tile == m_tileMap[isam-1].end()){ // New tile found
- float x = (first->second).x();
- float y = (first->second).y();
- unsigned int ntubes = 1;
- TilePosition tp(x, y, ntubes);
- m_tileMap[isam-1][tileName] = tp;
- }
- else{ // Existing tile
- float x = (tile->second).x() + (first->second).x();
- float y = (tile->second).y() + (first->second).y();
- unsigned int ntubes = (tile->second).ntubes() + 1;
- TilePosition tp(x, y, ntubes);
- m_tileMap[isam-1][tileName] = tp;
- }
- ++first;
- }
-
- //
- // List the number of tubes and tiles
- //
- // std::cout << "FCAL::create_tilemap: FCAL # " << isam
- // << " Number of tiles = " << m_tileMap[isam-1].size()
- // << " Number of tubes = " << m_tubeMap[isam-1].size()
- // << std::endl;
-
- // this->print_tubemap(isam);
-
-
- //
- // loop over tiles and set (x,y) to average tile positions
- //
- tileMap_const_iterator tilefirst = m_tileMap[isam-1].begin();
- tileMap_const_iterator tilelast = m_tileMap[isam-1].end();
- while (tilefirst != tilelast) {
- tileName_t tileName = tilefirst->first;
- unsigned int ntubes = (tilefirst->second).ntubes();
- float xtubes = (float) ntubes;
- float x = (tilefirst->second).x() / xtubes;
- float y = (tilefirst->second).y() / xtubes;
- TilePosition tp(x, y, ntubes);
- m_tileMap[isam-1][tileName] = tp;
- ++tilefirst;
- }
-
-}
-
-//---------- for New LArFCAL_ID ------------------------
-
-// *********************************************************************
-// get Tile ID
-//
-// Original code: Stefan Simion, Randy Sobie
-// -------------------------------------------------------------------
-// This function computes the tile identifier for any given position
-// Inputs:
-// - isam the sampling number, from G3 data;
-// - x the tube x position, in CLHEP::cm, from G3 data;
-// - y the tube y position, in CLHEP::cm, from G3 data.
-// Outputs:
-// - pair of indices eta, phi
-//
-// Attention side-effect: x is changed by this function.
-// --------------------------------------------------------------------
-// June 2002 HMA
-// ***********************************************************************
-bool
-FCAL_ChannelMap::getTileID(int isam, float x_orig, float y_orig,
- int& eta, int& phi) const throw (std::range_error)
-{
-
-// /* ### MIRROR for compatibility between G3 and ASCII files ### */
-
- float x = x_orig;
- float y = y_orig;
-
- if(m_invert_xy){
- x = y_orig;
- y = x_orig;
- }
-
- if(m_invert_x) x = -x;
-
- /* === Compute the tubeID */
- int ktx = (int) (x / m_tubeDx[isam-1]);
- int kty = (int) (y / m_tubeDy[isam-1]);
- if (x < 0.) ktx--;
- if (y < 0.) kty--;
-
- // S.M.: in case we lookup real positions inside the Tile (not only
- // integer grids for the tubes) half of the time we are outisde a
- // tube bin since the integer pattern of the tubes is like in this
- // sketch:
- //
- // # # # #
- // # # # #
- // # # # #
- //
- // in order to avoid this problem we have to make sure the integer
- // indices for x and y have either both to be even or both to be odd
- // (For Module 0 one has to be odd the other even ...). We take the
- // y-index and check for odd/even and change the x-index in case
- // it's different from the first tube in the current sampling ...
- //
- // S.M. update: in case we are in a hole of the integer grid the
- // relative x and y w.r.t to the original tube are used to assign a
- // tube according to the hexagonal pattern.
-
- tubemap_const_iterator it = m_tubeMap[isam-1].begin();
- unsigned int firstId = it->first;
-
- // take offset from actual map
- int ix = ktx+((int)((firstId&0xffff)-it->second.x()/m_tubeDx[isam-1]))+1;
- int iy = kty+((int)((firstId>>16)-it->second.y()/m_tubeDy[isam-1]))+1;
-
- int isOddEven = (((firstId>>16)%2)+(firstId%2))%2;
- bool movex = false;
-
- if ( (iy%2) != ((ix+isOddEven)%2) ) {
- double yc = y/m_tubeDy[isam-1] - kty - 0.5;
- if ( fabs(yc) > 0.5/sqrt(3) ) {
- double xk = x/m_tubeDx[isam-1] - ktx;
- if ( xk > 0.5 ) {
- xk = 1 - xk;
- }
- double yn = 0.5-xk/3;
- if ( fabs(yc) > fabs(yn) ) {
- if ( yc > 0 )
- iy++;
- else
- iy--;
- }
- else
- movex = true;
- }
- else
- movex = true;
- if ( movex ) {
- if ( x/m_tubeDx[isam-1] - ktx > 0.5 )
- ix++;
- else
- ix--;
- }
- }
-
- tubeID_t tubeID = (iy << 16) + ix;
-
- it = m_tubeMap[isam-1].find(tubeID);
- if (it != m_tubeMap[isam-1].end()){
- tileName_t tilename = (it->second).get_tileName();
- phi = tilename & 0xffff;
- eta = tilename >> 16;
- return true ;
- }
- // reach here only if it failed the second time.
-
- return false;
-
-}
-
-
-
-/* ----------------------------------------------------------------------
- To decode the tile x position from the tile identifier
- ---------------------------------------------------------------------- */
-float
-FCAL_ChannelMap::x(int isam, int eta, int phi) const
- throw(std::range_error)
-{
- if(m_invert_xy){
- // temp turn off the flag
- m_invert_xy=false;
- float y1 = y(isam,eta,phi);
- m_invert_xy=true;
- return y1;
- }
- float x;
-
- tileName_t tilename = (eta << 16) + phi ;
-
- tileMap_const_iterator it = m_tileMap[isam-1].find(tilename);
- if(it != m_tileMap[isam-1].end())
- {
- x = (it->second).x();
- } else
- { // can't find the tile, throw exception.
- char l_str[200] ;
- snprintf(l_str, sizeof(l_str),
- "Error in FCAL_ChannelMap::x, wrong tile,phi= %d ,eta=: %d ",phi,eta);
- std::string errorMessage(l_str);
- throw std::range_error(errorMessage.c_str());
- }
-
- if(m_invert_x) {
- return -x;
- }
- else {
- return x;
- }
-
- return x;
-
-}
-
-
-/* ----------------------------------------------------------------------
- To decode the tile y position from the tile identifier
- ---------------------------------------------------------------------- */
-float
-FCAL_ChannelMap::y(int isam, int eta, int phi) const
- throw(std::range_error)
-{
- if(m_invert_xy){
-
- // temp turn off the flag
- m_invert_xy=false;
- float x1 = x(isam,eta,phi);
- m_invert_xy=true;
- return x1;
-
- }
-
- float y;
-
- tileName_t tilename = (eta << 16) + phi ;
-
- tileMap_const_iterator it = m_tileMap[isam-1].find(tilename);
- if(it != m_tileMap[isam-1].end())
- {
- y = (it->second).y();
- } else
- { // can't find the tile, throw exception.
- char l_str[200] ;
- snprintf(l_str, sizeof(l_str),
- "Error in FCAL_ChannelMap::x, wrong tile,phi= %d ,eta=: %d",phi,eta);
- std::string errorMessage(l_str);
- throw std::range_error(errorMessage.c_str());
- }
-
- return y;
-}
-
-/* ----------------------------------------------------------------------
- To decode the tile dx size from the tile identifier
- ---------------------------------------------------------------------- */
-
-void FCAL_ChannelMap::tileSize(int sam, int ntubes, float &dx, float &dy) const {
-
- dx = m_tubeDx[sam-1];
- dy = m_tubeDy[sam-1];
- // float ntubes = (it->second).ntubes();
- if(sam == 1 || sam == 3) {
- float scale =sqrt(ntubes);
- dx = dx * scale;
- dy = dy * scale;
- }
- else {
- float scale = sqrt(ntubes/1.5);
- dx = dx * scale;
- dy = dy * scale * 1.5 ;
- }
-
-
- // There is a fundamental discrepancy between dx and dy. A cell will
- // contain twice as many distinct x-positions as y-positions. Diagram:
-
- // . . . . -
- //. . . . -
- // . . . . - 4 x dy
- // . . . . -
- // ||||||||
- // 8 x dx
-
- dx = 2*dx;
-
- if(m_invert_xy){
- // switch xy
- float temp = dx;
- dx = dy;
- dy = temp;
- }
-
-}
-
-void FCAL_ChannelMap::tileSize(int sam, int eta, int phi,
- float& dx, float& dy ) const throw(std::range_error)
-{
-
- tileName_t tilename = (eta << 16) + phi ;
-
- tileMap_const_iterator it = m_tileMap[sam-1].find(tilename);
- if(it != m_tileMap[sam-1].end()) {
- int ntubes = (it->second).ntubes();
- tileSize(sam,ntubes,dx,dy);
- return ;
- }
- else {
- // can't find the tile, throw exception.
- char l_str[200] ;
- snprintf(l_str, sizeof(l_str),
- "Error in FCAL_ChannelMap::tilesize, wrong tile,phi= %d ,eta=: %d ",phi,eta);
- std::string errorMessage(l_str);
- throw std::range_error(errorMessage.c_str());
- }
-}
-
-
-// ********************** print_tubemap *************************************
-void
-FCAL_ChannelMap::print_tubemap( int imap) const
-{
- FCAL_ChannelMap::tubemap_const_iterator it = m_tubeMap[imap-1].begin();
-
- //boost::io::ios_all_saver ias (std::cout);
- std::cout << "First 10 elements of the New FCAL tube map : " << imap << std::endl;
- std::cout.precision(5);
- for ( int i=0; i<10; i++, it++)
- std::cout << std::hex << it->first << "\t"
- << (it->second).get_tileName() << std::dec <<"\t"
- << (it->second).x() <<"\t"
- << (it->second).y() << std::endl;
-
-}
-
-
-// ********************** tubemap_begin **************************************
-FCAL_ChannelMap::tubemap_const_iterator
-FCAL_ChannelMap::tubemap_begin (int imap ) const
-{
- return m_tubeMap[imap-1].begin();
-}
-
-
-// ********************** tubemap_end ****************************************
-FCAL_ChannelMap::tubemap_const_iterator
-FCAL_ChannelMap::tubemap_end (int imap ) const
-{
- return m_tubeMap[imap-1].end();
-}
-
-// ********************** tubemap_size ***************************************
-FCAL_ChannelMap::tubemap_sizetype
-FCAL_ChannelMap::tubemap_size (int imap) const
-{
- return m_tubeMap[imap-1].size();
-}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FastCaloSimParamAlg.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FastCaloSimParamAlg.cxx
index 3d57daf1ab87..0e98d1b15762 100755
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FastCaloSimParamAlg.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/FastCaloSimParamAlg.cxx
@@ -193,7 +193,7 @@ StatusCode FastCaloSimParamAlg::execute()
}
//clusterize(eventSteps);
std::cout << "Size after clusterization: " << MergeeventSteps->size() << std::endl;
- StatusCode sc = evtStore()->record(MergeeventSteps,"MergedEventSteps");
+ StatusCode sc = evtStore()->record(MergeeventSteps,"ZHMergedEventSteps");
if (sc.isFailure())
{
std::cout<<"Coudn't resave merged collection "<<std::endl;
@@ -219,7 +219,7 @@ const HepMC::GenParticle* FastCaloSimParamAlg::getParticleFromMC()
const ISF_FCS_Parametrization::FCS_StepInfoCollection* FastCaloSimParamAlg::getFCS_StepInfo()
{
const ISF_FCS_Parametrization::FCS_StepInfoCollection* eventStepsES;
- StatusCode sc = evtStore()->retrieve(eventStepsES, "EventSteps");
+ StatusCode sc = evtStore()->retrieve(eventStepsES, "ZHEventSteps");
if (sc.isFailure()) return NULL;
return eventStepsES;
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 d1152721b01c..40083e8c512f 100755
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/ISF_HitAnalysis.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/ISF_HitAnalysis.cxx
@@ -779,7 +779,7 @@ StatusCode ISF_HitAnalysis::execute()
//Get the FastCaloSim step info collection from store
const ISF_FCS_Parametrization::FCS_StepInfoCollection* eventStepsES;
- StatusCode sc = evtStore()->retrieve(eventStepsES, "MergedEventSteps");
+ StatusCode sc = evtStore()->retrieve(eventStepsES, "ZHMergedEventSteps");
if (sc.isFailure())
{
ATH_MSG_WARNING( "No FastCaloSim steps read from StoreGate?" );
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloGeometryFromFile.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloGeometryFromFile.cxx
index a1c38db61ff2..28dd5b70d6de 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloGeometryFromFile.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloGeometryFromFile.cxx
@@ -5,7 +5,7 @@
#include "CaloGeometryFromFile.h"
#include <TTree.h>
#include <TFile.h>
-#include "CaloDetDescr/CaloDetDescrElement.h"
+#include "ISF_FastCaloSimParametrization/CaloGeoDetDescrElement.h"
#include <fstream>
#include <sstream>
#include <TGraph.h>
@@ -13,8 +13,28 @@
using namespace std;
-CaloGeometryFromFile::CaloGeometryFromFile() : CaloGeometry()
+map<unsigned long long, unsigned long long> g_cellId_vs_cellHashId_map;
+
+CaloGeometryFromFile::CaloGeometryFromFile() : CaloGeoGeometry()
{
+ ifstream textfile("/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/cellId_vs_cellHashId_map.txt");
+ unsigned long long id, hash_id;
+ cout << "Loading cellId_vs_cellHashId_map" << endl;
+ int i=0;
+ string line;
+ stringstream s;
+ while(1){
+ //getline(textfile,line);
+ s.str(line);
+ textfile >> id;
+ if(textfile.eof())break;
+ textfile >> hash_id;
+ g_cellId_vs_cellHashId_map[id]=hash_id;
+ if(i%10000==0)cout << "Line: " << i << " line " << line << " id " << hex << id << " hash_id " << dec << hash_id << endl;
+ i++;
+ }
+ cout << "Done." << endl;
+
}
CaloGeometryFromFile::~CaloGeometryFromFile()
@@ -83,14 +103,20 @@ bool CaloGeometryFromFile::LoadGeometryFromFile(TString filename,TString treenam
if (ientry < 0) break;
fChain->GetEntry(jentry);
- CaloGeoDetDescrElement* pcell=new CaloGeoDetDescrElement(cell);
+ if (g_cellId_vs_cellHashId_map.find(cell.m_identify)!=g_cellId_vs_cellHashId_map.end()) cell.m_hash_id=g_cellId_vs_cellHashId_map[cell.m_identify];
+ else cout << endl << "ERROR: Cell id not found in the cellId_vs_cellHashId_map!!!" << endl << endl;
+
+ const CaloGeoDetDescrElement* pcell=new CaloGeoDetDescrElement(cell);
int sampling = pcell->getSampling();
+
+
if(sampling >20) continue;
- addcell(pcell);
-
- //if(jentry%25000==0) {
- if(jentry==nentries-1) {
- //cout<<jentry<<" : "<<cell.getSampling()<<", "<<cell.identify()<<endl;
+ this->addcell(pcell);
+
+ if(jentry%25000==0) {
+ //if(jentry==nentries-1) {
+ cout << "Checking loading cells from file" << endl;
+ cout<<jentry<<" : "<<pcell->getSampling()<<", "<<pcell->identify()<<endl;
@@ -99,25 +125,41 @@ bool CaloGeometryFromFile::LoadGeometryFromFile(TString filename,TString treenam
}
}
//cout<<"all : "<<m_cells.size()<<endl;
- Long64_t max(0);
- Long64_t min_id=m_cells_in_sampling[0].begin()->first;
- for(int i=0;i<21;++i) {
- //cout<<" cells sampling "<<i<<" : "<<m_cells_in_sampling[i].size()<<" cells";
- //cout<<", "<<m_cells_in_regions[i].size()<<" lookup map(s)"<<endl;
- for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
- //cout << it->second->getSampling() << " " << it->first << endl;
- if(min_id/10 >= it->first){
- //cout << "Warning: Identifiers are not in increasing order!!!!" << endl;
- //cout << min_id << " " << it->first << endl;
+ //unsigned long long max(0);
+ //unsigned long long min_id=m_cells_in_sampling[0].begin()->first;
+ //for(int i=0;i<21;++i) {
+ ////cout<<" cells sampling "<<i<<" : "<<m_cells_in_sampling[i].size()<<" cells";
+ ////cout<<", "<<m_cells_in_regions[i].size()<<" lookup map(s)"<<endl;
+ //for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
+ ////cout << it->second->getSampling() << " " << it->first << endl;
+ //if(min_id/10 >= it->first){
+ ////cout << "Warning: Identifiers are not in increasing order!!!!" << endl;
+ ////cout << min_id << " " << it->first << endl;
- }
- if(min_id > it->first)min_id = it->first;
- }
- }
+ //}
+ //if(min_id > it->first)min_id = it->first;
+ //}
+ //}
//cout << "Min id for samplings < 21: " << min_id << endl;
delete f;
//return true;
- return PostProcessGeometry();
+ bool ok = PostProcessGeometry();
+
+ cout << "Result of PostProcessGeometry(): " << ok << endl;
+ const CaloGeoDetDescrElement* mcell=0;
+ unsigned long long cellid64(3179554531063103488);
+ Identifier cellid(cellid64);
+ mcell=this->getDDE(cellid); //This is working also for the FCal
+
+ std::cout << "\n \n";
+ std::cout << "Testing whether CaloGeoGeometry is loaded properly" << std::endl;
+ if(!mcell)std::cout << "Cell is not found" << std::endl;
+ std::cout << "Identifier " << mcell->identify() <<" sampling " << mcell->getSampling() << " eta: " << mcell->eta() << " phi: " << mcell->phi() << std::endl<< std::endl;
+
+ const CaloGeoDetDescrElement* mcell2 = this->getDDE(mcell->getSampling(),mcell->eta(),mcell->phi());
+ std::cout << "Identifier " << mcell2->identify() <<" sampling " << mcell2->getSampling() << " eta: " << mcell2->eta() << " phi: " << mcell2->phi() << std::endl<< std::endl;
+
+ return ok;
}
@@ -134,12 +176,12 @@ bool CaloGeometryFromFile::LoadFCalGeometryFromFiles(TString filename1,TString f
int thisTubeI;
int thisTubeJ;
//int thisTubeID;
- int thisTubeMod;
+ //int thisTubeMod;
double thisTubeX;
double thisTubeY;
TString tubeName;
- int second_column;
+ //int second_column;
string seventh_column;
string eight_column;
int ninth_column;
@@ -178,7 +220,7 @@ bool CaloGeometryFromFile::LoadFCalGeometryFromFiles(TString filename1,TString f
if (tileStream2) tileStream2 >> a2;
if (tileStream3) tileStream3 >> a3;
- unsigned int tileName= (a3 << 16) + a2;
+ //unsigned int tileName= (a3 << 16) + a2;
stringstream s;
@@ -192,14 +234,16 @@ bool CaloGeometryFromFile::LoadFCalGeometryFromFiles(TString filename1,TString f
i++;
}
}
+
+
m_FCal_ChannelMap.finish(); // Creates maps
- Long64_t phi_index,eta_index;
+ unsigned long long phi_index,eta_index;
float x,y,dx,dy;
long id;
//auto it =m_cells_in_sampling[20].rbegin();
//it--;
- //Long64_t identify=it->first;
+ //unsigned long long identify=it->first;
//for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
//
//}
@@ -240,28 +284,49 @@ bool CaloGeometryFromFile::LoadFCalGeometryFromFiles(TString filename1,TString f
DDE->m_identify=(id << 44);
//DDE->m_identify=(id << 12);
//cout << DDE->identify() << endl;
+ if (g_cellId_vs_cellHashId_map.find(DDE->m_identify)!=g_cellId_vs_cellHashId_map.end()) DDE->m_hash_id=g_cellId_vs_cellHashId_map[DDE->m_identify];
+ else cout << endl << "ERROR: Cell id not found in the cellId_vs_cellHashId_map!!!" << endl << endl;
+
addcell(DDE);
id=(mask2[imap-1]<<12) + (eta_index << 5) +2*phi_index;
if(imap==2) id+= (8<<8);
DDE2->m_identify=(id << 44);
//DDE->m_identify=(id << 12);
+ if (g_cellId_vs_cellHashId_map.find(DDE2->m_identify)!=g_cellId_vs_cellHashId_map.end()) DDE2->m_hash_id=g_cellId_vs_cellHashId_map[DDE2->m_identify];
+ else cout << endl << "ERROR: Cell id not found in the cellId_vs_cellHashId_map!!!" << endl << endl;
addcell(DDE2);
+ delete DDE;
+ delete DDE2;
}
}
- auto it =m_cells_in_sampling[0].begin();
+ //auto it =m_cells_in_sampling[0].begin();
//it--;
- Long64_t identify=it->first;
- //cout << "min identifier from sampling < 21: " << identify << endl;
+ //unsigned long long identify=it->first;
+ ////cout << "min identifier from sampling < 21: " << identify << endl;
- for(int i=21;i<MAX_SAMPLING;++i) {
+ //for(int i=0;i<21;++i) {
+ ////for(int i=21;i<MAX_SAMPLING;++i) {
//cout<<" cells sampling "<<i<<" : "<<m_cells_in_sampling[i].size()<<" cells";
//cout<<", "<<m_cells_in_regions[i].size()<<" lookup map(s)"<<endl;
- for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
- //cout << it->second->getSampling() << " " << it->first << endl;
- }
+ //int k=0;
+ //for(auto it=m_cells_in_sampling[i].begin(); it!=m_cells_in_sampling[i].end();it++){
+ //if(k<10){
+ //cout << " cells sampling "<< it->second->getSampling() << " : " << it->first << endl;
+ //const CaloGeoDetDescrElement *DDE = this->getDDE(it->second->getSampling(),it->second->eta(),it->second->phi());
+ //cout << " cells sampling "<< DDE->getSampling() << " : " << DDE->identify() << endl;
+ //DDE = this->getDDE(it->first);
+ //cout << " cells sampling "<< DDE->getSampling() << " : " << DDE->identify() << endl << endl;
+ //}
+ //k++;
+ //}
- }
+ //}
+
+
+ for(int imodule=1;imodule<=3;imodule++) delete electrodes[imodule-1];
+ electrodes.clear();
+
return true;
}
@@ -274,7 +339,7 @@ void CaloGeometryFromFile::DrawFCalGraph(int isam,int color){
const int size=m_cells_in_sampling[isam].size();
double x[size];
double y[size];
- const CaloGeoDetDescrElement* cell;
+ //const CaloGeoDetDescrElement* cell;
int i=0;
for(auto it=m_cells_in_sampling[isam].begin();it!=m_cells_in_sampling[isam].end();it++){
x[i]=it->second->x();
@@ -291,12 +356,13 @@ void CaloGeometryFromFile::DrawFCalGraph(int isam,int color){
graph->SetMarkerStyle(21);
graph->SetMarkerColor(color);
graph->SetMarkerSize(0.5);
- graph->GetXaxis()->SetTitle("x");
- graph->GetYaxis()->SetTitle("y");
+ graph->GetXaxis()->SetTitle("x [mm]");
+ graph->GetYaxis()->SetTitle("y [mm]");
graph->Draw("AP");
}
+
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C
index f3c901f18325..3a71c69b2933 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/init_geo.C
@@ -8,7 +8,7 @@
#include "../ISF_FastCaloSimParametrization/MeanAndRMS.h"
#include "Identifier/Identifier.h"
-#include "CaloDetDescr/CaloDetDescrElement.h"
+//#include "CaloDetDescr/CaloDetDescrElement.h"
void init_geo();
@@ -21,9 +21,9 @@ void init_geo()
gInterpreter->AddIncludePath("..");
gInterpreter->AddIncludePath("../../ISF_FastCaloSimEvent");
- gROOT->LoadMacro("CaloSampling.cxx+");
- gROOT->LoadMacro("../src/CaloGeometry.cxx+");
- gROOT->LoadMacro("CaloGeometryFromFile.cxx+");
+ //gROOT->LoadMacro("CaloSampling.cxx+");
+ gROOT->LoadMacro("../src/CaloGeoGeometry.cxx+");
+ gROOT->LoadMacro("../src/CaloGeometryFromFile.cxx+");
gROOT->LoadMacro("../src/FCAL_ChannelMap.cxx+");
cout<<"init geometry done"<<endl;
cout << "running run_geo.C" << endl;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C
index b6ec56c5da96..ca975dbf459a 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/run_geo.C
@@ -7,6 +7,10 @@
#include "TGraph.h"
void run_geo();
+
+void WriteGeneralInfo(TString /*cut_label*/, TString lumi, float size, float x, float y);
+void ATLASLabel(Double_t x,Double_t y,const char* text, float tsize, Color_t color=1);
+void WriteInfo(TString info, float size, float x, float y, int color=1);
void run_geo()
{
@@ -54,15 +58,27 @@ void run_geo()
geo->DrawGeoForPhi0();
canvas->SaveAs("Calorimeter.png");*/
- TCanvas* canvas = new TCanvas("FCal1_xy","FCal1_xy");
+ float xInfo = 0.18;
+ float yInfo = 0.9;
+ float sizeInfo=0.035;
+
+
+
+ TCanvas* canvas = new TCanvas("FCal1_xy","FCal1_xy",600,600);
geo->DrawFCalGraph(21,1);
- canvas->SaveAs("FCal1Geometry.png");
- TCanvas* canvas2 = new TCanvas("FCal2_xy","FCal2_xy");
+ WriteGeneralInfo("","",sizeInfo,xInfo,yInfo);
+ WriteInfo("FCal1", 0.05, 0.20, 0.21);
+ canvas->SaveAs("FCal1Geometry.png");
+ TCanvas* canvas2 = new TCanvas("FCal2_xy","FCal2_xy",600,600);
geo->DrawFCalGraph(22,1);
- canvas2->SaveAs("FCal2Geometry.png");
- TCanvas* canvas3 = new TCanvas("FCal3_xy","FCal3_xy");
+ WriteGeneralInfo("","",sizeInfo,xInfo,yInfo);
+ WriteInfo("FCal2", 0.05, 0.20, 0.21);
+ canvas2->SaveAs("FCal2Geometry.png");
+ TCanvas* canvas3 = new TCanvas("FCal3_xy","FCal3_xy",600,600);
geo->DrawFCalGraph(23,1);
- canvas3->SaveAs("FCal3Geometry.png");
+ WriteGeneralInfo("","",sizeInfo,xInfo,yInfo);
+ WriteInfo("FCal3", 0.05, 0.20, 0.21);
+ canvas3->SaveAs("FCal3Geometry.png");
@@ -190,4 +206,58 @@ void run_geo()
}
-
+void ATLASLabel(Double_t x,Double_t y,const char* text, float tsize, Color_t color){
+ TLatex l; //l.SetTextAlign(12);
+ if (tsize>0) l.SetTextSize(tsize);
+ l.SetNDC();
+ l.SetTextFont(72);
+ l.SetTextColor(color);
+
+ //double delx = 0.115*696*gPad->GetWh()/(472*gPad->GetWw());
+ double delx = 0.14;
+
+ l.DrawLatex(x,y,"ATLAS");
+ if (text) {
+ TLatex p;
+ p.SetNDC();
+ if (tsize>0) p.SetTextSize(tsize);
+ p.SetTextFont(42);
+ p.SetTextColor(color);
+ p.DrawLatex(x+delx,y,text);
+ // p.DrawLatex(x,y,"#sqrt{s}=900GeV");
+ }
+}
+void WriteGeneralInfo(TString /*cut_label*/, TString lumi, float size, float x, float y){
+ TString label="";
+ if (lumi=="") label+=" Simulation";
+ //label+=" Internal";
+ label+=" Preliminary";
+ ATLASLabel(x,y,label.Data(),size*1.15);
+ TString ToWrite="";
+ TLatex l;
+ l.SetNDC();
+ l.SetTextFont(42);
+ l.SetTextSize(size*0.9);
+ l.SetTextColor(1);
+ //l.DrawLatex(x-0.005,y-0.07,cut_label.Data());
+
+ double shift=0.55;
+ //ToWrite="#sqrt{s}=13 TeV";
+ //if (lumi!=""){
+ ////ToWrite="L_{int}=";
+ //ToWrite+=", ";
+ //ToWrite+=lumi;
+ //ToWrite+=" fb^{-1}";
+ //shift=0.43;
+ //}
+ l.DrawLatex(x+shift,y,ToWrite.Data());
+
+}
+void WriteInfo(TString info, float size, float x, float y, int color){
+ TLatex l;
+ l.SetNDC();
+ l.SetTextFont(42);
+ l.SetTextSize(size);
+ l.SetTextColor(color);
+ l.DrawLatex(x,y,info.Data());
+}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt
index d222b374cfb6..c0dbcc846ea0 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt
@@ -9,8 +9,6 @@ atlas_subdir( ISF_FastCaloSimServices )
atlas_depends_on_subdirs( PUBLIC
GaudiKernel
PRIVATE
- Control/AthenaKernel
- DetectorDescription/IdDictParser
Control/AthenaBaseComps
Simulation/Barcode/BarcodeEvent
Simulation/ISF/ISF_Core/ISF_Interfaces
@@ -18,14 +16,14 @@ atlas_depends_on_subdirs( PUBLIC
Tracking/TrkExtrapolation/TrkExInterfaces
Calorimeter/CaloInterface
Calorimeter/CaloEvent
+ Calorimeter/CaloDetDescr
Control/StoreGate
Event/NavFourMom
Generators/GeneratorObjects
Simulation/FastShower/FastCaloSim
Simulation/ISF/ISF_Core/ISF_Event
Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent
- Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimInterfaces
- Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization )
+ Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimInterfaces )
# External dependencies:
find_package( CLHEP )
@@ -36,7 +34,7 @@ atlas_add_component( ISF_FastCaloSimServices
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS}
- LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaBaseComps AthenaKernel GaudiKernel IdDictParser ISF_Interfaces TrkEventPrimitives TrkExInterfaces CaloEvent StoreGateLib SGtests NavFourMom GeneratorObjects FastCaloSimLib ISF_Event ISF_FastCaloSimEvent ISF_FastCaloSimInterfaces ISF_FastCaloSimParametrizationLib )
+ LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaBaseComps GaudiKernel ISF_Interfaces TrkEventPrimitives TrkExInterfaces CaloEvent StoreGateLib SGtests NavFourMom GeneratorObjects FastCaloSimLib ISF_Event ISF_FastCaloSimEvent ISF_FastCaloSimInterfaces )
# Install files from the package:
atlas_install_headers( ISF_FastCaloSimServices )
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/AdditionalConfig.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/AdditionalConfig.py
index 3cacab794fb4..7e7d36b7dd02 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/AdditionalConfig.py
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/AdditionalConfig.py
@@ -13,11 +13,9 @@ from AthenaCommon.Constants import * # FATAL,ERROR etc.
from AthenaCommon.SystemOfUnits import *
from AthenaCommon.DetFlags import DetFlags
+from ISF_Config.ISF_jobProperties import ISF_Flags # IMPORTANT: Flags must be set before tools are retrieved
from ISF_FastCaloSimServices.ISF_FastCaloSimJobProperties import ISF_FastCaloSimFlags
-from ISF_Algorithms.collection_merger_helpers import generate_mergeable_collection_name
-
-
def getAdditionalParticleParametrizationFileNames():
return [
"DB=/GLOBAL/AtlfastII/FastCaloSimParam:2:EnergyResults/pdgid_211/EN_1000/eta_central",
@@ -674,7 +672,7 @@ def getPunchThroughTool(name="ISF_PunchThroughTool", **kwargs):
kwargs.setdefault("MinEnergy" , [ 938.3, 135.6, 50., 50., 105.7 ] )
kwargs.setdefault("MaxNumParticles" , [ -1, -1, -1, -1, -1 ] )
kwargs.setdefault("EnergyFactor" , [ 1., 1., 1., 1., 1. ] )
- kwargs.setdefault("BarcodeSvc" , simFlags.TruthStrategy.BarcodeServiceName() )
+ kwargs.setdefault("BarcodeSvc" , ISF_Flags.BarcodeService() )
kwargs.setdefault("EnvelopeDefSvc" , getService('AtlasGeometry_EnvelopeDefSvc') )
kwargs.setdefault("BeamPipeRadius" , 500. )
@@ -703,6 +701,16 @@ def getNITimedExtrapolator(name="ISF_NITimedExtrapolator", **kwargs):
from TrkExTools.TrkExToolsConf import Trk__TimedExtrapolator as TimedExtrapolator
return TimedExtrapolator(name, **kwargs )
+
+def getTimedExtrapolator(name="TimedExtrapolator", **kwargs):
+ kwargs.setdefault("MaterialEffectsUpdators" , [ 'ISF_NIMatEffUpdator' ])
+ kwargs.setdefault("ApplyMaterialEffects" , False )
+ kwargs.setdefault("STEP_Propagator" , 'ISF_NIPropagator' )
+
+ from TrkExTools.TrkExToolsConf import Trk__TimedExtrapolator as TimedExtrapolator
+ return TimedExtrapolator(name, **kwargs )
+
+
## FastShowerCellBuilderTool
def getDefaultFastShowerCellBuilderTool(name, **kwargs):
@@ -791,50 +799,7 @@ def getPileupFastShowerCellBuilderTool(name="ISF_PileupFastShowerCellBuilderTool
kwargs.setdefault("sampling_energy_reweighting", weightsfcs )
return getFastShowerCellBuilderTool(name, **kwargs)
-def getFastHitConvertTool(name="ISF_FastHitConvertTool", **kwargs):
- mergeable_collection_suffix = "_FastCaloSim"
-
- EMB_hits_bare_collection_name = "LArHitEMB"
- EMB_hits_merger_input_property = "LArEMBHits"
- EMB_hits_collection_name = generate_mergeable_collection_name(
- EMB_hits_bare_collection_name,
- mergeable_collection_suffix,
- EMB_hits_merger_input_property)
-
- EMEC_hits_bare_collection_name = "LArHitEMEC"
- EMEC_hits_merger_input_property = "LArEMECHits"
- EMEC_hits_collection_name = generate_mergeable_collection_name(
- EMEC_hits_bare_collection_name,
- mergeable_collection_suffix,
- EMEC_hits_merger_input_property)
-
- FCAL_hits_bare_collection_name = "LArHitFCAL"
- FCAL_hits_merger_input_property = "LArFCALHits"
- FCAL_hits_collection_name = generate_mergeable_collection_name(
- FCAL_hits_bare_collection_name,
- mergeable_collection_suffix,
- FCAL_hits_merger_input_property)
-
- HEC_hits_bare_collection_name = "LArHitHEC"
- HEC_hits_merger_input_property = "LArHECHits"
- HEC_hits_collection_name = generate_mergeable_collection_name(
- HEC_hits_bare_collection_name,
- mergeable_collection_suffix,
- HEC_hits_merger_input_property)
-
- tile_hits_bare_collection_name = "TileHitVec"
- tile_hits_merger_input_property = "TileHits"
- tile_hits_collection_name = generate_mergeable_collection_name(
- tile_hits_bare_collection_name,
- mergeable_collection_suffix,
- tile_hits_merger_input_property)
-
- kwargs.setdefault('embHitContainername', EMB_hits_collection_name)
- kwargs.setdefault('emecHitContainername', EMEC_hits_collection_name)
- kwargs.setdefault('fcalHitContainername', FCAL_hits_collection_name)
- kwargs.setdefault('hecHitContainername', HEC_hits_collection_name)
- kwargs.setdefault('tileHitContainername', tile_hits_collection_name)
-
+def getFastHitConvertTool(name="ISF_FastHitConvertTool",**kwargs):
from FastCaloSimHit.FastCaloSimHitConf import FastHitConvertTool
return FastHitConvertTool(name,**kwargs)
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py
index 48621e9ef886..2797c9e57cba 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py
@@ -96,14 +96,33 @@ def getFastCaloSimSvcV2(name="ISF_FastCaloSimSvcV2", **kwargs):
kwargs.setdefault("CaloCellMakerTools_release" , [ 'ISF_CaloCellContainerFinalizerTool',
'ISF_FastHitConvertTool' ])
kwargs.setdefault("DoRandomFluctuations" , False )
+ kwargs.setdefault("useOneDShapeParametrisation" , False )
+ kwargs.setdefault("nHits" , 100 )
kwargs.setdefault("ParamsInputFilename" , ISF_FastCaloSimFlags.ParamsInputFilename())
+ kwargs.setdefault("Extrapolator" , 'TimedExtrapolator')
+ kwargs.setdefault("FastCaloSimCaloExtrapolation" , 'FastCaloSimCaloExtrapolation')
+ kwargs.setdefault("useFastCaloSimCaloExtrapolation" , True )
# register the FastCaloSim random number streams
from G4AtlasApps.SimFlags import simFlags
if not simFlags.RandomSeedList.checkForExistingSeed(ISF_FastCaloSimFlags.RandomStreamName()):
simFlags.RandomSeedList.addSeed( ISF_FastCaloSimFlags.RandomStreamName(), 98346412, 12461240 )
-
+
kwargs.setdefault("RandomStream" , ISF_FastCaloSimFlags.RandomStreamName())
kwargs.setdefault("RandomSvc" , simFlags.RandomSvc.get_Value() )
-
+
return CfgMgr.ISF__FastCaloSimSvcV2(name, **kwargs )
+
+
+def getFastCaloSimCaloExtrapolation(name="FastCaloSimCaloExtrapolation", **kwargs):
+ from ISF_FastCaloSimParametrization.ISF_FastCaloSimParametrizationConf import FastCaloSimCaloExtrapolation
+ kwargs.setdefault("CaloBoundaryR" , 1148.0 )
+ kwargs.setdefault("CaloBoundaryZ" , 3549.5 )
+ kwargs.setdefault("CaloMargin" , 100 )
+ kwargs.setdefault("Extrapolator" , "TimedExtrapolator" )
+ kwargs.setdefault("CaloSurfaceHelper" , "CaloSurfaceHelper" )
+ kwargs.setdefault("CaloGeometryHelper" , "FastCaloSimGeometryHelper" )
+ kwargs.setdefault("CaloEntrance" , "InDet::Containers::InnerDetector" )
+ return CfgMgr.FastCaloSimCaloExtrapolation(name, **kwargs)
+
+
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py
index 01a9b5a603ea..80b88a40c31f 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py
@@ -10,6 +10,7 @@ from AthenaCommon.CfgGetter import addTool, addService, addAlgorithm
addTool("ISF_FastCaloSimServices.AdditionalConfig.getNIMatEffUpdator", "ISF_NIMatEffUpdator")
addTool("ISF_FastCaloSimServices.AdditionalConfig.getNIPropagator", "ISF_NIPropagator")
addTool("ISF_FastCaloSimServices.AdditionalConfig.getNITimedExtrapolator", "ISF_NITimedExtrapolator")
+addTool("ISF_FastCaloSimServices.AdditionalConfig.getTimedExtrapolator", "TimedExtrapolator")
addTool("ISF_FastCaloSimServices.AdditionalConfig.getPunchThroughTool", "ISF_PunchThroughTool")
addTool("ISF_FastCaloSimServices.AdditionalConfig.getEmptyCellBuilderTool", "ISF_EmptyCellBuilderTool")
addTool("ISF_FastCaloSimServices.AdditionalConfig.getFastShowerCellBuilderTool", "ISF_FastShowerCellBuilderTool")
@@ -30,5 +31,4 @@ addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getFastCaloSim
addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getLegacyAFIIFastCaloSimSvc", "ISF_LegacyAFIIFastCaloSimSvc")
addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getFastHitConvAlgLegacyAFIIFastCaloSimSvc", "ISF_FastHitConvAlgLegacyAFIIFastCaloSimSvc")
addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getFastCaloSimSvcV2", "ISF_FastCaloSimSvcV2")
-
-
+addTool("ISF_FastCaloSimParametrization.ISF_FastCaloSimParametrizationConfig.getFastCaloSimCaloExtrapolation" , "FastCaloSimCaloExtrapolation" )
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx
index c2d9290a5a43..57d242f64f58 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx
@@ -9,10 +9,19 @@
// class header include
#include "FastCaloSimSvcV2.h"
+
// FastCaloSim includes
#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
#include "ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h"
+#include "ISF_FastCaloSimEvent/TFCSTruthState.h"
+#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
+//!
+#include "AtlasDetDescr/AtlasDetectorID.h"
+#include "IdDictParser/IdDictParser.h"
+#include "CaloIdentifier/LArEM_ID.h"
+#include "CaloIdentifier/TileID.h"
+//!
// StoreGate
#include "StoreGate/StoreGateSvc.h"
@@ -22,8 +31,8 @@
#include "CLHEP/Random/RandFlat.h"
#include "CaloEvent/CaloCellContainer.h"
-#include "ISF_FastCaloSimParametrization/CaloGeometryFromFile.h"
-#include "ISF_FastCaloSimParametrization/CaloDetDescrElement.h"
+#include "CaloDetDescr/CaloDetDescrElement.h"
+#include "CaloDetDescr/CaloDetDescrManager.h"
#include "TCanvas.h" // Just for testing
#include <fstream>
@@ -36,13 +45,18 @@ ISF::FastCaloSimSvcV2::FastCaloSimSvcV2(const std::string& name, ISvcLocator* sv
BaseSimulationSvc(name, svc),
m_rndGenSvc("AtRndmGenSvc", name)
{
- declareProperty("ParamsInputFilename" , m_paramsFilename);
- declareProperty("CaloCellsOutputName" , m_caloCellsOutputName) ;
- declareProperty("CaloCellMakerTools_setup" , m_caloCellMakerToolsSetup) ;
- declareProperty("CaloCellMakerTools_release" , m_caloCellMakerToolsRelease) ;
- declareProperty("DoRandomFluctuations" , m_doRandomFluctuations) ;
- declareProperty("RandomSvc" , m_rndGenSvc );
- declareProperty("RandomStream" , m_randomEngineName );
+ declareProperty("ParamsInputFilename" , m_paramsFilename);
+ declareProperty("CaloCellsOutputName" , m_caloCellsOutputName) ;
+ declareProperty("CaloCellMakerTools_setup" , m_caloCellMakerToolsSetup) ;
+ declareProperty("CaloCellMakerTools_release" , m_caloCellMakerToolsRelease) ;
+ declareProperty("DoRandomFluctuations" , m_doRandomFluctuations) ;
+ declareProperty("RandomSvc" , m_rndGenSvc );
+ declareProperty("RandomStream" , m_randomEngineName );
+ declareProperty("Extrapolator" , m_extrapolator );
+ declareProperty("FastCaloSimCaloExtrapolation" , m_FastCaloSimCaloExtrapolation );
+ declareProperty("useFastCaloSimCaloExtrapolation", m_useFastCaloSimCaloExtrapolation );
+ declareProperty("useOneDShapeParametrisation" , m_useOneDShapeParametrisation );
+ declareProperty("nHits" , m_nHits );
}
ISF::FastCaloSimSvcV2::~FastCaloSimSvcV2()
@@ -51,47 +65,93 @@ ISF::FastCaloSimSvcV2::~FastCaloSimSvcV2()
/** framework methods */
StatusCode ISF::FastCaloSimSvcV2::initialize()
{
- ATH_MSG_DEBUG(m_screenOutputPrefix << "Initializing ...");
+ ATH_MSG_INFO(m_screenOutputPrefix << "Initializing ...");
+
+ double wiggleLayer1[]={0.0110626,0.0242509,0.0398173,0.055761,0.0736173,0.0938847,0.115154,0.13639,0.157644,0.178934,0.200182,0.221473,0.242745,0.264019,0.285264,0.306527,0.327811,0.349119,0.370387,0.391668,0.412922,0.434208,0.45546,0.476732,0.498023,0.51931,0.540527,0.561799,0.583079,0.604358,0.625614,0.646864,0.668112,0.689351,0.710629,0.731894,0.75318,0.774426,0.795695,0.81699,0.838258,0.859528,0.880783,0.90202,0.922515,0.941276,0.958477,0.975062,0.988922,1};
+ double wiggleLayer2[]={0.0127507,0.0255775,0.0395137,0.0542644,0.0695555,0.0858206,0.102274,0.119653,0.137832,0.156777,0.176938,0.197727,0.217576,0.236615,0.256605,0.277766,0.2995,0.321951,0.344663,0.367903,0.392401,0.417473,0.443514,0.470867,0.498296,0.52573,0.553114,0.57921,0.604326,0.628822,0.652191,0.674853,0.697268,0.718983,0.739951,0.759866,0.778877,0.798762,0.819559,0.839789,0.858923,0.877327,0.894831,0.911693,0.92821,0.94391,0.959156,0.973593,0.986752,1};
+ double wiggleLayer3[]={0.0217932,0.0438502,0.0670992,0.091085,0.11651,0.143038,0.169524,0.196205,0.222944,0.249703,0.276629,0.303559,0.33034,0.356842,0.383579,0.410385,0.437272,0.464214,0.49118,0.518202,0.545454,0.572667,0.600037,0.627544,0.655072,0.6826,0.709824,0.733071,0.754764,0.775672,0.793834,0.810904,0.828219,0.844119,0.858339,0.871248,0.882485,0.894889,0.907955,0.920289,0.931136,0.941039,0.949844,0.957641,0.965787,0.97392,0.981706,0.988892,0.994527,1};
+
+ for(int i=0;i<50;i++)
+ {
+ m_wiggleLayer1[i]=wiggleLayer1[i];
+ m_wiggleLayer2[i]=wiggleLayer2[i];
+ m_wiggleLayer3[i]=wiggleLayer3[i];
+ }
+
+ cellcheck=0;
+
+ if(cellcheck)
+ {
+ detID=new AtlasDetectorID();
+ IdDictParser* parser = new IdDictParser;
+ IdDictMgr& idd = parser->parse ("IdDictParser/ATLAS_IDS.xml");
+ detID->initialize_from_dictionary(idd);
+
+ larID=new LArEM_ID();
+ larID->initialize_from_dictionary(idd);
+
+ tileID=new TileID();
+ tileID->set_do_neighbours(false);
+ tileID->initialize_from_dictionary(idd);
+
+ hecID=new LArHEC_ID();
+ hecID->initialize_from_dictionary(idd);
+
+ fcalID=new LArFCAL_ID();
+ fcalID->set_do_neighbours(false);
+ fcalID->initialize_from_dictionary(idd);
+
+ }
+
m_paramsFile = TFile::Open(m_paramsFilename.c_str());
- ATH_MSG_DEBUG("------> file = " << m_paramsFile);
-
- TFCSPCAEnergyParametrization* epara_pions=new TFCSPCAEnergyParametrization("pions","pions");
- TFCSPCAEnergyParametrization* epara_photons=new TFCSPCAEnergyParametrization("photons","photons");
+ ATH_MSG_INFO("------> file = " << m_paramsFile);
+ m_paramsFile_photons = TFile::Open("/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/secondPCA_photons_july17.root");
+ ATH_MSG_INFO("------> photon file = " << m_paramsFile_photons);
+ m_photonFile = TFile::Open("/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/photons_test_norm1.root");
+ ATH_MSG_INFO("------> file = " << m_photonFile);
+ m_elFile = TFile::Open("/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/InputDistribution_el_50.000000GeV_eta_0.200000_0.250000.root");
+ ATH_MSG_INFO("------> file = " << m_elFile);
+ m_pionFile = TFile::Open("/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/InputDistribution_pion_50.000000GeV_eta_0.200000_0.250000.root");
+ ATH_MSG_INFO("------> file = " << m_pionFile);
+
+ TFCSPCAEnergyParametrization* epara_pions =new TFCSPCAEnergyParametrization("pions","pions");
+ TFCSPCAEnergyParametrization* epara_photons =new TFCSPCAEnergyParametrization("photons","photons");
TFCSPCAEnergyParametrization* epara_electrons=new TFCSPCAEnergyParametrization("electrons","electrons");
-
+
m_eparas.push_back(epara_pions);
m_eparas.push_back(epara_photons);
m_eparas.push_back(epara_electrons);
-
- ATH_MSG_DEBUG("loading inputs pions");
+
+ ATH_MSG_INFO("loading inputs pions");
m_eparas[0]->loadInputs(m_paramsFile,"EnergyParams/pdgid_211/EN_50000/eta_0_20");
- ATH_MSG_DEBUG("loading inputs photons");
- m_eparas[1]->loadInputs(m_paramsFile,"EnergyParams/pdgid_22/EN_50000/eta_0_20");
- ATH_MSG_DEBUG("loading inputs electrons");
+ ATH_MSG_INFO("loading inputs photons");
+ m_eparas[1]->loadInputs(m_paramsFile_photons,"");
+ ATH_MSG_INFO("loading inputs electrons");
m_eparas[2]->loadInputs(m_paramsFile,"EnergyParams/pdgid_11/EN_65000/eta_0_20");
- ATH_MSG_DEBUG("loading electrons done!");
-
+ ATH_MSG_INFO("loading electrons done!");
+
/** commenting out close and delete: accessing the file in ::simulate */
- // file->Close();
- //delete file;
-
+ m_paramsFile->Close();
+ delete m_paramsFile;
+
ATH_CHECK(m_rndGenSvc.retrieve());
m_randomEngine = m_rndGenSvc->GetEngine( m_randomEngineName);
if(!m_randomEngine)
- {
- ATH_MSG_ERROR("Could not get random number engine from RandomNumberService. Abort.");
- return StatusCode::FAILURE;
- }
-
- m_caloGeo = new CaloGeometryFromFile();
+ {
+ ATH_MSG_ERROR("Could not get random number engine from RandomNumberService. Abort.");
+ return StatusCode::FAILURE;
+ }
+
+ const CaloDetDescrManager* calo_dd_man = CaloDetDescrManager::instance();
+ m_caloGeo = new CaloGeometryFromCaloDDM();
+ m_caloGeo->LoadGeometryFromCaloDDM(calo_dd_man);
TString path_to_fcal_geo_files = "/afs/cern.ch/atlas/groups/Simulation/FastCaloSimV2/";
- m_caloGeo->LoadGeometryFromFile("/afs/cern.ch/atlas/groups/Simulation/FastCaloSim/ATLAS-GEO-20-00-01.root", "ATLAS-GEO-20-00-01");
m_caloGeo->LoadFCalGeometryFromFiles(path_to_fcal_geo_files + "FCal1-electrodes.sorted.HV.09Nov2007.dat", path_to_fcal_geo_files + "FCal2-electrodes.sorted.HV.April2011.dat", path_to_fcal_geo_files + "FCal3-electrodes.sorted.HV.09Nov2007.dat");
-
+
for(unsigned int i=0;i<24;i++)
- m_rlayer.push_back(999);
-
+ m_rlayer.push_back(999);
+
m_rlayer[0]=1455.98;
m_rlayer[1]=1542.83;
m_rlayer[2]=1757.11;
@@ -117,391 +177,768 @@ StatusCode ISF::FastCaloSimSvcV2::initialize()
m_rlayer[22]=999;
m_rlayer[23]=999;
+ // Get TimedExtrapolator
+ if(!m_extrapolator.empty() && m_extrapolator.retrieve().isFailure())
+ return StatusCode::FAILURE;
+
+
+ // Get FastCaloSimCaloExtrapolation
+ if(m_FastCaloSimCaloExtrapolation.retrieve().isFailure())
+ {
+ ATH_MSG_ERROR("FastCaloSimCaloExtrapolation not found ");
+ return StatusCode::FAILURE;
+ }
+
+ rnd = new TRandom();
+
return StatusCode::SUCCESS;
}
/** framework methods */
StatusCode ISF::FastCaloSimSvcV2::finalize()
{
- ATH_MSG_DEBUG(m_screenOutputPrefix << "Finalizing ...");
+ ATH_MSG_INFO(m_screenOutputPrefix << "Finalizing ...");
return StatusCode::SUCCESS;
}
StatusCode ISF::FastCaloSimSvcV2::setupEvent()
{
- ATH_MSG_DEBUG(m_screenOutputPrefix << "setupEvent NEW EVENT!");
-
+ ATH_MSG_INFO(m_screenOutputPrefix << "setupEvent NEW EVENT!");
+
m_theContainer = new CaloCellContainer(SG::VIEW_ELEMENTS);
- // Only record if the container has not already been recorded
- if (!evtStore()->contains<CaloCellContainer>(m_caloCellsOutputName)) {
- if (evtStore()->record(m_theContainer, m_caloCellsOutputName).isFailure()) {
- ATH_MSG_FATAL( m_screenOutputPrefix << "cannot record CaloCellContainer " << m_caloCellsOutputName );
- return StatusCode::FAILURE;
- }
+ StatusCode sc = evtStore()->record(m_theContainer, m_caloCellsOutputName);
+ if (sc.isFailure())
+ {
+ ATH_MSG_FATAL( m_screenOutputPrefix << "cannot record CaloCellContainer " << m_caloCellsOutputName );
+ return StatusCode::FAILURE;
}
- ATH_CHECK( m_caloCellMakerToolsSetup.retrieve() );
+ CHECK( m_caloCellMakerToolsSetup.retrieve() );
ATH_MSG_DEBUG( "Successfully retrieve CaloCellMakerTools: " << m_caloCellMakerToolsSetup );
ToolHandleArray<ICaloCellMakerTool>::iterator itrTool = m_caloCellMakerToolsSetup.begin();
ToolHandleArray<ICaloCellMakerTool>::iterator endTool = m_caloCellMakerToolsSetup.end();
for (; itrTool != endTool; ++itrTool)
+ {
+ StatusCode sc = (*itrTool)->process(m_theContainer);
+ if (sc.isFailure())
{
- StatusCode sc = (*itrTool)->process(m_theContainer);
- if (sc.isFailure())
- {
- ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() );
- return StatusCode::FAILURE;
- }
+ ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() );
+ return StatusCode::FAILURE;
}
+ }
return StatusCode::SUCCESS;
}
StatusCode ISF::FastCaloSimSvcV2::releaseEvent()
{
-
- ATH_MSG_DEBUG(m_screenOutputPrefix << "release Event");
-
- ATH_CHECK( m_caloCellMakerToolsRelease.retrieve() );
- //run release tools in a loop
- ToolHandleArray<ICaloCellMakerTool>::iterator itrTool = m_caloCellMakerToolsRelease.begin();
- ToolHandleArray<ICaloCellMakerTool>::iterator endTool = m_caloCellMakerToolsRelease.end();
- for (; itrTool != endTool; ++itrTool)
- {
- ATH_MSG_DEBUG( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
-
- StatusCode sc = (*itrTool)->process(m_theContainer);
-
- if (sc.isFailure())
- {
- ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() );
- }
- }
-
- return StatusCode::SUCCESS;
-
+
+ ATH_MSG_INFO(m_screenOutputPrefix << "release Event");
+
+ CHECK( m_caloCellMakerToolsRelease.retrieve() );
+
+ //run release tools in a loop
+ ToolHandleArray<ICaloCellMakerTool>::iterator itrTool = m_caloCellMakerToolsRelease.begin();
+ ToolHandleArray<ICaloCellMakerTool>::iterator endTool = m_caloCellMakerToolsRelease.end();
+ for (; itrTool != endTool; ++itrTool)
+ {
+ ATH_MSG_INFO( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
+
+ StatusCode sc = (*itrTool)->process(m_theContainer);
+
+ if (sc.isFailure())
+ {
+ ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() );
+ }
+ }
+
+ return StatusCode::SUCCESS;
+
}
/** Simulation Call */
StatusCode ISF::FastCaloSimSvcV2::simulate(const ISF::ISFParticle& isfp)
{
- ATH_MSG_DEBUG("NEW PARTICLE! FastCaloSimSvcV2 called with ISFParticle: " << isfp);
-
- int pdgid = fabs(isfp.pdgCode());
- int barcode=isfp.barcode();
-
- if(barcode!=10001)
- {
- ATH_MSG_DEBUG("ISF particle barcode is barcode "<<barcode<<". Go to next Particle.");
- return StatusCode::SUCCESS;
- }
-
- Amg::Vector3D particle_position = isfp.position();
- double eta(0.), phi(0.);
- eta = particle_position.eta();
- phi = particle_position.phi();
- if(abs(eta) > 0.3 || abs(eta) < 0.15) //somewhat enlarged to not cut off too many particles
- {
- ATH_MSG_DEBUG("ISF particle is out of eta range: "<<eta<<". Go to next Particle.");
- return StatusCode::SUCCESS;
- }
- if(pdgid==22 || pdgid==211 || pdgid==11)
- {
-
- int index_epara=0;
- if(pdgid==22) index_epara=1;
- if(pdgid==11) index_epara=2;
-
- /** get the number of pca bins from the epara: **/
-
- int n_pcabins = m_eparas[index_epara]->n_pcabins();
-
- //determine the pca bin randomly:
- int pcabin = 1;
- float uniform=CLHEP::RandFlat::shoot(m_randomEngine);
-
- for (int n = 0; n < n_pcabins; n++)
- {
- if (uniform > n * 1.0 / (double)n_pcabins && uniform < (n + 1.) * 1.0 / (double)n_pcabins)
- pcabin = n + 1;
- }
-
- ATH_MSG_DEBUG("pca bin "<<pcabin);
-
- //-----ENERGY:-----
-
- TFCSSimulationState simulstate;
- simulstate.set_Ebin(pcabin);
-
- m_eparas[index_epara]->simulate(simulstate, nullptr, nullptr);
-
- ATH_MSG_DEBUG("Energy returned: " << simulstate.E());
- ATH_MSG_DEBUG("Energy fraction for layer: ");
- for (int s = 0; s < CaloCell_ID_FCS::MaxSample; s++)
- ATH_MSG_VERBOSE(" Sampling " << s << " energy " << simulstate.E(s));
-
- //-----SHAPE:-----
-
- if (m_doRandomFluctuations == true && pdgid == 211)
- {
- ATH_MSG_DEBUG("\n\n*******************************************\n"
- << "RANDOM FLUCTUATIONS ARE ON"
- << "\n*******************************************\n\n");
- }
-
- else {
-
- ATH_MSG_DEBUG("\n\n*******************************************\n"
- << "RANDOM FLUCTUATIONS ARE OFF"
- << "\n*******************************************\n\n");
- }
-
- /** get the relevant calo layers */
- IntArray *m_layers = m_eparas[index_epara]->get_layers();
-
- /** get the appropriate input histogram */
- std::string inputHisto = "";
-
- for ( int ilayer = 0; ilayer < m_layers->GetSize(); ilayer++)
- {
-
- /** get the layer ID and associated energy */
- int layer = m_layers->GetAt(ilayer);
- double layerE = simulstate.E(layer);
- int nHits=-1;
-
- ATH_MSG_DEBUG("NOW RUNNING ON LAYER "<<layer);
-
- inputHisto = "ShapeParams/pdgid_" + std::to_string(pdgid) + "/EN_50000/eta_0_20/hEnergyDensity_layer" + std::to_string(layer) + "_pca" + std::to_string(pcabin);
- m_histEnergyDensity2D = (TH2F*)m_paramsFile->Get(inputHisto.c_str());
-
- if(!m_paramsFile)
- {
- ATH_MSG_FATAL("No paramFile!");
- return StatusCode::FAILURE;
- }
- if(!m_histEnergyDensity2D)
- {
- ATH_MSG_FATAL("Problem. No histogram!! :-(");
- return StatusCode::FAILURE;
- }
-
- /** set number of hits due to random fluctuations if particle is a pion and flag is set */
- if (m_doRandomFluctuations == true && pdgid == 211) nHits = 10;
-
- /** otherwise get the number of hits to account for the resolution of the layer */
- else nHits = hitCalc(layerE, layer, pdgid);
-
- ATH_MSG_DEBUG("number of HITS = "<<nHits);
-
- /** simulate random hits from the input histogram */
- for (int i = 0; i < nHits; i++)
- {
- //std::cout<<"Hit nr "<<i<<std::endl;
-
- double r, alpha;
- double delta_phi,delta_phi_mm;
-
- const CaloGeoDetDescrElement* mcell = 0;
- m_histEnergyDensity2D->GetRandom2(alpha, r);
-
- //std::cout<<"got a hit positon from the histogram!"<<" r "<<r<<" alpha "<<alpha<<std::endl;
-
- //double r_layer=m_rlayers[ilayer*n_pcabins+pcabin-1];
- double r_layer=m_rlayer[layer];
- double hit_eta=findHitEta(alpha,r,r_layer,particle_position.z(),eta);
-
- delta_phi_mm = r * sin(alpha);
- delta_phi=delta_phi_mm/r_layer;
- double hit_phi = phi+delta_phi;
-
- if (layer < 21){
- double wiggle = 0.0;
- if (layer < 4 && layer > 0){
- wiggle = doWiggle(layer);
- //std::cout << "wiggle is " << wiggle << std::endl;
- }
- mcell = m_caloGeo->getDDE(layer, hit_eta, hit_phi-wiggle);
- }
-
- else
- {
-
- double theta=2.*atan(exp(-hit_eta));
- double rT=particle_position.z()*tan(theta);
- double hit_x=rT*cos(hit_phi);
- double hit_y=rT*sin(hit_phi);
-
- mcell = m_caloGeo->getFCalDDE(layer, hit_x, hit_y, particle_position.z());
- }
-
- if (!mcell) continue;
-
-
- CaloCell* theCell = (CaloCell*)m_theContainer->findCell(mcell->calo_hash());
-
- if(i<10)
- {
- ATH_MSG_DEBUG("Hit nr "<<i<<" eta: " << hit_eta << " phi: " << hit_phi << " Particle eta: " << eta << " phi: " << phi << " delta_eta: " << hit_eta - eta << " delta_phi: " << hit_phi - phi);
- ATH_MSG_DEBUG(" Cell from CaloGeometry: eta: " << mcell->eta() << " phi: " << mcell->phi() << " |CELL_eta - HIT_eta| " << abs(mcell->eta() - hit_eta) << " |CELL_phi - HIT_phi| " << abs(mcell->phi() - hit_phi));
- ATH_MSG_DEBUG(" energy input into cell: "<<layerE / nHits);
- }
-
- theCell->addEnergy(layerE / nHits);
-
- }
- }
-
- } //pdgid must be 211 or 22 or 11
- else
- {
- ATH_MSG_DEBUG("Oh no! ISF particle has pdgid "<<pdgid<<" . That's not supported yet :(");
- }
-
- return StatusCode::SUCCESS;
-
+
+ ATH_MSG_INFO("NEW PARTICLE! FastCaloSimSvcV2 called with ISFParticle: " << isfp);
+
+ int pdgid = fabs(isfp.pdgCode());
+ int barcode=isfp.barcode();
+
+ if(barcode!=10001)
+ {
+ ATH_MSG_INFO("ISF particle barcode is barcode "<<barcode<<". Go to next Particle.");
+ return StatusCode::SUCCESS;
+ }
+
+ particle_position = isfp.position();
+ eta = 0;
+ phi = 0;
+ eta_isfp = 0;
+ phi_isfp = 0;
+
+ eta_isfp = particle_position.eta(); // isfp eta and phi, in case we need them
+ phi_isfp = particle_position.phi();
+ if(abs(eta_isfp) > 0.3 || abs(eta_isfp) < 0.15) //somewhat enlarged to not cut off too many particles
+ {
+ ATH_MSG_INFO("ISF particle is out of eta range: "<<eta_isfp<<". Go to next Particle.");
+ return StatusCode::SUCCESS;
+ }
+
+ if(!(pdgid==22 || pdgid==211 || pdgid==11))
+ {
+ ATH_MSG_INFO("Oh no! ISF particle has pdgid "<<pdgid<<" . That's not supported yet :(");
+ return StatusCode::SUCCESS;
+ }
+
+ //Get extrapolation result
+
+ //commenting out next line, because it fails for some events (where the truth is not saved, this apparanetly can happen)
+ //HepMC::GenParticle* isfpTruth = isfp.getTruthBinding()->getTruthParticle();
+ //will use the ISFP momentum instead:
+ TFCSTruthState truth(isfp.momentum().x(),isfp.momentum().y(),isfp.momentum().z(),sqrt(pow(isfp.ekin(),2)+pow(isfp.mass(),2)),isfp.pdgCode());
+ TFCSExtrapolationState result;
+ m_FastCaloSimCaloExtrapolation->extrapolate(result,&truth);
+
+ int index_epara=0;
+ if(pdgid==22) index_epara=1;
+ if(pdgid==11) index_epara=2;
+
+ /** get the number of pca bins from the epara: **/
+ int n_pcabins = m_eparas[index_epara]->n_pcabins();
+
+ //determine the pca bin randomly:
+ int pcabin = 1;
+ float uniform=CLHEP::RandFlat::shoot(m_randomEngine);
+
+ for (int n = 0; n < n_pcabins; n++)
+ {
+ if(uniform > n * 1.0 / (double)n_pcabins && uniform < (n + 1.) * 1.0 / (double)n_pcabins)
+ pcabin = n + 1;
+ }
+ ATH_MSG_INFO("pca bin "<<pcabin);
+
+ //-----ENERGY:-----
+ simulstate.set_Ebin(pcabin);
+ m_eparas[index_epara]->simulate(simulstate, nullptr, nullptr);
+
+ ATH_MSG_INFO("Energy returned: " << simulstate.E());
+ ATH_MSG_INFO("Energy fraction for layer: ");
+ for (int s = 0; s < CaloCell_ID_FCS::MaxSample; s++)
+ ATH_MSG_INFO(" Sampling " << s << " energy " << simulstate.E(s));
+
+ //-----SHAPE:-----
+
+ if(m_doRandomFluctuations == true && pdgid == 211)
+ {
+ ATH_MSG_INFO("\n\n*******************************************\n"
+ << "RANDOM FLUCTUATIONS ARE ON"
+ << "\n*******************************************\n\n");
+ }
+ else
+ {
+ ATH_MSG_INFO("\n\n*******************************************\n"
+ << "RANDOM FLUCTUATIONS ARE OFF"
+ << "\n*******************************************\n\n");
+ }
+
+ if (m_useOneDShapeParametrisation == true )
+ {
+ ATH_MSG_INFO("\n\n*******************************************\n"
+ << "USING 1D SHAPE PARAMETRISATION\n"
+ << "NHITS=" << m_nHits
+ << "\n*******************************************\n\n");
+ }
+ else
+ {
+ ATH_MSG_INFO("\n\n*******************************************\n"
+ << "USING 1D SHAPE PARAMETRISATION"
+ << "\n*******************************************\n\n");
+ }
+
+ /** get the relevant calo layers */
+ m_layers = m_eparas[index_epara]->get_layers();
+
+ /** get the appropriate input histogram */
+ std::string inputHistoName = "";
+
+ esum=0.0;
+ nCells=0;
+ nCells_Tile=0;
+
+ for ( ilayer = 0; ilayer < m_layers->GetSize(); ilayer++)
+ {
+ /** get the layer ID and associated energy */
+ layer = m_layers->GetAt(ilayer);
+ layerE = simulstate.E(layer);
+
+ // get eta and phi from extrapolator or isfp
+ if(m_useFastCaloSimCaloExtrapolation)
+ {
+ ATH_MSG_INFO("Using FastCaloSimCaloExtrapolation result for eta, phi, r and z!");
+ eta = result.eta(layer, 0);
+ phi = result.phi(layer, 0);
+ r_layer = result.r(layer,0);
+ z_particle = result.z(layer,0);
+ }
+ if(!m_useFastCaloSimCaloExtrapolation)
+ {
+ ATH_MSG_INFO("Using isfp for eta, phi, r and z!");
+ eta = eta_isfp;
+ phi = phi_isfp;
+ r_layer = m_rlayer[layer];
+ z_particle = particle_position.z();
+ }
+
+ ATH_MSG_INFO("NOW RUNNING ON LAYER "<<layer);
+
+ inputHistoName = "hEnergy_layer"+std::to_string(layer)+"_pca"+std::to_string(pcabin);
+
+ if(pdgid==11) //el
+ m_histEnergyDensity2D = (TH2F*)m_elFile->Get(inputHistoName.c_str());
+ if(pdgid==211) //pion
+ m_histEnergyDensity2D = (TH2F*)m_pionFile->Get(inputHistoName.c_str());
+ if(pdgid==22) //photons
+ m_histEnergyDensity2D = (TH2F*)m_photonFile->Get(inputHistoName.c_str());
+
+ if(!m_photonFile)
+ {
+ ATH_MSG_FATAL("No photonFile!");
+ return StatusCode::FAILURE;
+ }
+ if(!m_elFile)
+ {
+ ATH_MSG_FATAL("No elFile!");
+ return StatusCode::FAILURE;
+ }
+ if(!m_pionFile)
+ {
+ ATH_MSG_FATAL("No pionFile!");
+ return StatusCode::FAILURE;
+ }
+ if(!m_paramsFile)
+ {
+ ATH_MSG_FATAL("No paramFile!");
+ return StatusCode::FAILURE;
+ }
+ if(!m_histEnergyDensity2D)
+ {
+ ATH_MSG_FATAL("Problem. No histogram!! :-(");
+ return StatusCode::FAILURE;
+ }
+
+ if(m_useOneDShapeParametrisation)
+ {
+ if (m_doRandomFluctuations == true && pdgid == 211) m_nHits = 5; //to be optimised
+ }
+ if(!m_useOneDShapeParametrisation)
+ {
+ /** set number of hits due to random fluctuations if particle is a pion and flag is set */
+ if (m_doRandomFluctuations == true && pdgid == 211) m_nHits = 10;
+ /** otherwise get the number of hits to account for the resolution of the layer */
+ else m_nHits = hitCalc(layerE, layer, pcabin, pdgid);
+ }
+
+ ATH_MSG_INFO("number of HITS = "<<m_nHits);
+
+ R=sqrt(r_layer*r_layer + z_particle*z_particle);
+ findEtaMMRange(R,eta,deltaEtaMMmin,deltaEtaMMmax);
+ deltaEtaMMmin*=0.999;
+ deltaEtaMMmax*=0.999;
+
+ /*
+ std::cout<<"Range for Delta eta[mm]: " << "[" << deltaEtaMMmin << "," << deltaEtaMMmax << "]" << std::endl;
+ std::cout << "Particle position: R: " << R << " eta: " << eta << std::endl;
+ std::cout<<"layer "<<layer<<" number of HITS = "<<m_nHits<<std::endl;
+ */
+
+ if(m_useOneDShapeParametrisation)
+ {
+ TH1D* proj = m_histEnergyDensity2D->ProjectionY("proj",0,-1,"e");
+
+ for (int i = 1; i <= proj->GetNbinsX(); ++i)
+ {
+ double mean, err, energyDensityInBin, energyInBin;
+ mean = proj->GetBinContent(i);
+ err = proj->GetBinError(i);
+ //std::cout << "Bin: " << i << " mean: " << mean << " err: " << err << std::endl;
+
+ energyDensityInBin = rnd->Gaus(mean, err);
+ energyInBin = energyDensityInBin * layerE;
+
+ double minR = proj->GetBinLowEdge(i);
+ double maxR = proj->GetBinLowEdge(i) + proj->GetBinWidth(i);
+
+ //std::cout << "Energy in layer: " << layerE << " energyDensityInBin: " << energyDensityInBin << " energyInBin: " << energyInBin << std::endl;
+ LoopOverHits(energyInBin, minR, maxR);
+ if(cellcheck) TestCell();
+ }
+ }
+ if(!m_useOneDShapeParametrisation)
+ {
+ LoopOverHits(layerE);
+ if(cellcheck) TestCell();
+ }
+
+ } //loop layer
+
+ if(cellcheck) std::cout<<"ECHECK etot "<<esum<<" epara "<<simulstate.E()<<" nCells "<<nCells<<" nCells_Tile "<<nCells_Tile<<std::endl;
+
+ delete m_histEnergyDensity2D;
+ return StatusCode::SUCCESS;
+
}
-
+void ISF::FastCaloSimSvcV2::LoopOverHits(double totalEnergy, double minR, double maxR)
+{
+
+ for (int hit = 0; hit < m_nHits; hit++)
+ {
+ double alpha, r;
+ double delta_eta_mm;
+
+ const CaloDetDescrElement* mcell = 0;
+ do
+ {
+ if(m_useOneDShapeParametrisation)
+ {
+ r = rnd->Uniform(minR, maxR);
+ alpha = rnd->Uniform(2*TMath::Pi());
+ }
+ if(!m_useOneDShapeParametrisation)
+ {
+ m_histEnergyDensity2D->GetRandom2(alpha, r);
+ }
+ delta_eta_mm=r * cos(alpha);
+ }
+ while(delta_eta_mm<deltaEtaMMmin || delta_eta_mm>deltaEtaMMmax);
+
+ double delta_phi_mm = r * sin(alpha);
+ //std::cout<<"got a hit positon from the histogram!"<<" r "<<r<<" alpha "<<alpha<<" r_layer "<<r_layer<<" z_particle "<<z_particle<<" eta "<<eta<<std::endl;
+ //double r_layer=m_rlayers[ilayer*n_pcabins+pcabin-1];
+
+ double hit_eta=findHitEta(delta_eta_mm,R,eta);
+ //std::cout<<"hit_eta "<<hit_eta<<std::endl;
+
+ double delta_phi=delta_phi_mm/r_layer;
+ double hit_phi = phi+delta_phi;
+ double hit_phi_shifted=hit_phi;
+
+ double wiggle = 0.0;
+ if(layer < 21)
+ {
+ if(layer < 4 && layer > 0)
+ {
+ wiggle = doWiggle(layer);
+ //std::cout << "wiggle is " << wiggle << std::endl;
+ }
+ hit_phi_shifted=hit_phi-wiggle;
+ hit_phi_shifted=TVector2::Phi_mpi_pi(hit_phi_shifted);
+ mcell = m_caloGeo->getDDE(layer, hit_eta, hit_phi_shifted);
+ }
+ if(layer >= 21) //FCAL
+ {
+ double theta=2.*atan(exp(-hit_eta));
+ double rT=z_particle*tan(theta);
+ double hit_x=rT*cos(hit_phi);
+ double hit_y=rT*sin(hit_phi);
+ mcell = m_caloGeo->getFCalDDE(layer, hit_x, hit_y, z_particle);
+ }
+
+ if (!mcell) continue;
+
+ CaloCell* theCell = (CaloCell*)m_theContainer->findCell(mcell->calo_hash());
+ theCell->addEnergy(totalEnergy / m_nHits);
+
+ /*
+ if(hit<10)
+ {
+ ATH_MSG_INFO("Hit nr "<<i<<" eta: " << hit_eta << " phi: " << hit_phi << " Particle eta: " << eta << " phi: " << phi << " delta_eta: " << hit_eta - eta << " delta_phi: " << hit_phi - phi);
+ ATH_MSG_INFO(" Cell from CaloGeometry: eta: " << mcell->eta() << " phi: " << mcell->phi() << " |CELL_eta - HIT_eta| " << abs(mcell->eta() - hit_eta) << " |CELL_phi - HIT_phi| " << abs(mcell->phi() - hit_phi));
+ ATH_MSG_INFO(" energy density in layer: "<<energyDensityInBin);
+ ATH_MSG_INFO(" energy input into cell: "<<energyInBin / m_nHits);
+ }
+ */
+
+ if(cellcheck){
+ const CaloDetDescrElement* DDE = theCell->caloDDE();
+
+ // Comparison CaloDetDescrElement vs CaloGeometry element
+ if(abs(DDE->eta() - mcell->eta())>0.0001 || abs(DDE->phi() - mcell->phi())>0.0001 ){
+ ATH_MSG_ERROR("Error in CaloGeometry: Cell with same ID has different position.");
+ ATH_MSG_ERROR("CaloCell: " << "eta: " << DDE->eta() << " phi: " << DDE->phi() << " ID: " << DDE->identify());
+ ATH_MSG_ERROR("CaloGeometry: " << "eta: " << mcell->eta() << " phi: " << mcell->phi() << " ID: " << mcell->identify());
+
+ }
+
+ // Hit-to-cell assignment check
+
+ //double delta_phi_check=abs(std::fmod(DDE->phi() - hit_phi_shifted,TMath::Pi()*2.));
+ double delta_phi_check=abs(TVector2::Phi_mpi_pi(DDE->phi() - hit_phi_shifted));
+ //if(delta_phi_check>TMath::Pi())delta_phi_check=TMath::Pi()*2.-delta_phi_check;
+ if(abs(DDE->eta() - hit_eta) > DDE->deta()/2*1.01 || abs(delta_phi_check) > DDE->dphi()/2*1.01) {
+ ATH_MSG_ERROR("Error in hit-to-cell assignment: Hit is too far from the cell.");
+ ATH_MSG_ERROR(" Cell from CaloGeometry: eta: " << DDE->eta() << " deta: " << DDE->deta() << " phi: " << DDE->phi() << " dphi: " << DDE->dphi() << " |CELL_eta - HIT_eta|/(0.5*deta) " << abs(DDE->eta() - hit_eta)*2./ DDE->deta()
+ << " HIT_eta " << hit_eta << " |CELL_phi - HIT_phi|/(0.5*dphi) " << abs(delta_phi_check)*2./DDE->dphi() << " HIT_phi: " << hit_phi_shifted );
+ }
+
+
+ if(hit<10)
+ {
+ ATH_MSG_INFO("Hit nr "<<hit<<" eta: " << hit_eta << " phi: " << hit_phi << " Particle eta: " << eta << " phi: " << phi << " delta_eta: " << hit_eta - eta << " delta_phi: " << hit_phi - phi);
+ ATH_MSG_INFO(" Cell from CaloGeometry: eta: " << mcell->eta() << " phi: " << mcell->phi() << " |CELL_eta - HIT_eta|/(0.5*CELL_deta) " << abs(mcell->eta() - hit_eta)*2/mcell->deta() << " |CELL_phi - HIT_phi|/(0.5*CELL_dphi) " << abs(mcell->phi() - hit_phi_shifted)*2./mcell->dphi());
+ ATH_MSG_INFO(" energy input into cell: "<<layerE / m_nHits);
+ }
+} //if cellcheck
+
+ } //for hit
+
+} //LoopOverHits
+
+void ISF::FastCaloSimSvcV2::TestCell()
+{
+ esum=0.0;
+ nCells=0;
+ nCells_Tile=0;
+
+ double cut=0.1;
+
+ CaloCellContainer::iterator it =m_theContainer->begin();
+ CaloCellContainer::iterator it_e=m_theContainer->end();
+
+ for(;it!=it_e;++it)
+ {
+ CaloCell* theCell=(*it);
+ //const unsigned int hash_id=theCell->caloDDE()->calo_hash();
+ const Identifier cell_id=theCell->ID();
+ esum+=theCell->energy();
+ if(theCell->energy()>cut)
+ {
+ int isLAr=0; int isTile=0; int isHEC=0; int isFCAL=0;
+
+ if(detID->is_lar_em(cell_id)==1) isLAr=1;
+ if(detID->is_lar_hec(cell_id)==1) isHEC=1;
+ if(detID->is_lar_fcal(cell_id)==1) isFCAL=1;
+ if(detID->is_tile(cell_id)==1) isTile=1;
+
+ std::cout<<"cell_id "<<cell_id<<" energy "<<theCell->energy()<<" eta "<<theCell->eta()<<" phi "<<phi;
+ if(isLAr) { std::cout<<" LArEM sampling "<<larID->sampling(cell_id); }
+ if(isHEC) { std::cout<<" HEC sampling "<<hecID->sampling(cell_id); }
+ if(isFCAL) { std::cout<<" FCAL module "<<fcalID->module(cell_id); }
+ if(isTile) { std::cout<<" TILE sampling "<<tileID->sampling(cell_id); nCells_Tile++; }
+
+ std::cout<<" "<<std::endl;
+ nCells++;
+ }
+ }
+
+ double sum_epara=0.0;
+ for(int a=0;a<=ilayer;a++) sum_epara+=simulstate.E(m_layers->GetAt(a));
+
+ std::cout<<"ECHECK layer "<<layer<<" esum "<<esum<<" epara "<<sum_epara<<" (this layer: "<<simulstate.E(layer)<<") nCells "<<nCells<<" nCells_Tile "<<nCells_Tile<<std::endl;
+}
//-----HITS:-----
-int ISF::FastCaloSimSvcV2::hitCalc(double energy, int calolayer, int pdgid){
-
- int hits=1;
-
- /** Calculate estimates for cell energy fluctuations. Should be updated with better numbers */
-
- if(energy>0.0001)
- {
-
- double stochastic=0.1;
- double constant=0.0;
-
- if(pdgid == 211){
-
- if(calolayer >= 0 && calolayer<= 7){
- stochastic=0.101; //LAr 10.1%/sqrt(E)
- constant=0.002;
- }
-
- else if(calolayer >= 8 && calolayer <= 11){
- stochastic=0.706; //HadEC 70.6%/sqrt(E) for pions
- constant=0.058;
- }
-
- else if(calolayer >= 12 && calolayer <= 20){
- stochastic=0.564; //TileCal 56.4%/sqrt(E)
- constant=0.055;
- }
-
- else if(calolayer >= 21 && calolayer <= 23){
- stochastic=0.942; //FCAL 94.2%/sqrt(E) for pions
- constant=0.075;
- }
- }
-
- else {
- if(calolayer >= 0 && calolayer<= 7){
- stochastic=0.101; //LAr 10.1%/sqrt(E)
- constant=0.002;
- }
-
- else if(calolayer >= 8 && calolayer <= 11){
- stochastic=0.214; //HadEC 21.4%/sqrt(E)
- constant=0.0;
- }
-
- else if(calolayer >= 12 && calolayer <= 20){
- stochastic=0.564; //TileCal 56.4%/sqrt(E)
- constant=0.055;
- }
-
- else if(calolayer >= 21 && calolayer <= 23){
- stochastic=0.285; //FCAL 28.5%/sqrt(E)
- constant=0.035;
- }
- }
-
- hits = 1.0 / ((stochastic/sqrt(energy/1000.0))*(stochastic/sqrt(energy/1000.0)) + constant*constant);
-
- }
- else
- ATH_MSG_DEBUG("Something is wrong. Layer "<<calolayer<<" has very little energy "<<energy);
+int ISF::FastCaloSimSvcV2::hitCalc(double energy, int calolayer,int pcabin, int pdgid)
+{
+
+ int hits=1;
+
+ std::vector<std::vector<int> > hitVectorPions(24);
+ //for(int i=0;i<24;i++)hitVector[i].resize(5);
+ hitVectorPions[0]={105,62,30,9,4};
+ hitVectorPions[1]={569,283,117,31,7};
+ hitVectorPions[2]={2320,1859,1189,380,35};
+ hitVectorPions[3]={103,128,136,50,3};
+ hitVectorPions[12]={18,28,37,49,36};
+ hitVectorPions[13]={6,14,26,44,70};
+ hitVectorPions[14]={0,0,4,0,8};
+
+ std::vector<std::vector<int> > hitVectorPhotons(24);
+ hitVectorPhotons[0]={8,10,5,5,8};
+ hitVectorPhotons[1]={5,5,5,5,8};
+ hitVectorPhotons[2]={25,55,60,150,150};
+ hitVectorPhotons[3]={50,50,50,50,50};
+ hitVectorPhotons[12]={30,40,50,25,30};
+ hitVectorPhotons[13]={0,0,0,0,0};
+ hitVectorPhotons[14]={0,0,0,0,0};
+ /** Calculate estimates for cell energy fluctuations. Should be updated with better numbers */
+
+ if(energy>0.0001)
+ {
+
+ double stochastic=0.1;
+ double constant=0.0;
+
+ if(pdgid == 211)
+ {
+ return hitVectorPions[calolayer][pcabin];
+
+ /*
+ if(calolayer >= 0 && calolayer<= 7)
+ {
+ stochastic=0.101; //LAr 10.1%/sqrt(E)
+ constant=0.002;
+ }
+ if(calolayer >= 8 && calolayer <= 11)
+ {
+ stochastic=0.706; //HadEC 70.6%/sqrt(E) for pions
+ constant=0.058;
+ }
+ if(calolayer >= 12 && calolayer <= 20)
+ {
+ stochastic=0.564; //TileCal 56.4%/sqrt(E)
+ constant=0.055;
+ }
+ if(calolayer >= 21 && calolayer <= 23)
+ {
+ stochastic=0.942; //FCAL 94.2%/sqrt(E) for pions
+ constant=0.075;
+ }
+ */
+ }
- return hits;
+ if(pdgid != 211)
+ {
+ if(calolayer >= 0 && calolayer<= 7)
+ {
+ stochastic=0.101; //LAr 10.1%/sqrt(E)
+ constant=0.002;
+ }
+
+ if(calolayer >= 8 && calolayer <= 11)
+ {
+ stochastic=0.214; //HadEC 21.4%/sqrt(E)
+ constant=0.0;
+ }
+
+ if(calolayer >= 12 && calolayer <= 20)
+ {
+ stochastic=0.564; //TileCal 56.4%/sqrt(E)
+ constant=0.055;
+ }
+
+ if(calolayer >= 21 && calolayer <= 23)
+ {
+ stochastic=0.285; //FCAL 28.5%/sqrt(E)
+ constant=0.035;
+ }
+
+ } //if pdgid!=211
+
+ hits = 1.0 / ((stochastic/sqrt(energy/1000.0))*(stochastic/sqrt(energy/1000.0)) + constant*constant);
+
+ } //if pass energy
+
+ else
+ ATH_MSG_INFO("Something is wrong. Layer "<<calolayer<<" has very little energy "<<energy);
+
+ return hits;
}
-double ISF::FastCaloSimSvcV2::findHitEta(const double alpha,const double r, const double r_layer,const double z_particle,const double eta_particle){
- const double delta_eta_mm = r * cos(alpha);
- const double R=sqrt(r_layer*r_layer + z_particle*z_particle);
- double x=eta_particle; // approximation of hit_eta, starting point for iterations
+double ISF::FastCaloSimSvcV2::findHitEta(const double delta_eta_mm,const double R,const double eta_particle)
+{
+ double x=eta_particle; // approximation of hit_eta, starting point for iterations
+
+ //std::cout<<"in findHitEta: delta_eta_mm "<<delta_eta_mm<<" R "<<R<<" x "<<x<<std::endl;
+
double a,b;
double c,d;
double delta;
- double e;
-
- const double epsilon=0.0001;
- do{
- c=exp(-x);
- d=1./(1.+(c*c));
- e=(2.*R)*(c*d);
+ double sech,tanh;
+
+ int count=0;
+ const double epsilon=0.001;
+ do
+ {
+ c=exp(x);
+ d=1./c;
+ sech=2./(c+d);
+ tanh=(c-d)/2.*sech;
delta=x-eta_particle;
- a = e*delta - delta_eta_mm; // transformation
- b = e*(1. - (delta*d)); // first derivative of the transformation
+ a = transformationEtaToMM(R,delta,sech) - delta_eta_mm;// transformation
+ b = transformationEtaToMMFirstDerivative(R,delta,sech,tanh);// first derivative of the transformation
x = x - a/b;
- } while(abs(a) > epsilon);
- if(x!=x) ATH_MSG_ERROR("Error: Hit eta not defined.");
+ count++;
+
+ //std::cout<<"count "<<count<<" a "<<a<<" b "<<b<<" c "<<c<<" d "<<d<<" delta "<<delta<<" x "<<x<<std::endl;
+ }
+ while((abs(a) > epsilon) && count<1000);
+ if(x!=x || count>=1000) ATH_MSG_ERROR("Error: Hit eta not found.");
//std::cout << "hit_eta: " << x << " error: " << a/b << std::endl;
return x;
}
-double ISF::FastCaloSimSvcV2::doWiggle(int layer)
+double ISF::FastCaloSimSvcV2::findStartingPoint(const double R,const double eta_particle,bool forMaximum)
{
-
- double wiggle = 0.0;
- double cell_dphi = 0.0;
-
- //Define cell size in phi
- if (layer == 1)
- cell_dphi = 0.0981748;
- else if (layer == 2 || layer == 3)
- cell_dphi = 0.0245437;
- else{
- ATH_MSG_ERROR("I am inside the wiggle calculation, but the layer is " << layer << "!");
- return 0.0;
- }
-
- //Set random numbers
- double searchRand = CLHEP::RandFlat::shoot(m_randomEngine);
-
- //Now for layer dependant approach
- if (layer == 1){
- double m_wiggleLayer1[]={0.0110626,0.0242509,0.0398173,0.055761,0.0736173,0.0938847,0.115154,0.13639,0.157644,0.178934,0.200182,0.221473,0.242745,0.264019,0.285264,0.306527,0.327811,0.349119,0.370387,0.391668,0.412922,0.434208,0.45546,0.476732,0.498023,0.51931,0.540527,0.561799,0.583079,0.604358,0.625614,0.646864,0.668112,0.689351,0.710629,0.731894,0.75318,0.774426,0.795695,0.81699,0.838258,0.859528,0.880783,0.90202,0.922515,0.941276,0.958477,0.975062,0.988922,1};
-
- int chosenBin = (Int_t) TMath::BinarySearch(50, m_wiggleLayer1, searchRand);
- double x_wigg = ((-0.98)+(chosenBin+1)*0.04)/2;
- wiggle = x_wigg*cell_dphi/4;
- }
-
- if (layer == 2){
- double m_wiggleLayer2[]={0.0127507,0.0255775,0.0395137,0.0542644,0.0695555,0.0858206,0.102274,0.119653,0.137832,0.156777,0.176938,0.197727,0.217576,0.236615,0.256605,0.277766,0.2995,0.321951,0.344663,0.367903,0.392401,0.417473,0.443514,0.470867,0.498296,0.52573,0.553114,0.57921,0.604326,0.628822,0.652191,0.674853,0.697268,0.718983,0.739951,0.759866,0.778877,0.798762,0.819559,0.839789,0.858923,0.877327,0.894831,0.911693,0.92821,0.94391,0.959156,0.973593,0.986752,1};
-
- int chosenBin = (Int_t) TMath::BinarySearch(50, m_wiggleLayer2, searchRand);
- double x_wigg = ((-0.98)+(chosenBin+1)*0.04)/2;
- wiggle = x_wigg*cell_dphi;
- }
-
- if (layer == 3){
- double m_wiggleLayer3[]={0.0217932,0.0438502,0.0670992,0.091085,0.11651,0.143038,0.169524,0.196205,0.222944,0.249703,0.276629,0.303559,0.33034,0.356842,0.383579,0.410385,0.437272,0.464214,0.49118,0.518202,0.545454,0.572667,0.600037,0.627544,0.655072,0.6826,0.709824,0.733071,0.754764,0.775672,0.793834,0.810904,0.828219,0.844119,0.858339,0.871248,0.882485,0.894889,0.907955,0.920289,0.931136,0.941039,0.949844,0.957641,0.965787,0.97392,0.981706,0.988892,0.994527,1};
-
- int chosenBin = (Int_t) TMath::BinarySearch(50, m_wiggleLayer3, searchRand);
- double x_wigg = ((-0.98)+(chosenBin+1)*0.04)/2;
- wiggle = x_wigg*cell_dphi;
+ double x=eta_particle;
+ double delta,sech;
+ double c,d;
+ const double sign= forMaximum ? 1. : -1.;
+
+ double y=0,max=0;
+
+ const double epsilon=sign*0.2;
+ //std::cout << "Serching for starting point" << std::endl;
+ //std::cout << "Epsilon" << epsilon << std::endl;
+ int count=0;
+ for(;count<1000;count++)
+ {
+ x+=epsilon;
+ delta=x-eta_particle;
+ c=exp(x);
+ d=1./c;
+ sech=2./(c+d);
+ y=transformationEtaToMM(R,delta,sech);
+ //std::cout << "y " << y << " max " << max << " sign*(y - max) " << sign*(y - max) << std::endl;
+ if(sign*(y - max)>0)max=y;
+ else break;
}
+ x-=epsilon;
+
+ return x;
+}
- return wiggle;
+void ISF::FastCaloSimSvcV2::findEtaMMRange(const double R,const double eta_particle,double & etaMin,double & etaMax)
+{
+
+ // const double epsilon1=0.1;
+ double delta;
+ double x=eta_particle;
+ // double secondDerivative(0);
+ double c,d,sech,tanh;
+
+ double startingPointForMaximum=findStartingPoint(R,eta_particle,true);
+ double startingPointForMinimum=findStartingPoint(R,eta_particle,false);
+
+
+ c=exp(startingPointForMaximum);
+ d=1./c;
+ sech=2./(c+d);
+ tanh=(c-d)/2.*sech;
+ delta=startingPointForMaximum-eta_particle;
+ //std::cout<< "Starting point for maximum: " << startingPointForMaximum << std::endl;
+ //std::cout<< "Function evaluated in starting point for maximum: " << transformationEtaToMM(R,delta,sech) << std::endl;
+ //std::cout<< "First derivative in starting point for maximum: " << transformationEtaToMMFirstDerivative(R,delta,sech,tanh) << std::endl;
+ //std::cout<< "Second derivative in starting point for maximum: " << transformationEtaToMMSecondDerivative(R,delta,sech,tanh) << std::endl;
+
+
+ c=exp(startingPointForMinimum);
+ d=1./c;
+ sech=2./(c+d);
+ tanh=(c-d)/2.*sech;
+ delta=startingPointForMinimum-eta_particle;
+ //std::cout<< "Starting point for minimum: " << startingPointForMinimum << std::endl;
+ //std::cout<< "Function evaluated in starting point for minimum: " << transformationEtaToMM(R,delta,sech) << std::endl;
+ //std::cout<< "First derivative in starting point for minimum: " << transformationEtaToMMFirstDerivative(R,delta,sech,tanh) << std::endl;
+ //std::cout<< "Second derivative in starting point for minimum: " << transformationEtaToMMSecondDerivative(R,delta,sech,tanh) << std::endl;
+
+ x=startingPointForMaximum;
+ const double epsilon=0.001;
+ int count =0;
+ double a,b;
+ //double e;
+ double gamma=1;
+
+ //std::cout << "Searching for maximum" << std::endl;
+ do
+ {
+ c=exp(x);
+ d=1./c;
+ sech=2./(c+d);
+ tanh=(c-d)/2.*sech;
+ delta=x-eta_particle;
+ a = transformationEtaToMMFirstDerivative(R,delta,sech,tanh);// transformation
+ b = transformationEtaToMMSecondDerivative(R,delta,sech,tanh);// first derivative of the transformation
+ //e = transformationEtaToMM(R,delta,sech);
+ x = x - gamma*a/b;
+ count++;
+
+ //std::cout<<"count "<<count<<" a "<<a<<" b "<<b<< " value " << e <<" c "<<c<<" d "<<d<<" delta "<<delta<<" x "<<x<<std::endl;
+ } while((abs(a) > epsilon) && count<1000);
+ if(x!=x || count>=1000) ATH_MSG_ERROR("Error: Maximum for delta eta[mm] range not found.");
+
+ delta=x-eta_particle;
+ c=exp(x);
+ d=1./c;
+ sech=2./(c+d);
+ etaMax=transformationEtaToMM(R,delta,sech);
+ //std::cout << " Position of maximum " << x << " Eta particle " << eta_particle << " Delta Eta max[mm] " << etaMax << std::endl;
+
+
+ x=startingPointForMinimum;
+ count =0;
+ //std::cout << "Searching for minimum" << std::endl;
+ do
+ {
+ c=exp(x);
+ d=1./c;
+ sech=2./(c+d);
+ tanh=(c-d)/2.*sech;
+ delta=x-eta_particle;
+ a = transformationEtaToMMFirstDerivative(R,delta,sech,tanh);// transformation
+ b = transformationEtaToMMSecondDerivative(R,delta,sech,tanh);// first derivative of the transformation
+ //e = transformationEtaToMM(R,delta,sech);
+ x = x - gamma*a/b;
+ count++;
+
+ //std::cout<<"count "<<count<<" a "<<a<<" b "<<b<< " value " << e <<" c "<<c<<" d "<<d<<" delta "<<delta<<" x "<<x<<std::endl;
+ } while((abs(a) > epsilon) && count<1000);
+ if(x!=x || count>=1000) ATH_MSG_ERROR("Error: Minimum for delta eta[mm] range not found.");
+ delta=x-eta_particle;
+ c=exp(x);
+ d=1./c;
+ sech=2./(c+d);
+ etaMin=transformationEtaToMM(R,delta,sech);
+ //std::cout << " Position of minimum " << x << " Eta particle " << eta_particle << " Delta Eta min[mm] " << etaMin << std::endl;
+
+}
+double ISF::FastCaloSimSvcV2::doWiggle(int layer)
+{
+
+ double wiggle = 0.0;
+ double cell_dphi = 0.0;
+
+ //Define cell size in phi
+ if(layer == 1)
+ cell_dphi = 0.0981748;
+ if(layer == 2 || layer == 3)
+ cell_dphi = 0.0245437;
+ if(layer==0 || layer>3)
+ {
+ ATH_MSG_ERROR("I am inside the wiggle calculation, but the layer is " << layer << "!");
+ return 0.0;
+ }
+
+ //Set random numbers
+ double searchRand = CLHEP::RandFlat::shoot(m_randomEngine);
+
+ //Now for layer dependant approach
+ if(layer == 1)
+ {
+ int chosenBin = (Int_t) TMath::BinarySearch(50, m_wiggleLayer1, searchRand);
+ double x_wigg = ((-0.98)+(chosenBin+1)*0.04)/2;
+ wiggle = x_wigg*cell_dphi/4;
+ }
+
+ if(layer == 2)
+ {
+ int chosenBin = (Int_t) TMath::BinarySearch(50, m_wiggleLayer2, searchRand);
+ double x_wigg = ((-0.98)+(chosenBin+1)*0.04)/2;
+ wiggle = x_wigg*cell_dphi;
+ }
+
+ if(layer == 3)
+ {
+ int chosenBin = (Int_t) TMath::BinarySearch(50, m_wiggleLayer3, searchRand);
+ double x_wigg = ((-0.98)+(chosenBin+1)*0.04)/2;
+ wiggle = x_wigg*cell_dphi;
+ }
+
+ return wiggle;
+
}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.h
index 5c993ea97417..52c7d31f99ee 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.h
@@ -14,17 +14,27 @@
// FastCaloSim includes
#include "ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h"
+#include "ISF_FastCaloSimParametrization/CaloGeometryFromCaloDDM.h"
+#include "ISF_FastCaloSimParametrization/IFastCaloSimCaloExtrapolation.h"
+#include "TrkExInterfaces/ITimedExtrapolator.h"
#include "CaloInterface/ICaloCellMakerTool.h"
#include "ISF_FastCaloSimEvent/IntArray.h"
#include "TH2D.h"
+#include "TRandom3.h"
#include "AthenaKernel/IAtRndmGenSvc.h"
+#include "AtlasDetDescr/AtlasDetectorID.h"
+#include "CaloIdentifier/LArEM_ID.h"
+#include "CaloIdentifier/LArHEC_ID.h"
+#include "CaloIdentifier/LArFCAL_ID.h"
+#include "CaloIdentifier/TileID.h"
+
namespace CLHEP
{
- class HepRandomEngine;
+ class HepRandomEngine;
}
//forward declarations
@@ -35,72 +45,134 @@ class CaloGeometryFromFile;
namespace ISF {
/** @class FastCaloSimSvcV2
- @author Elmar.Ritsch -at- cern.ch, Geraldine.Conti -at- cern.ch, Flavia.Dias -at- cern.ch
- */
- class FastCaloSimSvcV2 : public BaseSimulationSvc {
- public:
- /** Constructor with parameters */
- FastCaloSimSvcV2(const std::string& name, ISvcLocator* pSvcLocator);
-
- /** Destructor */
- virtual ~FastCaloSimSvcV2();
-
- /** Athena algorithm's interface methods */
- virtual StatusCode initialize() override final;
- virtual StatusCode finalize() override final;
-
- /** Simulation Call */
- virtual StatusCode simulate(const ISFParticle& isp) override final;
-
- /** Setup Event chain - in case of a begin-of event action is needed */
- virtual StatusCode setupEvent() override final;
-
- /** Release Event chain - in case of an end-of event action is needed */
- virtual StatusCode releaseEvent() override final;
-
- private:
-
- int hitCalc(double energy, int calolayer, int pdgid);
- double findHitEta(const double alpha,const double r, const double r_layer,const double z_particle,const double eta_particle);
-
- //** Calculate the wiggle for the hit-to-cell assignment on layers with accordion shape **//
- double doWiggle(int layer);
-
- /** FCSParam file **/
- TFile* m_paramsFile;
- TFile* m_energyDensity2D;
-
- std::string m_paramsFilename;
-
- //TFCSPCAEnergyParametrization m_epara;
-
- std::vector<TFCSPCAEnergyParametrization*> m_eparas;
-
- ToolHandleArray<ICaloCellMakerTool> m_caloCellMakerToolsSetup ;
- ToolHandleArray<ICaloCellMakerTool> m_caloCellMakerToolsRelease ;
-
- CaloCellContainer * m_theContainer;
-
- TH2F* m_histEnergyDensity2D;
-
- ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
- CLHEP::HepRandomEngine* m_randomEngine;
- std::string m_randomEngineName;
-
- CaloGeometryFromFile* m_caloGeo;
-
- std::string m_caloCellsOutputName;
-
- std::vector<double> m_rlayer;
-
- /** enable/disable random fluctuations in calorimeter layers */
- bool m_doRandomFluctuations;
-
- //** Array for the hit-to-cell assignment accordion structure fix (wiggle) **//
- //** To be moved to the conditions database at some point **//
- static const double m_wiggleLayer1[50];
- static const double m_wiggleLayer2[50];
- static const double m_wiggleLayer3[50];
+ @author Elmar.Ritsch -at- cern.ch, Geraldine.Conti -at- cern.ch, Flavia.Dias -at- cern.ch
+ */
+
+ class FastCaloSimSvcV2 : public BaseSimulationSvc
+ {
+ public:
+ /** Constructor with parameters */
+ FastCaloSimSvcV2(const std::string& name, ISvcLocator* pSvcLocator);
+
+ /** Destructor */
+ virtual ~FastCaloSimSvcV2();
+
+ /** Athena algorithm's interface methods */
+ StatusCode initialize();
+ StatusCode finalize();
+
+ /** Simulation Call */
+ StatusCode simulate(const ISFParticle& isp);
+
+ /** Setup Event chain - in case of a begin-of event action is needed */
+ StatusCode setupEvent();
+
+ /** Release Event chain - in case of an end-of event action is needed */
+ StatusCode releaseEvent();
+
+ int hitCalc(double energy, int calolayer,int pcabin, int pdgid);
+ double findHitEta(const double delta_eta_mm,const double R,const double eta_particle);
+ inline double transformationEtaToMM(const double R,const double deltaEta, const double sech)
+ {
+ return R*sech*deltaEta;
+ }
+ inline double transformationEtaToMMFirstDerivative(const double R,const double deltaEta, const double sech,const double tanh)
+ {
+ return R*sech*(1.-(deltaEta)*tanh);
+ }
+ inline double transformationEtaToMMSecondDerivative(const double R,const double deltaEta, const double sech,const double tanh)
+ {
+ return R*sech*((deltaEta)*(1-2.*sech*sech) - 2.*tanh);
+ }
+ void findEtaMMRange(const double R,const double eta_particle,double & etaMin,double & etaMax);
+ double findStartingPoint(const double R,const double eta_particle,bool forMaximum);
+
+ void LoopOverHits(const double totalEnergy, double minR = 0, double maxR = 0);
+ void TestCell();
+
+ //** Calculate the wiggle for the hit-to-cell assignment on layers with accordion shape **//
+ double doWiggle(int layer);
+
+ AtlasDetectorID* detID;
+ LArEM_ID* larID;
+ TileID* tileID;
+ LArHEC_ID* hecID;
+ LArFCAL_ID* fcalID;
+
+ int cellcheck;
+
+ /** FCSParam file **/
+ TFile* m_paramsFile;
+ TFile* m_photonFile;
+ TFile* m_elFile;
+ TFile* m_pionFile;
+ TFile* m_paramsFile_photons;
+
+ std::string m_paramsFilename;
+
+ //TFCSPCAEnergyParametrization m_epara;
+
+ std::vector<TFCSPCAEnergyParametrization*> m_eparas;
+
+ ToolHandleArray<ICaloCellMakerTool> m_caloCellMakerToolsSetup ;
+ ToolHandleArray<ICaloCellMakerTool> m_caloCellMakerToolsRelease ;
+
+ ToolHandle<IFastCaloSimCaloExtrapolation> m_FastCaloSimCaloExtrapolation;
+ ToolHandle<Trk::ITimedExtrapolator> m_extrapolator;
+
+ CaloCellContainer * m_theContainer;
+ Amg::Vector3D particle_position;
+ TH2F* m_histEnergyDensity2D;
+
+ ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
+ CLHEP::HepRandomEngine* m_randomEngine;
+ std::string m_randomEngineName;
+
+ CaloGeometryFromCaloDDM* m_caloGeo;
+
+ std::string m_caloCellsOutputName;
+
+ std::vector<double> m_rlayer;
+
+ TRandom3* m_random;
+
+ /** enable/disable random fluctuations in calorimeter layers */
+ bool m_doRandomFluctuations;
+ bool m_useFastCaloSimCaloExtrapolation;
+
+ // Settings for new shape parametrisation
+ bool m_useOneDShapeParametrisation;
+ int m_nHits;
+ TRandom *rnd;
+
+ double eta;
+ double phi;
+ double eta_isfp;
+ double phi_isfp;
+ double R;
+
+ IntArray *m_layers;
+ double r_layer;
+ double z_particle;
+ int layer;
+ double layerE;
+
+ double deltaEtaMMmin;
+ double deltaEtaMMmax;
+
+ double esum;
+ int nCells;
+ int nCells_Tile;
+
+ TFCSSimulationState simulstate;
+ // iterator over layers
+ int ilayer;
+
+ //** Array for the hit-to-cell assignment accordion structure fix (wiggle) **//
+ //** To be moved to the conditions database at some point **//
+ double m_wiggleLayer1[50];
+ double m_wiggleLayer2[50];
+ double m_wiggleLayer3[50];
};
--
GitLab