From 4ffe7b9ab71be93c096ae7c6508436253736fd6b Mon Sep 17 00:00:00 2001 From: Denis Oliveira Damazio <Denis.Oliveira.Damazio@cern.ch> Date: Tue, 13 May 2014 18:26:12 +0200 Subject: [PATCH] Fix for the python (LArL1Sim-00-10-33) --- .../LArL1Sim/LArL1Sim/LArSCL1Maker.h | 233 ++ .../LArL1Sim/LArL1Sim/LArSCSimpleMaker.h | 70 + .../LArL1Sim/LArL1Sim/LArTTL1Calib.h | 113 + .../LArL1Sim/LArL1Sim/LArTTL1Maker.h | 265 ++ LArCalorimeter/LArL1Sim/cmt/requirements | 49 + LArCalorimeter/LArL1Sim/doc/mainpage.h | 22 + .../LArL1Sim/python/LArSCL1Getter.py | 109 + .../LArL1Sim/python/LArTTL1Getter.py | 100 + LArCalorimeter/LArL1Sim/python/__init__.py | 6 + .../LArL1Sim/share/Fcal_ptweights_table7.data | 65 + .../LArL1Sim/share/LArL1Sim_G4_jobOptions.py | 12 + .../LArL1Sim/share/LArL1Sim_MC_jobOptions.py | 45 + .../share/LArL1Sim_calib_jobOptions.py | 69 + .../share/LArL1Sim_poolTestRead_jobOptions.py | 80 + .../LArL1Sim_poolTestWrite_jobOptions.py | 38 + .../share/LArL1Sim_readRDO_jobOptions.py | 100 + .../share/LArL1Sim_testCosmic_jobOptions.py | 53 + .../share/LArL1Sim_testG4_jobOptions.py | 88 + .../share/LArSCSimpleMaker_jobOptions.py | 2 + .../LArL1Sim/share/SimpleSC_From_ESD.py | 191 ++ LArCalorimeter/LArL1Sim/src/LArSCL1Maker.cxx | 864 ++++++ .../LArL1Sim/src/LArSCSimpleMaker.cxx | 210 ++ LArCalorimeter/LArL1Sim/src/LArTTL1Calib.cxx | 452 +++ LArCalorimeter/LArL1Sim/src/LArTTL1Maker.cxx | 2649 +++++++++++++++++ .../src/components/LArL1Sim_entries.cxx | 23 + .../LArL1Sim/src/components/LArL1Sim_load.cxx | 8 + 26 files changed, 5916 insertions(+) create mode 100755 LArCalorimeter/LArL1Sim/LArL1Sim/LArSCL1Maker.h create mode 100755 LArCalorimeter/LArL1Sim/LArL1Sim/LArSCSimpleMaker.h create mode 100755 LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Calib.h create mode 100755 LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Maker.h create mode 100755 LArCalorimeter/LArL1Sim/cmt/requirements create mode 100644 LArCalorimeter/LArL1Sim/doc/mainpage.h create mode 100644 LArCalorimeter/LArL1Sim/python/LArSCL1Getter.py create mode 100644 LArCalorimeter/LArL1Sim/python/LArTTL1Getter.py create mode 100644 LArCalorimeter/LArL1Sim/python/__init__.py create mode 100755 LArCalorimeter/LArL1Sim/share/Fcal_ptweights_table7.data create mode 100755 LArCalorimeter/LArL1Sim/share/LArL1Sim_G4_jobOptions.py create mode 100755 LArCalorimeter/LArL1Sim/share/LArL1Sim_MC_jobOptions.py create mode 100755 LArCalorimeter/LArL1Sim/share/LArL1Sim_calib_jobOptions.py create mode 100755 LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestRead_jobOptions.py create mode 100755 LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestWrite_jobOptions.py create mode 100644 LArCalorimeter/LArL1Sim/share/LArL1Sim_readRDO_jobOptions.py create mode 100755 LArCalorimeter/LArL1Sim/share/LArL1Sim_testCosmic_jobOptions.py create mode 100755 LArCalorimeter/LArL1Sim/share/LArL1Sim_testG4_jobOptions.py create mode 100644 LArCalorimeter/LArL1Sim/share/LArSCSimpleMaker_jobOptions.py create mode 100755 LArCalorimeter/LArL1Sim/share/SimpleSC_From_ESD.py create mode 100755 LArCalorimeter/LArL1Sim/src/LArSCL1Maker.cxx create mode 100755 LArCalorimeter/LArL1Sim/src/LArSCSimpleMaker.cxx create mode 100755 LArCalorimeter/LArL1Sim/src/LArTTL1Calib.cxx create mode 100755 LArCalorimeter/LArL1Sim/src/LArTTL1Maker.cxx create mode 100755 LArCalorimeter/LArL1Sim/src/components/LArL1Sim_entries.cxx create mode 100755 LArCalorimeter/LArL1Sim/src/components/LArL1Sim_load.cxx diff --git a/LArCalorimeter/LArL1Sim/LArL1Sim/LArSCL1Maker.h b/LArCalorimeter/LArL1Sim/LArL1Sim/LArSCL1Maker.h new file mode 100755 index 00000000000..aaf2383ef4e --- /dev/null +++ b/LArCalorimeter/LArL1Sim/LArL1Sim/LArSCL1Maker.h @@ -0,0 +1,233 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef LARL1SIM_LARSCL1MAKER_H +#define LARL1SIM_LARSCL1MAKER_H +// +======================================================================+ +// + + +// + Author ........: Denis O. Damazio + +// + Institut ......: BNL + +// + Creation date .: 18/11/2013 + +// + + +// +======================================================================+ +// +// ....... include +// + +#include "GaudiKernel/Algorithm.h" +#include "GaudiKernel/Property.h" +#include "GaudiKernel/ServiceHandle.h" + +#include "LArDigitization/LArHitEMap.h" +#include "LArElecCalib/ILArAutoCorrNoiseTool.h" +#include "LArElecCalib/ILArADC2MeVTool.h" + + +class StoreGateSvc; +class PileUpMergeSvc; +class IAtRndmGenSvc; +class ITriggerTime; +class CaloCell_SuperCell_ID; +class ICaloSuperCellIDTool; +class LArEM_ID; +class LArHEC_ID; +class LArFCAL_ID; +class LArOnline_SuperCellID; +class LArSuperCellCablingTool; +class ILArShape; +class ILArfSampl; +class ILArPedestal; +class ILArNoise; + +namespace CLHEP +{ + class HepRandomEngine; +} + +/** + @brief The aim of this algorithm is the simulation of the LAr analogue Super-Cell sums. + + It includes correct allocation of cells to Super-Cell, pulse profile as a function of time, saturation, appropriate noise addition, pile-up. <br> + The resulting signals are an input to the Preprocessor (which in turn performs digitization, filtering, bunch crossing id., noise suppression,...) <br> + Since it needs hits, the simulation only takes "simul" datasets as input, NOT digitized datasets. + + @warning although the output is not digitized, LArSCL1Maker is part of the digitization simulation. + + @author Denis O. Damazio (BNL) + */ + +class LArSCL1Maker : public Algorithm, + public IIncidentListener +{ +// +// >>>>>>>> public methods +// + public: + /** constructor */ + LArSCL1Maker(const std::string& name, ISvcLocator* pSvcLocator); + + /** destructor */ + ~LArSCL1Maker(); +// +// ..... Gaudi algorithm hooks +// + /** Read ascii files for auxiliary data (puslse shapes, noise, etc...) */ + virtual StatusCode initialize(); + /** Create LArSCL1 object + save in TES (2 containers: 1 EM, 1 hadronic) + */ + virtual StatusCode execute(); + + virtual StatusCode finalize(); + virtual void handle(const Incident&); + + + private: + +// +// >>>>>>>> private algorithm parts +// + /** initialize hit map + */ + + std::vector<float> computeSignal(const Identifier towerId, const int Ieta, const int specialCase, + std::vector<float> visEnergy, const int refTime) const; + + std::vector<float> computeNoise(const Identifier towerId, const int Ieta, + std::vector<float>& inputV) ; + + /** method called at the begining of execute() to fill the hit map */ + StatusCode fillEMap(int& totHit) ; + + /** method called at initialization to read auxiliary data from ascii files */ + StatusCode readAuxiliary(); + + int decodeInverse(int region, int eta); + + /** Method to update all conditions */ + StatusCode updateConditions(); + + /** Method to print SuperCell Conditions */ + void printConditions(const HWIdentifier& hwSC); + + /** Method for converting Hits from samples (simplified version + * of the same method in LarPileUpTool) */ + void ConvertHits2Samples(const HWIdentifier & hwSC, CaloGain::CaloGain igain, + const std::vector<std::pair<float,float> > *TimeE, + std::vector<float>& samples); + +// +// >>>>>>>> private data parts +// + + StoreGateSvc* m_detectorStore; + ServiceHandle<StoreGateSvc> m_storeGateSvc; + IChronoStatSvc* m_chronSvc; + PileUpMergeSvc* m_mergeSvc; + ServiceHandle<IAtRndmGenSvc> m_atRndmGenSvc; + std::string m_rndmEngineName; + CLHEP::HepRandomEngine* m_rndmEngine; + + /** Alorithm property: use trigger time or not*/ + bool m_useTriggerTime; + /** Alorithm property: name of the TriggerTimeTool*/ + StringProperty m_triggerTimeToolName; + /** pointer to the TriggerTimeTool */ + ITriggerTime* p_triggerTimeTool; + + int m_BeginRunPriority; + + LArSuperCellCablingTool* m_cablingSCSvc; + //CaloTriggerTowerService* m_ttSvc; + ICaloSuperCellIDTool* m_scidtool; + /** pointer to the offline TT helper */ + //const CaloLVL1_ID* m_lvl1Helper; + const CaloCell_SuperCell_ID* m_scHelper; + /** pointer to the online LAr helper */ + const LArOnline_SuperCellID* m_OnlSCHelper; + + + /** number of sampling (in depth) */ + static const short s_NBDEPTHS = 4 ; + /** number of samples in TTL1s */ + static const short s_NBSAMPLES = 7 ; + /** max number of samples in pulse shape */ + static const short s_MAXSAMPLES = 24 ; + /** peak position */ + static const short s_PEAKPOS = 3 ; + /** number of eta bins */ + static const short s_NBETABINS = 15 ; + /** number of energies at which saturation is described (em) */ + static const short s_NBENERGIES = 12 ; + + /** auxiliary HEC data: auto-correlation matrix */ + std::vector< std::vector<float> > m_autoCorrHec ; + + /** hit map */ + LArHitEMap* m_hitmap; // map of hits in cell + /** list of hit containers */ + std::vector <std::string> m_HitContainer; + +/** algorithm property: sub-detectors to be simulated */ + std::string m_SubDetectors; +/** algorithm property: container name for the EM TTL1s */ + std::string m_SCL1ContainerName; + +#ifdef DONTDO +/** algorithm property: container name of the EMB hits */ + std::string m_EmBarrelHitContainerName; +/** algorithm property: container name of the EMEC hits */ + std::string m_EmEndCapHitContainerName; +/** algorithm property: container name of the HEC hits */ + std::string m_HecHitContainerName; +/** algorithm property: container name of the FCAL hits */ + std::string m_ForWardHitContainerName; +#endif + + +/** algorithm property: noise (in all sub-detectors) is on if true */ + bool m_NoiseOnOff; +/** algorithm property: pile up or not */ + bool m_PileUp; +/** algorithm property: no calibration mode for EM towers */ + bool m_noEmCalibMode; +/** algorithm property: no calibration mode for had towers */ + bool m_noHadCalibMode; +/** algorithm property: debug threshold */ + float m_debugThresh; +/** algorithm property: switch chrono on */ + bool m_chronoTest; + +/** Conditions (shape) of SuperCell */ + const ILArShape* m_shapes; +/** Conditions (fSamples) of SuperCell */ + const ILArfSampl* m_fracS; +/** Conditions (LArPedestalSC) of SuperCell */ + const ILArPedestal* m_PedestalSC; +/** Conditions (LArNoiseSC) of SuperCell */ + const ILArNoise* m_NoiseSC; +/** Special Tool for AutoCorrelation sqrt */ + ToolHandle<ILArAutoCorrNoiseTool> m_autoCorrNoiseTool; +/** Special Tool for AutoCorrelation sqrt */ + ToolHandle<ILArADC2MeVTool> m_adc2mevTool; +/** key for fSampl conditions */ + std::string m_fSamplKey; +/** key for Shape conditions */ + std::string m_shapesKey; +/** key for Noise conditions */ + std::string m_noiseKey; +/** key for Pedestal conditions */ + std::string m_pedestalKey; + + /** output number of samples */ + unsigned int m_nSamples; + /** output first samples */ + unsigned int m_firstSample; + /** m_first */ + bool m_first; + +}; + +#endif + diff --git a/LArCalorimeter/LArL1Sim/LArL1Sim/LArSCSimpleMaker.h b/LArCalorimeter/LArL1Sim/LArL1Sim/LArSCSimpleMaker.h new file mode 100755 index 00000000000..b2c2349b623 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/LArL1Sim/LArSCSimpleMaker.h @@ -0,0 +1,70 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef LARL1SIM_LARSCSIMPLEMAKER_H +#define LARL1SIM_LARSCSIMPLEMAKER_H +// +======================================================================+ +// + + +// + Author ........: Denis Oliveira Damazio + +// + Institute ......: BNL + +// + Creation date .: 17/06/2012 + +// + + +// +======================================================================+ +// +// ....... include +// + +#include "AthenaBaseComps/AthAlgorithm.h" +#include "GaudiKernel/ToolHandle.h" +#include "DataModel/DataPool.h" +#include "AthenaKernel/IOVSvcDefs.h" + + +class CaloSuperCellDetDescrManager; +class ICaloSuperCellIDTool; +class CaloDetDescrElement; +class CaloCell; +class CaloIdManager; + +/** + * @brief Make a dummy supercell @c LArRawChannelContainer from a + * @c CaloCallContainer. + */ +class LArSCSimpleMaker : public AthAlgorithm +{ +public: + /// Standard Gaudi algorithm constructor. + LArSCSimpleMaker(const std::string& name, ISvcLocator* pSvcLocator); + + + /// Standard Gaudi initialize method. + virtual StatusCode initialize(); + + /// Algorithm execute method. + virtual StatusCode execute(); + + +private: + /// Property: Offline / supercell mapping tool. + ToolHandle<ICaloSuperCellIDTool> m_scidtool; + + /// Property: SG key for the input calorimeter cell container. + std::string m_cellContainer; + + /// Property: SG key for the output supercell LAr channel container. + std::string m_sCellContainer; + + /// Geometry manager. + const CaloSuperCellDetDescrManager* m_sem_mgr; + + /// Entry point for calorimeter ID helpers. + const CaloIdManager* m_calo_id_manager; + + /// Data Pool + DataPool<CaloCell> m_dataPool; + +}; + + +#endif // not LARL1SIM_LARSCSIMPLEMAKER_H diff --git a/LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Calib.h b/LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Calib.h new file mode 100755 index 00000000000..df450d848b6 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Calib.h @@ -0,0 +1,113 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef LARL1SIM_LARTTL1CALIB_H +#define LARL1SIM_LARTTL1CALIB_H +// +======================================================================+ +// + + +// + Author ........: F. Ledroit + +// + Institut ......: ISN Grenoble + +// + Creation date .: 20/04/2004 + +// + + +// +======================================================================+ +// +// ....... include +// + +#include "GaudiKernel/Algorithm.h" +#include "GaudiKernel/NTuple.h" + + +class StoreGateSvc; +class CaloLVL1_ID; + +/** + @brief This algorithm is meant to be run in standalone only for TTL1 calibration purposes + */ +class LArTTL1Calib : public Algorithm +{ +// +// >>>>>>>> public methods +// + public: + /** constructor */ + LArTTL1Calib(const std::string& name, ISvcLocator* pSvcLocator); + +/** destructor */ + ~LArTTL1Calib(); +// +// ..... Gaudi algorithm hooks +// + virtual StatusCode initialize(); + virtual StatusCode execute(); + virtual StatusCode finalize(); + + private: + +// +// >>>>>>>> private algorithm parts +// + + void decodeInverseTTChannel(int region, int layer, int eta, int & bec, int & emhf, int & Ieta) const; + +// +// >>>>>>>> private data parts +// + + float m_thresholdGeVTTs ; + int m_maxNtt ; + +// pointer to ntuple + NTuple::Tuple* m_pnt; + +// variables to be in ntuple. + + NTuple::Item<long> m_ntotem; + NTuple::Item<float> m_etotem; + NTuple::Item<long> m_n3x3em; + NTuple::Item<float> m_e3x3em; + NTuple::Item<float> m_emostem; + NTuple::Item<long> m_bec_mostEem; + NTuple::Item<long> m_emfcal_mostEem; + NTuple::Item<long> m_Ieta_mostEem; + NTuple::Item<long> m_Lphi_mostEem; + NTuple::Item<long> m_nttem; + NTuple::Array<float> m_ettem; + NTuple::Array<long> m_becttem; + NTuple::Array<long> m_emfcalttem; + NTuple::Array<long> m_Ietattem; + NTuple::Array<long> m_Lphittem; + + NTuple::Item<long> m_ntothad; + NTuple::Item<float> m_etothad; + NTuple::Item<long> m_n3x3had; + NTuple::Item<float> m_e3x3had; + NTuple::Item<float> m_emosthad; + NTuple::Item<long> m_hecfcal_mostEhad; + NTuple::Item<long> m_Ieta_mostEhad; + NTuple::Item<long> m_Lphi_mostEhad; + NTuple::Item<long> m_ntthad; + NTuple::Array<float> m_etthad; + NTuple::Array<long> m_hecfcaltthad; + NTuple::Array<long> m_Ietatthad; + NTuple::Array<long> m_Lphitthad; + + StoreGateSvc* m_storeGateSvc; + const CaloLVL1_ID* m_lvl1Helper; + +}; + +#endif + + + + + + + + + + + + diff --git a/LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Maker.h b/LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Maker.h new file mode 100755 index 00000000000..31f85b481a8 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/LArL1Sim/LArTTL1Maker.h @@ -0,0 +1,265 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef LARL1SIM_LARTTL1MAKER_H +#define LARL1SIM_LARTTL1MAKER_H +// +======================================================================+ +// + + +// + Author ........: F. Ledroit + +// + Institut ......: ISN Grenoble + +// + Creation date .: 09/01/2003 + +// + + +// +======================================================================+ +// +// ....... include +// + +#include "GaudiKernel/Algorithm.h" +#include "GaudiKernel/Property.h" +#include "GaudiKernel/ServiceHandle.h" + +#include "LArDigitization/LArHitEMap.h" +#include "LArElecCalib/ILArfSampl.h" + + +class StoreGateSvc; +class PileUpMergeSvc; +class IAtRndmGenSvc; +class ITriggerTime; +class CaloTriggerTowerService; +class CaloLVL1_ID; +class LArEM_ID; +class LArHEC_ID; +class LArFCAL_ID; + +namespace CLHEP +{ + class HepRandomEngine; +} + +/** + @brief The aim of this algorithm is the simulation of the LAr analogue trigger tower sums. + + It includes correct allocation of cells to towers, pulse profile as a function of time, saturation, appropriate noise addition, pile-up. <br> + The resulting signals are an input to the Preprocessor (which in turn performs digitization, filtering, bunch crossing id., noise suppression,...) <br> + Since it needs hits, the simulation only takes "simul" datasets as input, NOT digitized datasets. + + @warning although the output is not digitized, LArTTL1Maker is part of the digitization simulation. + + @author F. Ledroit (LPSC-Grenoble) + */ + +class LArTTL1Maker : public Algorithm, + public IIncidentListener +{ +// +// >>>>>>>> public methods +// + public: + /** constructor */ + LArTTL1Maker(const std::string& name, ISvcLocator* pSvcLocator); + + /** destructor */ + ~LArTTL1Maker(); +// +// ..... Gaudi algorithm hooks +// + /** Read ascii files for auxiliary data (puslse shapes, noise, etc...) */ + virtual StatusCode initialize(); + /** Create LArTTL1 object + save in TES (2 containers: 1 EM, 1 hadronic) + */ + virtual StatusCode execute(); + + virtual StatusCode finalize(); + virtual void handle(const Incident&); + + + private: + +// +// >>>>>>>> private algorithm parts +// + StatusCode retrieveDatabase(); + + /** initialize hit map + */ + + virtual StatusCode initHitMap(); + + std::vector<float> computeSignal(const Identifier towerId, const int Ieta, const int specialCase, + std::vector<float> visEnergy, const int refTime) const; + + std::vector<float> computeNoise(const Identifier towerId, const int Ieta, + std::vector<float>& inputV) ; + + /** method called at the begining of execute() to fill the hit map */ + StatusCode fillEMap(int& totHit) ; + + /** method called at initialization to read auxiliary data from ascii files */ + StatusCode readAuxiliary(); + + int decodeInverse(int region, int eta); + +// +// >>>>>>>> private data parts +// + + StoreGateSvc* m_detectorStore; + ServiceHandle<StoreGateSvc> m_storeGateSvc; + IChronoStatSvc* m_chronSvc; + PileUpMergeSvc* m_mergeSvc; + ServiceHandle<IAtRndmGenSvc> m_atRndmGenSvc; + std::string m_rndmEngineName; + CLHEP::HepRandomEngine* m_rndmEngine; + + /** Alorithm property: use trigger time or not*/ + bool m_useTriggerTime; + /** Alorithm property: name of the TriggerTimeTool*/ + StringProperty m_triggerTimeToolName; + /** pointer to the TriggerTimeTool */ + ITriggerTime* p_triggerTimeTool; + /** use HitEmap from detector store or no */ + bool m_useMapfromStore; + + int m_BeginRunPriority; + + LArCablingService* m_cablingSvc; + CaloTriggerTowerService* m_ttSvc; + /** pointer to the offline TT helper */ + const CaloLVL1_ID* m_lvl1Helper; + /** pointer to the offline EM helper */ + const LArEM_ID* m_emHelper; + /** pointer to the offline HEC helper */ + const LArHEC_ID* m_hecHelper; + /** pointer to the offline FCAL helper */ + const LArFCAL_ID* m_fcalHelper; + /** Sampling fractions retrieved from DB */ + const DataHandle<ILArfSampl> m_dd_fSampl; + + /** number of sampling (in depth) */ + static const short s_NBDEPTHS = 4 ; + /** number of samples in TTL1s */ + static const short s_NBSAMPLES = 7 ; + /** max number of samples in pulse shape */ + static const short s_MAXSAMPLES = 24 ; + /** peak position */ + static const short s_PEAKPOS = 3 ; + /** number of eta bins */ + static const short s_NBETABINS = 15 ; + /** number of energies at which saturation is described (em) */ + static const short s_NBENERGIES = 12 ; + + /** auxiliary EM data: reference energies for saturation simulation */ + std::vector<float> m_refEnergyEm ; + /** auxiliary EM data: pulse shapes */ + std::vector<std::vector<float> > m_pulseShapeEm ; + /** auxiliary EM data: pulse shape derivative */ + std::vector<std::vector<float> > m_pulseShapeDerEm ; + /** auxiliary EM data: auto-correlation matrix */ + std::vector<float> m_autoCorrEm ; + + /** auxiliary EMBarrel data: sin(theta) */ + std::vector<float> m_sinThetaEmb ; + /** auxiliary EMBarrel data: calibration coefficient */ + std::vector<float> m_calibCoeffEmb ; + /** auxiliary EMBarrel data: noise rms */ + std::vector<float> m_noiseRmsEmb ; + // if later we want to disentangle between pre-sum and summing electronic noise... + // std::vector<std::vector<float> > m_noiseRmsEmb ; + + /** auxiliary EMEC data: sin(theta) */ + std::vector<float> m_sinThetaEmec ; + /** auxiliary EMEC data: calibration coeeficient */ + std::vector<float> m_calibCoeffEmec ; + /** auxiliary EMEC data: noise rms */ + std::vector<float> m_noiseRmsEmec ; + // if later we want to disentangle between pre-sum and summing electronic noise... + // std::vector<std::vector<float> > m_noiseRmsEmec ; + + /** auxiliary HEC data: pulse shape */ + std::vector<float> m_pulseShapeHec ; + /** auxiliary HEC data: pulse shape derivative */ + std::vector<float> m_pulseShapeDerHec ; + /** auxiliary HEC data: calibration coefficients */ + std::vector<float> m_calibCoeffHec ; + /** auxiliary HEC data: sin(theta) */ + std::vector<float> m_sinThetaHec ; + /** auxiliary HEC data: saturation energy */ + std::vector<float> m_satEnergyHec ; + /** auxiliary HEC data: noise rms */ + std::vector<float> m_noiseRmsHec ; + /** auxiliary HEC data: auto-correlation matrix */ + std::vector< std::vector<float> > m_autoCorrHec ; + + /** auxiliary FCAL data: relative gains */ + std::vector<float> m_cellRelGainFcal ; + /** auxiliary FCAL data: pulse shapes */ + std::vector< std::vector<float> > m_pulseShapeFcal ; + /** auxiliary FCAL data: pulse shape derivatives */ + std::vector< std::vector<float> > m_pulseShapeDerFcal ; + /** auxiliary FCAL data: calibration coefficients */ + std::vector< std::vector<float> > m_calibCoeffFcal ; + /** auxiliary FCAL data: calibration coefficients */ + std::vector<float> m_calibCoeffFcalEm ; + /** auxiliary FCAL data: calibration coefficients */ + std::vector<float> m_calibCoeffFcalHad ; + /** auxiliary FCAL data: noise rms */ + // std::vector<float> m_noiseRmsFcal ; + std::vector< std::vector<float> > m_noiseRmsFcal ; + /** auxiliary FCAL data: auto-correlation matrix */ + std::vector<float> m_autoCorrFcal ; + + /** hit map */ + LArHitEMap* m_hitmap; // map of hits in cell + /** list of hit containers */ + std::vector <std::string> m_HitContainer; + + +/** algorithm property: sub-detectors to be simulated */ + std::string m_SubDetectors; +/** algorithm property: container name for the EM TTL1s */ + std::string m_EmTTL1ContainerName; +/** algorithm property: container name for the HAD TTL1s */ + std::string m_HadTTL1ContainerName; + +/** algorithm property: container name of the EMB hits */ + std::string m_EmBarrelHitContainerName; +/** algorithm property: container name of the EMEC hits */ + std::string m_EmEndCapHitContainerName; +/** algorithm property: container name of the HEC hits */ + std::string m_HecHitContainerName; +/** algorithm property: container name of the FCAL hits */ + std::string m_ForWardHitContainerName; + + +/** algorithm property: noise (in all sub-detectors) is on if true */ + bool m_NoiseOnOff; +/** algorithm property: pile up or not */ + bool m_PileUp; +/** algorithm property: no calibration mode for EM towers */ + bool m_noEmCalibMode; +/** algorithm property: no calibration mode for had towers */ + bool m_noHadCalibMode; +/** algorithm property: debug threshold */ + float m_debugThresh; +/** algorithm property: switch chrono on */ + bool m_chronoTest; +/** key to access of fSamplKey */ + std::string m_fSamplKey; +}; + +#endif + + + + + + + + + + + + diff --git a/LArCalorimeter/LArL1Sim/cmt/requirements b/LArCalorimeter/LArL1Sim/cmt/requirements new file mode 100755 index 00000000000..1944597dbc4 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/cmt/requirements @@ -0,0 +1,49 @@ +package LArL1Sim + +author Fabienne Ledroit <ledroit@lpsc.in2p3.fr> + +use AtlasPolicy AtlasPolicy-* +use GaudiInterface GaudiInterface-* External +use AthenaBaseComps AthenaBaseComps-* Control + + +use LArDigitization LArDigitization-* LArCalorimeter +use LArElecCalib LArElecCalib-* LArCalorimeter +use AthenaKernel AthenaKernel-* Control +use DataModel DataModel-* Control + +private +use AtlasCLHEP AtlasCLHEP-* External +use AtlasBoost AtlasBoost-* External +use AthenaKernel AthenaKernel-* Control +use StoreGate StoreGate-* Control +use PileUpTools PileUpTools-* Control +use PathResolver PathResolver-* Tools +use EventInfo EventInfo-* Event +use CaloIdentifier CaloIdentifier-* Calorimeter +use CaloTriggerTool CaloTriggerTool-* Calorimeter +use CaloDetDescr CaloDetDescr-* Calorimeter +use CaloEvent CaloEvent-* Calorimeter +use LArIdentifier LArIdentifier-* LArCalorimeter +use LArTools LArTools-* LArCalorimeter +use LArSimEvent LArSimEvent-* LArCalorimeter +use LArRawEvent LArRawEvent-* LArCalorimeter +use AtlasCLHEP_RandomGenerators AtlasCLHEP_RandomGenerators-* Simulation/Tools +use GeoModelInterfaces GeoModelInterfaces-* DetectorDescription/GeoModel + +end_private + + +library LArL1Sim "*.cxx -s=../src/components *.cxx " +apply_pattern component_library + +apply_pattern declare_joboptions files="*.py" + +apply_pattern declare_runtime files="Fcal_ptweights_table7.data" + +apply_pattern declare_python_modules files="*.py" + +private + +macro_append use_cppflags -ftemplate-depth-99 +#macro_remove cppflags -pedantic-errors diff --git a/LArCalorimeter/LArL1Sim/doc/mainpage.h b/LArCalorimeter/LArL1Sim/doc/mainpage.h new file mode 100644 index 00000000000..c9fbfec918c --- /dev/null +++ b/LArCalorimeter/LArL1Sim/doc/mainpage.h @@ -0,0 +1,22 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +/** + +@mainpage + +The aim of the main algorithm of the present package (LArTTL1Maker) is the simulation of the LAr analogue trigger tower sums. <br> + + It includes correct allocation of cells to towers, pulse profile as a function of time, saturation, appropriate noise addition, pile-up. <br> + The resulting signals are an input to the Preprocessor (which in turn performs digitization, filtering, bunch crossing id., noise suppression,...) <br> + Since it needs hits, the simulation only takes "simul" datasets as input, NOT digitized datasets. <br> + +There is a second algorithm, LArTTL1Calib, which is meant to be run in standalone in order to perform the energy calibration of the Trigger Towers. + + +@htmlinclude used_packages.html + +@include requirements + +*/ diff --git a/LArCalorimeter/LArL1Sim/python/LArSCL1Getter.py b/LArCalorimeter/LArL1Sim/python/LArSCL1Getter.py new file mode 100644 index 00000000000..985ecc747c2 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/python/LArSCL1Getter.py @@ -0,0 +1,109 @@ +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + +# LArSCL1 creation from LArHits with LArSCL1Maker algorithm + +from AthenaCommon.Logging import logging +from RecExConfig.Configured import Configured +from AthenaCommon.AlgSequence import AlgSequence +topSequence = AlgSequence() +import traceback + +class LArSCL1Getter ( Configured ) : + + + def configure(self): + mlog = logging.getLogger( 'LArSCL1Getter::configure:' ) + mlog.info ('entering') + + # get handle to upstream object + try: + from LArL1Sim.LArSCL1Getter import LArSCL1Getter + theLArSCL1Getter=LArSCL1Getter() + except: + mlog.error("could not get handle to LArSCL1Getter Quit") + print traceback.format_exc() + return False + + if not theLArSCL1Getter.usable(): + if not self.ignoreConfigError(): + mlog.error("LArSCL1Getter unusable. Quit.") + return False + else: + mlog.error("LArSCL1Getter unusable. Continue nevertheless") + + # Instantiation of the C++ algorithm + try: + from LArL1Sim.LArL1SimConf import LArSCL1Maker + except: + mlog.error("could not import LArL1Sim.LArSCL1Maker") + print traceback.format_exc() + return False + + theLArSCL1Maker=LArSCL1Maker() + from AthenaCommon.AppMgr import ToolSvc + from LArRecUtils.LArAutoCorrNoiseSCToolDefault import LArAutoCorrNoiseSCToolDefault + theLArAutoCorrNoiseSCTool = LArAutoCorrNoiseSCToolDefault() + ToolSvc+=theLArAutoCorrNoiseSCTool + theLArSCL1Maker.AutoCorrNoiseTool = theLArAutoCorrNoiseSCTool + + from LArRecUtils.LArADC2MeVSCToolDefault import LArADC2MeVSCToolDefault + theLArADC2MeVSCTool = LArADC2MeVSCToolDefault() + ToolSvc+=theLArADC2MeVSCTool + theLArSCL1Maker.ADC2MeVTool = theLArADC2MeVSCTool + + theLArSCL1Maker.SCL1ContainerName = "LArDigitSCL1" + + self._LArSCL1Maker = theLArSCL1Maker ; + + from AthenaCommon.AlgSequence import AlgSequence + topSequence = AlgSequence() + + # check if LArdigitization is run before. If yes, uses hit map from detector store produces from lardigitization + from Digitization.DigitizationFlags import jobproperties + from AthenaCommon.DetFlags import DetFlags + if DetFlags.digitize.LAr_on(): + mlog.info("Using hit map from LArDigitMaker algoritm") + else: + mlog.info("digitmaker1 not found in topSequence, using own map in LArSCL1Maker") + return False + + # now add algorithm to topSequence + # this should always come at the end + from IOVDbSvc.CondDB import conddb + if ( conddb.isMC and not conddb.folderRequested('/LAR/IdentifierOfl/OnOffIdMap_SC') ) : + conddb.addFolder("LAR_OFL","<tag>LARIdentifierOflOnOffIdMap_SC-000</tag>/LAR/IdentifierOfl/OnOffIdMap_SC") + if ( conddb.isMC and not conddb.folderRequested('/LAR/ElecCalibMCSC/fSampl') ) : + conddb.addFolder("LAR_OFL","<tag>LARElecCalibMCSCfSampl-000</tag>/LAR/ElecCalibMCSC/fSampl") + if ( conddb.isMC and not conddb.folderRequested('/LAR/ElecCalibMCSC/Pedestal') ) : + conddb.addFolder("LAR_OFL","<tag>LARElecCalibMCSCPedestal-000</tag>/LAR/ElecCalibMCSC/Pedestal") + if ( conddb.isMC and not conddb.folderRequested('/LAR/ElecCalibMCSC/Noise') ) : + conddb.addFolder("LAR_OFL","<tag>LARElecCalibMCSCNoise-000</tag>/LAR/ElecCalibMCSC/Noise") + if ( conddb.isMC and not conddb.folderRequested('/LAR/ElecCalibMCSC/Shape') ) : + conddb.addFolder("LAR_OFL","<tag>LARElecCalibMCSCShape-000</tag>/LAR/ElecCalibMCSC/Shape") + + from LArRecUtils.LArRecUtilsConf import LArFlatConditionSvc + from AthenaCommon.AppMgr import ServiceMgr as svcMgr + if not hasattr(svcMgr,"LArFlatConditionSvc"): + svcMgr+=LArFlatConditionSvc() + svcMgr.ProxyProviderSvc.ProviderNames += [ "LArFlatConditionSvc" ] + svcMgr.LArFlatConditionSvc.DoSuperCells=True + + mlog.info(" now adding to topSequence") + topSequence += theLArSCL1Maker + + + return True + + def LArSCL1Maker(self): + return self._LArSCL1Maker + + def outputKey(cls): + return cls._output[cls._outputType] + + def outputType(cls): + return cls._outputType + + def outputTypeKey(self): + return str(self.outputType()+"#"+self.outputKey()) + + diff --git a/LArCalorimeter/LArL1Sim/python/LArTTL1Getter.py b/LArCalorimeter/LArL1Sim/python/LArTTL1Getter.py new file mode 100644 index 00000000000..d2488d39c36 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/python/LArTTL1Getter.py @@ -0,0 +1,100 @@ +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + +# Author: F. Ledroit (ledroit@lpsc.in2p3.fr) +# LArTTL1 creation from LArHits with LArTTL1Maker algorithm + +from AthenaCommon.Logging import logging +from RecExConfig.Configured import Configured +from AthenaCommon.AlgSequence import AlgSequence +topSequence = AlgSequence() +import traceback + +class LArTTL1Getter ( Configured ) : + + + def configure(self): + mlog = logging.getLogger( 'LArTTL1Getter::configure:' ) + mlog.info ('entering') + + # get handle to upstream object + try: + from LArL1Sim.LArTTL1Getter import LArTTL1Getter + theLArTTL1Getter=LArTTL1Getter() + except: + mlog.error("could not get handle to LArTTL1Getter Quit") + print traceback.format_exc() + return False + + if not theLArTTL1Getter.usable(): + if not self.ignoreConfigError(): + mlog.error("LArTTL1Getter unusable. Quit.") + return False + else: + mlog.error("LArTTL1Getter unusable. Continue nevertheless") + + # Instantiation of the C++ algorithm + try: + from LArL1Sim.LArL1SimConf import LArTTL1Maker + except: + mlog.error("could not import LArL1Sim.LArTTL1Maker") + print traceback.format_exc() + return False + + theLArTTL1Maker=LArTTL1Maker() + self._LArTTL1Maker = theLArTTL1Maker ; + + # Configure LArTTL1Maker here + #theLArTTL1Maker.SubDetectors="LAr_All" + #theLArTTL1Maker.EmBarrelHitContainerName="LArHitEMB" + #theLArTTL1Maker.EmEndCapHitContainerName="LArHitEMEC" + #theLArTTL1Maker.HecHitContainerName="LArHitHEC" + #theLArTTL1Maker.ForWardHitContainerName="LArHitFCAL" + + #theLArTTL1Maker.EmTTL1ContainerName="LArTTL1EM" + #theLArTTL1Maker.HadTTL1ContainerName="LArTTL1HAD" + + #theLArTTL1Maker.NoiseOnOff=true + + #theLArTTL1Maker.PileUp=false + #theLArTTL1Maker.UseTriggerTime=false + #theLArTTL1Maker.TriggerTimeToolName="CosmicTriggerTimeTool" + + #theLArTTL1Maker.EmBarrelCalibrationCoeffs=[1.03, 1.024, 1.019, 1.02, 1.02, 1.024, 1.03, 1.046, 1.06, 1.053, 1.057, 1.062, 1.063, 1.076, 1.176] + #theLArTTL1Maker.EmEndCapCalibrationCoeffs=[1.176, 1.061, 1.087, 1.015, 1.019, 1.014, 1.014, 1.009, 1.01, 1.003, 1.016, 1.003, 0.993, 1.005, 0.963] + #theLArTTL1Maker.HECCalibrationCoeffs=[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.] + #theLArTTL1Maker.EmFcalCalibrationCoeffs=[1., 1., 1., 1.] + #theLArTTL1Maker.HadFcalCalibrationCoeffs=[1., 1., 1., 1.] + + #theLArTTL1Maker.NoEmCalibrationMode=false + #theLArTTL1Maker.NoHadCalibrationMode=false + #theLArTTL1Maker.ChronoTest=false + #theLArTTL1Maker.DebugThreshold=5000. + + from AthenaCommon.AlgSequence import AlgSequence + topSequence = AlgSequence() + + # check if LArdigitization is run before. If yes, uses hit map from detector store produces from lardigitization + from Digitization.DigitizationFlags import jobproperties + from AthenaCommon.DetFlags import DetFlags + if DetFlags.digitize.LAr_on(): + mlog.info("Using hit map from LArDigitMaker algoritm") + else: + mlog.info("digitmaker1 not found in topSequence, using own map in LArTTL1Maker") + theLArTTL1Maker.useMapFromStore = False + + + + # now add algorithm to topSequence + # this should always come at the end + + mlog.info(" now adding to topSequence") + topSequence += theLArTTL1Maker + + + return True + + def LArTTL1Maker(self): + return self._LArTTL1Maker + + + diff --git a/LArCalorimeter/LArL1Sim/python/__init__.py b/LArCalorimeter/LArL1Sim/python/__init__.py new file mode 100644 index 00000000000..2bdd27b8898 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/python/__init__.py @@ -0,0 +1,6 @@ +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + +# File: LArL1Sim/python/__init__.py +__version__ = '1.0.0' +__author__ = 'Fabienne (ledroit@lpsc.in2p3.fr) ' + diff --git a/LArCalorimeter/LArL1Sim/share/Fcal_ptweights_table7.data b/LArCalorimeter/LArL1Sim/share/Fcal_ptweights_table7.data new file mode 100755 index 00000000000..b71900982b8 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/Fcal_ptweights_table7.data @@ -0,0 +1,65 @@ + 1A_4 3.263 1C_4 3.208 1E_4 3.195 1G_4 3.195 1I_4 3.263 1K_4 3.208 1M_4 3.195 1O_4 3.195 2B_2 1.604 2F_2 1.604 2J_2 1.604 2N_2 1.604 4D 0.762 4L 0.762 + 1A_3 1.681 1C_3 3.333 1E_3 3.333 1G_3 3.333 1I_3 1.681 1K_3 3.333 1M_3 3.333 1O_3 3.333 2B_1 1.770 2F_1 1.770 2J_1 1.770 2N_1 1.770 3D 1.352 3L 1.352 + 1A_2 3.407 1C_2 3.407 1E_2 3.407 1G_2 3.407 1I_2 3.407 1K_2 3.407 1M_2 3.407 1O_2 3.407 1B_2 1.925 1F_2 1.925 1J_2 1.925 1N_2 1.925 2D 1.764 2L 1.764 + 1A_1 1.721 1C_1 1.721 1E_1 1.721 1G_1 1.721 1I_1 1.721 1K_1 1.721 1M_1 1.721 1O_1 1.721 1B_1 0.994 1F_1 0.994 1J_1 0.994 1N_1 0.994 1D 1.935 1L 1.935 + 2A_4 2.571 2C_4 2.525 2E_4 2.452 2G_4 2.525 2I_4 2.571 2K_4 2.525 2M_4 2.452 2O_4 2.525 4B_2 0.503 4F_2 0.503 4J_2 0.503 4N_2 0.503 4C 0.707 4K 0.707 + 2A_3 2.750 2C_3 2.688 2E_3 2.628 2G_3 2.745 2I_3 2.750 2K_3 2.688 2M_3 2.628 2O_3 2.745 4B_1 0.766 4F_1 0.766 4J_1 0.766 4N_1 0.766 3C 1.238 3K 1.238 + 2A_2 2.962 2C_2 2.868 2E_2 2.868 2G_2 2.868 2I_2 2.962 2K_2 2.868 2M_2 2.868 2O_2 2.868 3B_2 1.141 3F_2 1.141 3J_2 1.141 3N_2 1.141 2C 1.687 2K 1.687 + 2A_1 3.130 2C_1 3.062 2E_1 2.996 2G_1 3.062 2I_1 3.130 2K_1 3.062 2M_1 2.996 2O_1 3.062 3B_1 1.375 3F_1 1.375 3J_1 1.375 3N_1 1.375 1C 1.935 1K 1.935 + 3A_4 1.644 3C_4 1.644 3E_4 1.439 3G_4 1.604 3I_4 1.644 3K_4 1.644 3M_4 1.439 3O_4 1.604 2A_2 1.501 2E_2 1.501 2I_2 1.501 2M_2 1.501 4B 0.654 4J 0.654 + 3A_3 1.925 3C_3 1.925 3E_3 1.721 3G_3 1.883 3I_3 1.925 3K_3 1.925 3M_3 1.721 3O_3 1.883 2A_1 1.681 2E_1 1.681 2I_1 1.681 2M_1 1.681 3B 1.151 3J 1.151 + 3A_2 2.154 3C_2 2.110 3E_2 2.013 3G_2 2.154 3I_2 2.154 3K_2 2.110 3M_2 2.013 3O_2 2.154 1A_2 1.840 1E_2 1.840 1I_2 1.840 1M_2 1.840 2B 1.631 2J 1.631 + 3A_1 2.355 3C_1 2.302 3E_1 2.251 3G_1 2.355 3I_1 2.355 3K_1 2.302 3M_1 2.251 3O_1 2.355 1A_1 1.969 1E_1 1.969 1I_1 1.969 1M_1 1.969 1B 1.935 1J 1.935 + 4A_4 0.336 4C_4 0.345 4E_4 0.352 4G_4 0.336 4I_4 0.336 4K_4 0.345 4M_4 0.352 4O_4 0.336 4A_2 0.224 4E_2 0.224 4I_2 0.224 4M_2 0.224 4A 0.762 4I 0.762 + 4A_3 0.406 4C_3 0.406 4E_3 0.416 4G_3 0.406 4I_3 0.406 4K_3 0.406 4M_3 0.416 4O_3 0.406 4A_1 0.700 4E_1 0.700 4I_1 0.700 4M_1 0.700 3A 1.352 3I 1.352 + 4A_2 1.013 4C_2 1.013 4E_2 0.480 4G_2 0.480 4I_2 1.013 4K_2 1.013 4M_2 0.480 4O_2 0.480 3A_2 1.041 3E_2 1.041 3I_2 1.041 3M_2 1.041 2A 1.764 2I 1.764 + 4A_1 1.305 4C_1 1.339 4E_1 1.116 4G_1 1.222 4I_1 1.305 4K_1 1.339 4M_1 1.116 4O_1 1.222 3A_1 1.305 3E_1 1.305 3I_1 1.305 3M_1 1.305 1A 1.935 1I 1.935 + 1B_4 3.195 1D_4 3.195 1F_4 3.195 1H_4 3.263 1J_4 3.195 1L_4 3.195 1N_4 3.195 1P_4 3.263 2D_2 1.501 2H_2 1.501 2L_2 1.501 2P_2 1.501 4H 0.762 4P 0.762 + 1B_3 3.333 1D_3 3.333 1F_3 3.333 1H_3 1.681 1J_3 3.333 1L_3 3.333 1N_3 3.333 1P_3 1.681 2D_1 1.681 2H_1 1.681 2L_1 1.681 2P_1 1.681 3H 1.352 3P 1.352 + 1B_2 3.407 1D_2 3.407 1F_2 3.407 1H_2 3.407 1J_2 3.407 1L_2 3.407 1N_2 3.407 1P_2 3.407 1D_2 1.840 1H_2 1.840 1L_2 1.840 1P_2 1.840 2H 1.764 2P 1.764 + 1B_1 1.721 1D_1 1.721 1F_1 1.721 1H_1 1.721 1J_1 1.721 1L_1 1.721 1N_1 1.721 1P_1 1.721 1D_1 1.969 1H_1 1.969 1L_1 1.969 1P_1 1.969 1H 1.935 1P 1.935 + 2B_4 2.525 2D_4 2.452 2F_4 2.525 2H_4 2.571 2J_4 2.525 2L_4 2.452 2N_4 2.525 2P_4 2.571 4D_2 0.224 4H_2 0.224 4L_2 0.224 4P_2 0.224 4G 0.707 4O 0.707 + 2B_3 2.745 2D_3 2.628 2F_3 2.745 2H_3 2.750 2J_3 2.745 2L_3 2.628 2N_3 2.745 2P_3 2.750 4D_1 0.700 4H_1 0.700 4L_1 0.700 4P_1 0.700 3G 1.238 3O 1.238 + 2B_2 2.868 2D_2 2.868 2F_2 2.868 2H_2 2.962 2J_2 2.868 2L_2 2.868 2N_2 2.868 2P_2 2.962 3D_2 1.041 3H_2 1.041 3L_2 1.041 3P_2 1.041 2G 1.687 2O 1.687 + 2B_1 3.062 2D_1 2.996 2F_1 3.062 2H_1 3.130 2J_1 3.062 2L_1 2.996 2N_1 3.062 2P_1 3.130 3D_1 1.305 3H_1 1.305 3L_1 1.305 3P_1 1.305 1G 1.935 1O 1.935 + 3B_4 1.604 3D_4 1.439 3F_4 1.604 3H_4 1.644 3J_4 1.604 3L_4 1.439 3N_4 1.604 3P_4 1.644 2C_2 1.604 2G_2 1.604 2K_2 1.604 2O_2 1.604 4F 0.654 4N 0.654 + 3B_3 1.883 3D_3 1.721 3F_3 1.883 3H_3 1.925 3J_3 1.883 3L_3 1.721 3N_3 1.883 3P_3 1.925 2C_1 1.770 2G_1 1.770 2K_1 1.770 2O_1 1.770 3F 1.151 3N 1.151 + 3B_2 2.154 3D_2 2.013 3F_2 2.154 3H_2 2.154 3J_2 2.154 3L_2 2.013 3N_2 2.154 3P_2 2.154 1C_2 1.925 1G_2 1.925 1K_2 1.925 1O_2 1.925 2F 1.631 2N 1.631 + 3B_1 2.355 3D_1 2.251 3F_1 2.355 3H_1 2.355 3J_1 2.355 3L_1 2.251 3N_1 2.355 3P_1 2.355 1C_1 0.994 1G_1 0.994 1K_1 0.994 1O_1 0.994 1F 1.935 1N 1.935 + 4B_4 0.336 4D_4 0.352 4F_4 0.336 4H_4 0.336 4J_4 0.336 4L_4 0.352 4N_4 0.336 4P_4 0.336 4C_2 0.503 4G_2 0.503 4K_2 0.503 4O_2 0.503 4E 0.762 4M 0.762 + 4B_3 0.406 4D_3 0.416 4F_3 0.406 4H_3 0.406 4J_3 0.406 4L_3 0.416 4N_3 0.406 4P_3 0.406 4C_1 0.766 4G_1 0.766 4K_1 0.766 4O_1 0.766 3E 1.352 3M 1.352 + 4B_2 0.480 4D_2 0.480 4F_2 0.480 4H_2 1.013 4J_2 0.480 4L_2 0.480 4N_2 0.480 4P_2 1.013 3C_2 1.141 3G_2 1.141 3K_2 1.141 3O_2 1.141 2E 1.764 2M 1.764 + 4B_1 1.222 4D_1 1.116 4F_1 1.222 4H_1 1.305 4J_1 1.222 4L_1 1.116 4N_1 1.222 4P_1 1.305 3C_1 1.375 3G_1 1.375 3K_1 1.375 3O_1 1.375 1E 1.935 1M 1.935 + 1H_4 3.263 1F_4 3.195 1D_4 3.195 1B_4 3.195 1P_4 3.263 1N_4 3.195 1L_4 3.195 1J_4 3.195 2G_2 1.604 2C_2 1.604 2O_2 1.604 2K_2 1.604 4H 0.762 4P 0.762 + 1H_3 1.681 1F_3 3.333 1D_3 3.333 1B_3 3.333 1P_3 1.681 1N_3 3.333 1L_3 3.333 1J_3 3.333 2G_1 1.770 2C_1 1.770 2O_1 1.770 2K_1 1.770 3H 1.352 3P 1.352 + 1H_2 3.407 1F_2 3.407 1D_2 3.407 1B_2 3.407 1P_2 3.407 1N_2 3.407 1L_2 3.407 1J_2 3.407 1G_2 1.925 1C_2 1.925 1O_2 1.925 1K_2 1.925 2H 1.764 2P 1.764 + 1H_1 1.721 1F_1 1.721 1D_1 1.721 1B_1 1.721 1P_1 1.721 1N_1 1.721 1L_1 1.721 1J_1 1.721 1G_1 0.994 1C_1 0.994 1O_1 0.994 1K_1 0.994 1H 1.935 1P 1.935 + 2H_4 2.571 2F_4 2.525 2D_4 2.452 2B_4 2.525 2P_4 2.571 2N_4 2.525 2L_4 2.452 2J_4 2.525 4G_2 0.503 4C_2 0.503 4O_2 0.503 4K_2 0.503 4G 0.707 4O 0.707 + 2H_3 2.750 2F_3 2.745 2D_3 2.628 2B_3 2.745 2P_3 2.750 2N_3 2.745 2L_3 2.628 2J_3 2.745 4G_1 0.766 4C_1 0.766 4O_1 0.766 4K_1 0.766 3G 1.238 3O 1.238 + 2H_2 2.962 2F_2 2.868 2D_2 2.868 2B_2 2.868 2P_2 2.962 2N_2 2.868 2L_2 2.868 2J_2 2.868 3G_2 1.141 3C_2 1.141 3O_2 1.141 3K_2 1.141 2G 1.687 2O 1.687 + 2H_1 3.130 2F_1 3.062 2D_1 2.996 2B_1 3.062 2P_1 3.130 2N_1 3.062 2L_1 2.996 2J_1 3.062 3G_1 1.375 3C_1 1.375 3O_1 1.375 3K_1 1.375 1G 1.935 1O 1.935 + 3H_4 1.644 3F_4 1.604 3D_4 1.439 3B_4 1.604 3P_4 1.644 3N_4 1.604 3L_4 1.439 3J_4 1.604 2H_2 1.501 2D_2 1.501 2P_2 1.501 2L_2 1.501 4F 0.654 4N 0.654 + 3H_3 1.925 3F_3 1.883 3D_3 1.721 3B_3 1.883 3P_3 1.925 3N_3 1.883 3L_3 1.721 3J_3 1.883 2H_1 1.681 2D_1 1.681 2P_1 1.681 2L_1 1.681 3F 1.151 3N 1.151 + 3H_2 2.154 3F_2 2.154 3D_2 2.013 3B_2 2.154 3P_2 2.154 3N_2 2.154 3L_2 2.013 3J_2 2.154 1H_2 1.840 1D_2 1.840 1P_2 1.840 1L_2 1.840 2F 1.631 2N 1.631 + 3H_1 2.355 3F_1 2.355 3D_1 2.251 3B_1 2.355 3P_1 2.355 3N_1 2.355 3L_1 2.251 3J_1 2.355 1H_1 1.969 1D_1 1.969 1P_1 1.969 1L_1 1.969 1F 1.935 1N 1.935 + 4H_4 0.336 4F_4 0.336 4D_4 0.352 4B_4 0.336 4P_4 0.336 4N_4 0.336 4L_4 0.352 4J_4 0.336 4H_2 0.224 4D_2 0.224 4P_2 0.224 4L_2 0.224 4E 0.762 4M 0.762 + 4H_3 0.406 4F_3 0.406 4D_3 0.416 4B_3 0.406 4P_3 0.406 4N_3 0.406 4L_3 0.416 4J_3 0.406 4H_1 0.700 4D_1 0.700 4P_1 0.700 4L_1 0.700 3E 1.352 3M 1.352 + 4H_2 1.013 4F_2 0.480 4D_2 0.480 4B_2 0.480 4P_2 1.013 4N_2 0.480 4L_2 0.480 4J_2 0.480 3H_2 1.041 3D_2 1.041 3P_2 1.041 3L_2 1.041 2E 1.764 2M 1.764 + 4H_1 1.305 4F_1 1.222 4D_1 1.116 4B_1 1.222 4P_1 1.305 4N_1 1.222 4L_1 1.116 4J_1 1.222 3H_1 1.305 3D_1 1.305 3P_1 1.305 3L_1 1.305 1E 1.935 1M 1.935 + 1G_4 3.195 1E_4 3.195 1C_4 3.208 1A_4 3.263 1O_4 3.195 1M_4 3.195 1K_4 3.208 1I_4 3.263 2E_2 1.501 2A_2 1.501 2M_2 1.501 2I_2 1.501 4D 0.762 4L 0.762 + 1G_3 3.333 1E_3 3.333 1C_3 3.333 1A_3 1.681 1O_3 3.333 1M_3 3.333 1K_3 3.333 1I_3 1.681 2E_1 1.681 2A_1 1.681 2M_1 1.681 2I_1 1.681 3D 1.352 3L 1.352 + 1G_2 3.407 1E_2 3.407 1C_2 3.407 1A_2 3.407 1O_2 3.407 1M_2 3.407 1K_2 3.407 1I_2 3.407 1E_2 1.840 1A_2 1.840 1M_2 1.840 1I_2 1.840 2D 1.764 2L 1.764 + 1G_1 1.721 1E_1 1.721 1C_1 1.721 1A_1 1.721 1O_1 1.721 1M_1 1.721 1K_1 1.721 1I_1 1.721 1E_1 1.969 1A_1 1.969 1M_1 1.969 1I_1 1.969 1D 1.935 1L 1.935 + 2G_4 2.525 2E_4 2.452 2C_4 2.525 2A_4 2.571 2O_4 2.525 2M_4 2.452 2K_4 2.525 2I_4 2.571 4E_2 0.224 4A_2 0.224 4M_2 0.224 4I_2 0.224 4C 0.707 4K 0.707 + 2G_3 2.745 2E_3 2.628 2C_3 2.688 2A_3 2.750 2O_3 2.745 2M_3 2.628 2K_3 2.688 2I_3 2.750 4E_1 0.700 4A_1 0.700 4M_1 0.700 4I_1 0.700 3C 1.238 3K 1.238 + 2G_2 2.868 2E_2 2.868 2C_2 2.868 2A_2 2.962 2O_2 2.868 2M_2 2.868 2K_2 2.868 2I_2 2.962 3E_2 1.041 3A_2 1.041 3M_2 1.041 3I_2 1.041 2C 1.687 2K 1.687 + 2G_1 3.062 2E_1 2.996 2C_1 3.062 2A_1 3.130 2O_1 3.062 2M_1 2.996 2K_1 3.062 2I_1 3.130 3E_1 1.305 3A_1 1.305 3M_1 1.305 3I_1 1.305 1C 1.935 1K 1.935 + 3G_4 1.604 3E_4 1.439 3C_4 1.644 3A_4 1.644 3O_4 1.604 3M_4 1.439 3K_4 1.644 3I_4 1.644 2F_2 1.604 2B_2 1.604 2N_2 1.604 2J_2 1.604 4B 0.654 4J 0.654 + 3G_3 1.883 3E_3 1.721 3C_3 1.925 3A_3 1.925 3O_3 1.883 3M_3 1.721 3K_3 1.925 3I_3 1.925 2F_1 1.770 2B_1 1.770 2N_1 1.770 2J_1 1.770 3B 1.151 3J 1.151 + 3G_2 2.154 3E_2 2.013 3C_2 2.110 3A_2 2.154 3O_2 2.154 3M_2 2.013 3K_2 2.110 3I_2 2.154 1F_2 1.925 1B_2 1.925 1N_2 1.925 1J_2 1.925 2B 1.631 2J 1.631 + 3G_1 2.355 3E_1 2.251 3C_1 2.302 3A_1 2.355 3O_1 2.355 3M_1 2.251 3K_1 2.302 3I_1 2.355 1F_1 0.994 1B_1 0.994 1N_1 0.994 1J_1 0.994 1B 1.935 1J 1.935 + 4G_4 0.336 4E_4 0.352 4C_4 0.345 4A_4 0.336 4O_4 0.336 4M_4 0.352 4K_4 0.345 4I_4 0.336 4F_2 0.503 4B_2 0.503 4N_2 0.503 4J_2 0.503 4A 0.762 4I 0.762 + 4G_3 0.406 4E_3 0.416 4C_3 0.406 4A_3 0.406 4O_3 0.406 4M_3 0.416 4K_3 0.406 4I_3 0.406 4F_1 0.766 4B_1 0.766 4N_1 0.766 4J_1 0.766 3A 1.352 3I 1.352 + 4G_2 0.480 4E_2 0.480 4C_2 1.013 4A_2 1.013 4O_2 0.480 4M_2 0.480 4K_2 1.013 4I_2 1.013 3F_2 1.141 3B_2 1.141 3N_2 1.141 3J_2 1.141 2A 1.764 2I 1.764 + 4G_1 1.222 4E_1 1.116 4C_1 1.339 4A_1 1.305 4O_1 1.222 4M_1 1.116 4K_1 1.339 4I_1 1.305 3F_1 1.375 3B_1 1.375 3N_1 1.375 3J_1 1.375 1A 1.935 1I 1.935 + diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_G4_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_G4_jobOptions.py new file mode 100755 index 00000000000..0522217fdef --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_G4_jobOptions.py @@ -0,0 +1,12 @@ +from LArL1Sim.LArTTL1Getter import * +theLArTTL1Getter = LArTTL1Getter() + +from LArL1Sim.LArL1SimConf import * +theLArTTL1Maker=LArTTL1Maker() +from Digitization.DigitizationFlags import digitizationFlags +theLArTTL1Maker.RndmSvc = digitizationFlags.rndmSvc.get_Value() +digitizationFlags.rndmSeedList.addSeed("LArTTL1Maker", 335242, 7306589 ) +#include( "LArDetDescr/LArDetDescr_joboptions.py" ) + +#include( "CaloConditions/CaloConditions_jobOptions.py" ) + diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_MC_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_MC_jobOptions.py new file mode 100755 index 00000000000..1348871cd4b --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_MC_jobOptions.py @@ -0,0 +1,45 @@ +# +===========================================================================+ +# + + +# + Author ........: F. Ledroit + +# + Institut ......: ISN Grenoble + +# + Creation date .: 04/02/2003 + +# + Subject .......: Job Option file to be used with zebra input + +# + ------------------------------------------------------------------------- + +# This jobOptions file is for creating LArTTL1 from simulated Hits. +# Not to be used for TB. +include( "LArSim/LArSim_MC_jobOptions.py" ) + +# switch neighbour initialisation off (default = true = perform init.) +#DetDescrCnvSvc.DoInitNeighbours = FALSE + +# The following are properties for LArL1Sim: +# +# -------------------------------------------------------- +# ............ declare the used top algo. +# -------------------------------------------------------- +theApp.Dlls += ["LArL1Sim"] +theApp.TopAlg += [ "LArTTL1Maker"] +# +# -------------------------------------------------------- +# ............ Algorithms Private Options +# -------------------------------------------------------- +# +# +# .................. sub-detector hit containers +# +# +LArTTL1Maker = Algorithm( "LArTTL1Maker" ) +LArTTL1Maker.SubDetectors = "LAr_All" +# LArTTL1Maker.SubDetectors = "LAr_EmBarrel" +# LArTTL1Maker.SubDetectors = "LAr_EmEndCap" +# LArTTL1Maker.SubDetectors = "LAr_Em" +# LArTTL1Maker.SubDetectors = "LAr_HEC" +# +# .................. set the TTL1 container name +# +LArTTL1Maker.EmTTL1ContainerName = "LArTTL1EM" +LArTTL1Maker.HadTTL1ContainerName = "LArTTL1HAD" +# +# ............ if true, noise added in all subdetectors (default = true) +# +#LArTTL1Maker.NoiseOnOff = FALSE diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_calib_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_calib_jobOptions.py new file mode 100755 index 00000000000..728bb8cff87 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_calib_jobOptions.py @@ -0,0 +1,69 @@ +# +===========================================================================+ +# + + +# + Author ........: F. Ledroit + +# + Institut ......: LPSC Grenoble + +# + Creation date .: 21/04/2004 Version 0.01 + +# + Subject .......: Job Option file to make calibration of the TTL1s + +# +===========================================================================+ +include( "LArL1Sim/LArL1Sim_testG4_jobOptions.py" ) + +theApp.TopAlg += [ "LArTTL1Calib"] + +#-------------------------------------------------------------- +# Private Application Configuration options +#-------------------------------------------------------------- + +#-------------------------------------------------------------- +# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) +#-------------------------------------------------------------- +LArTTL1Calib = Algorithm( "LArTTL1Calib" ) +#LArTTL1Calib.OutputLevel = 2 +svcMgr.MessageSvc.useColors = TRUE +svcMgr.MessageSvc.OutputLevel = 3 +svcMgr.MessageSvc.debugLimit = 100000 + +theLArTTL1Maker=LArTTL1Maker() +theLArTTL1Maker.OutputLevel = 3 + +theApp.EvtMax = 100000 + +#-------------------------------------------------------------- +# Select HBOOK or ROOT persistency +#-------------------------------------------------------------- +theApp.HistogramPersistency = "HBOOK" +#-------------------------------------------------------------- +# Histogram output file +#-------------------------------------------------------------- +# Specify the appropriate output file type +HistogramPersistencySvc = Service( "HistogramPersistencySvc" ) +HistogramPersistencySvc.OutputFile = "histo.hbook" +#HistogramPersistencySvc.OutputFile = "histo.rt" +#-------------------------------------------------------------- +# Ntuples +#-------------------------------------------------------------- +NTupleSvc = Service( "NTupleSvc" ) +NTupleSvc.Output = [ "FILE1 DATAFILE='tuple1.hbook' OPT='NEW'" ] +# +# +# ............ if true, noise added in all subdetectors (default = true) +# +theLArTTL1Maker.NoiseOnOff = FALSE +# +# ....... if true, NO calibration is applied to EM towers (for calibration coeffs computation !) +#.................(default = false ) +# +#theLArTTL1Maker.NoEmCalibrationMode = TRUE +# +# ....... if true, NO calibration is applied to HEC towers (for calibration coeffs computation !) +#.................(default = false ) +# +#theLArTTL1Maker.NoHadCalibrationMode = TRUE +# +# ............ transverse energy threshold for debug printout (default = 5000 MeV) +# +#theLArTTL1Maker.DebugThreshold = 1000. +#============================================================== +# +# End of job options file +# +############################################################### diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestRead_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestRead_jobOptions.py new file mode 100755 index 00000000000..4e552ad6ed1 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestRead_jobOptions.py @@ -0,0 +1,80 @@ +#*************************************************************************************+ +# + + +# + Author ........: F. Ledroit + +# + Institut ......: LPSC Grenoble + +# + Creation date .: 17/12/2004 Version 0.01 + +# + Subject .......: Job Option file to test the reading of pool persistified LArTTL1 + +# The test runs the calib algo and produces the calib ntuple + +#=====================================================================================+ + + +#setup GeoModel +DetDescrVersion = "ATLAS-GEO-02-01-00" +#DetDescrVersion = "ATLAS-GEO-03-00-00" + +from AthenaCommon.GlobalFlags import globalflags +# --- default is atlas geometry +#globalflags.DetGeo = 'atlas' +# --- set defaults +globalflags.DataSource = 'geant4' +#globalflags.InputFormat = 'pool' +# --- set geometry version +globalflags.DetDescrVersion = DetDescrVersion +# --- printout +globalflags.print_JobProperties() + +#-------------------------------------------------------------- +# Load POOL support +#-------------------------------------------------------------- +include( "AthenaPoolCnvSvc/ReadAthenaPool_jobOptions.py" ) + +# for ddcnvsvc +#include( "LArDetDescr/LArDetDescr_joboptions.py" ) +#include( "LArAthenaPool/LArAthenaPool_joboptions.py" ) +#theApp.Dlls += [ "LArAthenaPoolPoolCnv" ] + +# Define input +#EventSelector = Service( "EventSelector" ) +#EventSelector.InputCollections = [ "TTL1PoolFile2.root" ] +#EventSelector.InputCollections = [ "/afs/cern.ch/user/d/droussea/scratch0/g50GeV_rel_0_CSC-00-01-00.RDO.pool.root" ] +#EventSelector.InputCollections = [ "atlasMis_MyOutputFile.digit.RDO.00002.root" ] +ServiceMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/RDO/valid1.007003.singlepart_e_Et25.digit.RDO.e322_s484_tid027349/RDO.027349._00001.pool.root.1"] +#ServiceMgr.EventSelector.InputCollections = [ "rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/RDO/valid1.107204.singlepart_mu4.digit.RDO.e347_s462_d145_tid029243/RDO.029243._00400.pool.root.1" ] +#ServiceMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/RDO/valid1.107204.singlepart_mu4.digit.RDO.e347_s462_d144_tid029220/RDO.029220._00001.pool.root.1"] +#ServiceMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/RDO/valid1.107204.singlepart_mu4.digit.RDO.e347_s462_d147_tid030079/RDO.030079._00001.pool.root.1"] +#ServiceMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/RDO/valid1.107204.singlepart_mu4.digit.RDO.e347_s462_d133_tid027180/RDO.027180._00001.pool.root.1"] + + + +# ............ declare the used top algo. +# -------------------------------------------------------- +#load relevant libraries +theApp.Dlls += [ "LArL1Sim" ] +theApp.TopAlg += [ "LArTTL1Calib"] +#-------------------------------------------------------------- +# Private Application Configuration options +#-------------------------------------------------------------- +# Select the appropriate shared library +theApp.Dlls += [ "HbookCnv" ] +# Select HBOOK or ROOT persistency (HBOOK is default) +theApp.HistogramPersistency = "HBOOK" + +#-------------------------------------------------------------- +# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) +#-------------------------------------------------------------- +MessageSvc = Service( "MessageSvc" ) +MessageSvc.OutputLevel = 2 +MessageSvc.useColors = TRUE + +#-------------------------------------------------------------- +# Histogram output file +#-------------------------------------------------------------- +# Specify the appropriate output file type +HistogramPersistencySvc = Service( "HistogramPersistencySvc" ) +HistogramPersistencySvc.OutputFile = "histo.hbook" +#-------------------------------------------------------------- +# Ntuples +#-------------------------------------------------------------- +NTupleSvc = Service( "NTupleSvc" ) +NTupleSvc.Output = [ "FILE1 DATAFILE='tuple1.hbook' OPT='NEW'" ] + diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestWrite_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestWrite_jobOptions.py new file mode 100755 index 00000000000..34e432f79ad --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_poolTestWrite_jobOptions.py @@ -0,0 +1,38 @@ +# +===================================================================================+ +# + + +# + Author ........: F. Ledroit + +# + Institut ......: LPSC Grenoble + +# + Creation date .: 16/12/2004 Version 0.01 + +# + Subject .......: Job Option file to test the writing of pool persistified LArTTL1 + +# +===================================================================================+ + +# to read pool Data (hits) +include( "AthenaPoolCnvSvc/ReadAthenaPool_jobOptions.py" ) + +# G4 Pool input +# make sure the file listed is in PoolFileCatalog.xml in run directory +# AND that you have a copy or soft link of this file in the run directory +# it is possible to specify a list of files to be processed consecutively +EventSelector = Service( "EventSelector" ) + +#EventSelector.InputCollection = "$INFN"; +#EventSelector.InputCollections = ["rfio:/castor/cern.ch/grid/atlas/datafiles/dc2/simul/dc2.003026.simul.A0_top/dc2.003026.simul.A0_top._00001.pool.root.1","rfio:/castor/cern.ch/grid/atlas/datafiles/dc2/simul/dc2.003026.simul.A0_top/dc2.003026.simul.A0_top._00002.pool.root.1"] +#EventSelector.InputCollections = ["rfio:/castor/cern.ch/user/c/costanzo/validation/data/dc2val.e_et25_eta25.006002.g4sim/dc2val.e_et25_eta25.006002.g4sim.0000.root","rfio:/castor/cern.ch/user/c/costanzo/validation/data/dc2val.e_et25_eta25.006002.g4sim/dc2val.e_et25_eta25.006002.g4sim.0001.root"] +EventSelector.InputCollections = ["/afs/cern.ch/user/f/fledroit/scratch0/castor/dc2val.e_et25_eta25.006002.g4sim.0000.root"] + +# Number of events to be processed (default is 10) +theApp.EvtMax = 3 + +include( "LArAthenaPool/LArAthenaPool_joboptions.py" ) + +include( "LArL1Sim/LArL1Sim_G4_jobOptions.py" ) + +# Pool Output +Stream1 = Algorithm( "Stream1" ) +theApp.OutStream +=["Stream1"] +theApp.OutStreamType ="AthenaOutputStream" +Stream1.EvtConversionSvc ="AthenaPoolCnvSvc" +Stream1.OutputFile = "TTL1PoolFile2.root" +Stream1.ItemList+=["EventInfo#*"] +Stream1.ItemList+=["LArTTL1Container#*"] + diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_readRDO_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_readRDO_jobOptions.py new file mode 100644 index 00000000000..3db53c55eab --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_readRDO_jobOptions.py @@ -0,0 +1,100 @@ +# +==========================================================================================+ +# + + +# + Author ........: F. Ledroit + +# + Institut ......: LPSC Grenoble + +# + Creation date .: 06/08/2009 Version 0.01 + +# + Subject .......: Job Option file to read TTL1s in a RDO file by means of the calib algo + +# +==========================================================================================+ +#import AthenaCommon.AtlasUnixStandardJob +#import AthenaCommon.AtlasUnixGeneratorJob +from AthenaCommon.Logging import logging +from RecExConfig.Configured import Configured +from AthenaCommon.AlgSequence import AlgSequence +topSequence = AlgSequence() +import traceback + +from AthenaCommon.GlobalFlags import GlobalFlags +GlobalFlags.DetGeo.set_atlas() +GlobalFlags.DataSource.set_geant4() + +from AthenaCommon.DetFlags import DetFlags +# - Select detectors +DetFlags.ID_setOff() +DetFlags.Calo_setOn() +DetFlags.Muon_setOff() +DetFlags.Truth_setOff() +DetFlags.LVL1_setOff() +DetFlags.digitize.all_setOff() + +# following line to be cleaned away as soon as all subsequent jO have migrated +include( "LArDetDescr/LArDetDescr_joboptions.py" ) + +#setup GeoModel +DetDescrVersion="ATLAS-GEO-06-00-00" +from AthenaCommon.JobProperties import jobproperties +jobproperties.Global.DetDescrVersion = DetDescrVersion +from AtlasGeoModel import SetGeometryVersion +from AtlasGeoModel import GeoModelInit + + +from AthenaCommon.ConfigurableDb import getConfigurable +from AthenaCommon.AppMgr import ServiceMgr as svcMgr +svcMgr += getConfigurable( "ProxyProviderSvc" )() +svcMgr.ProxyProviderSvc.ProviderNames += [ "CondProxyProvider" ] +svcMgr += getConfigurable( "CondProxyProvider" )() +# la ligne suivante est necessaire, sinon le job plante !?!?!?!!!!! +svcMgr.CondProxyProvider.InputCollections = [ "LArTTCellMap-HadFcalFix.pool.root" ] + + +# Number of events to be processed (default is 10) +theApp.EvtMax = 10 + +# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) +svcMgr.MessageSvc.OutputLevel = 3 +svcMgr.MessageSvc.useColors = TRUE +svcMgr.MessageSvc.debugLimit = 100000 + + +import AthenaPoolCnvSvc.ReadAthenaPool + +#svcMgr.IOVDbSvc.GlobalTag="OFLCOND-CSC-01-00-00" +svcMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/RDO/valid1.105144.PythiaZee.digit.RDO.e380_s523_tid046366/RDO.046366._00001.pool.root.1"] +#PoolRDOInput = ["rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/RDO/valid1.105144.PythiaZee.digit.RDO.e380_s523_tid046366/RDO.046366._00001.pool.root.1"] + +try: + from LArL1Sim.LArL1SimConf import LArTTL1Calib +except: + mlog.error("could not import LArL1Sim.LArTTL1Calib") + +theLArTTL1Calib=LArTTL1Calib() +topSequence += theLArTTL1Calib + +#theApp.TopAlg += [ "LArTTL1Calib"] +#LArTTL1Calib = Algorithm( "LArTTL1Calib" ) +#-------------------------------------------------------------- +# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) +#-------------------------------------------------------------- +#LArTTL1Calib.OutputLevel = 2 + + +#-------------------------------------------------------------- +# Select HBOOK or ROOT persistency +#-------------------------------------------------------------- +theApp.HistogramPersistency = "HBOOK" +#-------------------------------------------------------------- +# Histogram output file +#-------------------------------------------------------------- +# Specify the appropriate output file type +HistogramPersistencySvc = Service( "HistogramPersistencySvc" ) +HistogramPersistencySvc.OutputFile = "histoRDO.hbook" +#HistogramPersistencySvc.OutputFile = "histo.rt" +#-------------------------------------------------------------- +# Ntuples +#-------------------------------------------------------------- +NTupleSvc = Service( "NTupleSvc" ) +NTupleSvc.Output = [ "FILE1 DATAFILE='tuple1RDO.hbook' OPT='NEW'" ] +# +# +# End of job options file +# +############################################################### diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_testCosmic_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_testCosmic_jobOptions.py new file mode 100755 index 00000000000..458222d1a7f --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_testCosmic_jobOptions.py @@ -0,0 +1,53 @@ +# +===========================================================================+ +# + + +# + Author ........: F. Ledroit + +# + Institut ......: ISN Grenoble + +# + Creation date .: 01/12/2005 Version 0.01 + +# + Subject .......: Job Option file to test LArTTL1Maker on Cosmics + +# +===========================================================================+ +# Top Level JobOptions to test LArTTl1Maker algorithm. +# Input is from LArHits (simulation) + +# to read pool Data +include( "AthenaPoolCnvSvc/ReadAthenaPool_jobOptions.py" ) + +# G4 Pool input +# make sure the file listed is in PoolFileCatalog.xml in run directory +# AND that you have a copy or soft link of this file in the run directory +# it is possible to specify a list of files to be processed consecutively +EventSelector = Service( "EventSelector" ) + +EventSelector.InputCollections = ["rfio:/castor/cern.ch/user/h/hma/g4_cosmic/comm.004920/simul/comm.004920.simul.cosmic._00000.pool.root"] +#EventSelector.InputCollections = ["rfio:/castor/cern.ch/user/f/fledroit/e-gamma_private/rome.007777.simul.piminus_pt600_eta165._00000.pool.root"] + +# Number of events to be processed (default is 10) +theApp.EvtMax = 3 + +include( "LArAthenaPool/LArAthenaPool_joboptions.py" ) + +include( "LArL1Sim/LArL1Sim_G4_jobOptions.py" ) + +# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) +#MessageSvc.OutputLevel = 2 +#LArTTL1Maker.OutputLevel = 2 + +#MessageSvc.debugLimit = 100000 + +# +# ............ if true, noise added in all subdetectors (default = true) +# +#LArTTL1Maker.NoiseOnOff = FALSE +# +# ............ if true, chrono monitoring is enabled (default = false) +# +#LArTTL1Maker.ChronoTest = TRUE +# +# ............ transverse energy threshold for debug printout (default = 5000 MeV) +# +#LArTTL1Maker.DebugThreshold = 1000. +# +# ..................... if true, retrieve trigger time (and subtract) (default = false) +#.......................this option is inactive (because useless) when PileUp is true +# +LArTTL1Maker.UseTriggerTime = TRUE +theApp.Dlls += ["CommissionUtils"] diff --git a/LArCalorimeter/LArL1Sim/share/LArL1Sim_testG4_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArL1Sim_testG4_jobOptions.py new file mode 100755 index 00000000000..6e2db0f0418 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArL1Sim_testG4_jobOptions.py @@ -0,0 +1,88 @@ +# +===========================================================================+ +# + + +# + Author ........: F. Ledroit + +# + Institut ......: ISN Grenoble + +# + Creation date .: 15/09/2004 Version 0.01 + +# + Subject .......: Job Option file to test LArTTL1Maker on G4 + +# +===========================================================================+ +# Top Level JobOptions to test LArTTl1Maker algorithm. +# Input is from LArHits (simulation) + +DetDescrVersion = "ATLAS-GEO-03-00-00" +include( "RecExCond/RecExCommon_flags.py" ) + +#DetFlags.detdescr.ID_setOff() +#DetFlags.detdescr.Calo_setOff() +#DetFlags.detdescr.LAr_setOff() +#DetFlags.detdescr.Tile_setOff() +#DetFlags.detdescr.Muon_setOff() +DetFlags.detdescr.all_setOff() +DetFlags.detdescr.Calo_setOn() + +# set up all detector description +include ("RecExCond/AllDet_detDescr.py") +include( "IOVDbSvc/IOVRecExCommon.py" ) + +#globalflags.ConditionsTag = 'OFLCOND-SIM-00-00-03' +if len(globalflags.ConditionsTag())!=0: + from IOVDbSvc.CondDB import conddb + conddb.setGlobalTag(globalflags.ConditionsTag()) + +# Number of events to be processed (default is 10) +theApp.EvtMax = 10 +##svcMgr.EventSelector.FirstEvent = 1 + +# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) +svcMgr.MessageSvc.OutputLevel = 3 +#svcMgr.MessageSvc.verboseLimit = 100000 +#svcMgr.MessageSvc.debugLimit = 100000 +#ServiceMgr.IOVSvc.OutputLevel = VERBOSE + +# to read pool Data +import AthenaPoolCnvSvc.ReadAthenaPool + +# G4 Pool input +#svcMgr.EventSelector.InputCollections = ["LFN:input.root"]; +#svcMgr.EventSelector.InputCollections = "$INFN"; + +# csc datasets +#svcMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/user/f/fledroit/csc7051/calib0_mc12.007051.singlepart_gammaTT_emecIW.simul.HITS.v12003103_tid004119._00015.pool.root"] +#svcMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/user/f/fledroit/csc7003/calib0_csc11.007003.singlepart_e_Et25.simul.HITS.v12000301_tid003216._00001.pool.root"] +#svcMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/user/s/svahsen/digitization/RTT/calib1_csc11.005200.T1_McAtNlo_Jimmy.simul.HITS.v12003104_tid004131._00069.pool.root.10"] +#svcMgr.EventSelector.InputCollections = ["/afs/cern.ch/user/f/fledroit/scratch0/test_ntuples/calib0_csc11.007003.singlepart_e_Et25.simul.HITS.v12000301_tid003216._00001.pool.root"] +# +#svcMgr.EventSelector.InputCollections = ["rfio:/castor/cern.ch/grid/atlas/atlasmcdisk/valid1/HITS/valid1.007003.singlepart_e_Et25.simul.HITS.e322_s472_tid025400/HITS.025400._00002.pool.root.1"] +svcMgr.EventSelector.InputCollections = ["/afs/cern.ch/user/k/ketevi/w0/data/singleMuon100.HITS.020162._00001.pool.root"] + + + +include( "LArL1Sim/LArL1Sim_G4_jobOptions.py" ) +theLArTTL1Maker.OutputLevel = 2 + +# Configure LArTTL1Maker here +#theLArTTL1Maker.SubDetectors="LAr_All" +#theLArTTL1Maker.EmBarrelHitContainerName="LArHitEMB" +#theLArTTL1Maker.EmEndCapHitContainerName="LArHitEMEC" +#theLArTTL1Maker.HecHitContainerName="LArHitHEC" +#theLArTTL1Maker.ForWardHitContainerName="LArHitFCAL" + +#theLArTTL1Maker.EmTTL1ContainerName="LArTTL1EM" +#theLArTTL1Maker.HadTTL1ContainerName="LArTTL1HAD" + +#theLArTTL1Maker.NoiseOnOff=true + +#theLArTTL1Maker.PileUp=false +#theLArTTL1Maker.UseTriggerTime=false +#theLArTTL1Maker.TriggerTimeToolName="CosmicTriggerTimeTool" + +#theLArTTL1Maker.EmBarrelCalibrationCoeffs=[1.03, 1.024, 1.019, 1.02, 1.02, 1.024, 1.03, 1.046, 1.06, 1.053, 1.057, 1.062, 1.063, 1.076, 1.176] +#theLArTTL1Maker.EmEndCapCalibrationCoeffs=[1.176, 1.061, 1.087, 1.015, 1.019, 1.014, 1.014, 1.009, 1.01, 1.003, 1.016, 1.003, 0.993, 1.005, 0.963] +#theLArTTL1Maker.HECCalibrationCoeffs=[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.] +#theLArTTL1Maker.EmFcalCalibrationCoeffs=[1., 1., 1., 1.] +#theLArTTL1Maker.HadFcalCalibrationCoeffs=[1., 1., 1., 1.] + +#theLArTTL1Maker.NoEmCalibrationMode=false +#theLArTTL1Maker.NoHadCalibrationMode=false +#theLArTTL1Maker.ChronoTest=false +#theLArTTL1Maker.DebugThreshold=5000. + diff --git a/LArCalorimeter/LArL1Sim/share/LArSCSimpleMaker_jobOptions.py b/LArCalorimeter/LArL1Sim/share/LArSCSimpleMaker_jobOptions.py new file mode 100644 index 00000000000..6781d5e7a47 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/LArSCSimpleMaker_jobOptions.py @@ -0,0 +1,2 @@ +from LArL1Sim.LArL1SimConf import LArSCSimpleMaker +topSequence+=LArSCSimpleMaker() diff --git a/LArCalorimeter/LArL1Sim/share/SimpleSC_From_ESD.py b/LArCalorimeter/LArL1Sim/share/SimpleSC_From_ESD.py new file mode 100755 index 00000000000..c9c6ddee50e --- /dev/null +++ b/LArCalorimeter/LArL1Sim/share/SimpleSC_From_ESD.py @@ -0,0 +1,191 @@ + +EvtMax=10 +doRawData = False + + +from RecExConfig.RecFlags import rec +from AthenaCommon.BeamFlags import jobproperties +from AthenaCommon.AthenaCommonFlags import athenaCommonFlags +from AthenaCommon.GlobalFlags import globalflags +from D3PDMakerCoreComps.D3PDObject import D3PDObject + +bdir = "/afs/cern.ch/user/h/hma/work/public/data/" +PoolESDFile =bdir + "ESD.01204497._005853.pool.root.1" + + +athenaCommonFlags.PoolESDInput.set_Value_and_Lock([PoolESDFile]) +rec.readESD.set_Value_and_Lock(True) + +tuple_name = "caloD3PD_from_esd.root" + +# +#from MuonCnvExample.MuonCnvFlags import muonCnvFlags +#muonCnvFlags.RpcCablingMode="new" +# + +rec.AutoConfiguration=['everything'] +athenaCommonFlags.EvtMax.set_Value_and_Lock(EvtMax) + +rec.doHist.set_Value_and_Lock(False) +rec.doCBNT.set_Value_and_Lock(False) +rec.doWriteTAGCOM.set_Value_and_Lock(False) +rec.doWriteTAG.set_Value_and_Lock(False) +rec.doWriteAOD.set_Value_and_Lock(False) +rec.doWriteESD.set_Value_and_Lock(False) +rec.doAOD.set_Value_and_Lock(False) +rec.doMonitoring.set_Value_and_Lock(False) +rec.readAOD.set_Value_and_Lock(False) + +# RecExCommon +include ("RecExCommon/RecExCommon_topOptions.py") + +# D3PDMaker calo block + +from D3PDMakerCoreComps.MakerAlg import * +from CaloD3PDMaker.ClusterD3PDObject import * +from CaloD3PDMaker.CaloCellD3PDObject import * +from CaloD3PDMaker.CaloInfoD3PDObject import * +from CaloD3PDMaker.LArDigitD3PDObject import * +from CaloD3PDMaker.TileDigitD3PDObject import * +from EventCommonD3PDMaker.EventInfoD3PDObject import * +#from CaloD3PDMaker.TileRawChannelD3PDObject import * + +from LArL1Sim.LArL1SimConf import LArSCSimpleMaker +topSequence+=LArSCSimpleMaker() + +print 'CHECK topSequence' +print topSequence + +include( "LArDetDescr/LArDetDescr_joboptions.py" ) +include( "CaloConditions/CaloConditions_jobOptions.py" ) + +alg = MakerAlg("caloD3PD", topSequence, file = tuple_name , D3PDSvc = 'D3PD::RootD3PDSvc') +alg += EventInfoD3PDObject (10) + +alg += AllCaloCellD3PDObject (10) +#alg += SelCaloCellD3PDObject (10) +# alg += EMCellD3PDObject (10) +# alg += HECCellD3PDObject (10) +# alg += FCALCellD3PDObject (10) +# alg += TileCellD3PDObject (10) + +#alg += CaloInfoD3PDObject (10) + +#alg += ClusterD3PDObject (10) + +#alg += LArDigitD3PDObject (2) + +alg += AllCaloCellD3PDObject (1,sgkey='SCellContainer',prefix='scells_') +alg.scells_Filler.scells_Filler_Detail1.SaveId=True + +import CaloD3PDMaker +from D3PDMakerCoreComps.D3PDObject import make_SGDataVector_D3PDObject + +if doRawData : + + LArDigitD3PDObject = make_SGDataVector_D3PDObject( "LArDigitContainer", + "LArDigitContainer_MC", + "lardigit_", "LArDigitD3PDObject" ) + + LArDigitD3PDObject.defineBlock( 0, 'Digits', + CaloD3PDMaker.LArDigitFillerTool, + SaveDigit= True, + SaveId = True, + SaveSCAAddress= False, + DumpIterResults= False ) + alg+=LArDigitD3PDObject(0) + + + TileDigitD3PDObject = D3PDObject (makeTileDigitD3PDObject, 'tiledigit_', 'TileDigitD3PDObject') + + TileDigitD3PDObject.defineBlock (0, 'Digits', + CaloD3PDMaker.TileDigitFillerTool, + SaveOfflineInfo= True, + SaveHardwareInfo=True, + ) + alg+=TileDigitD3PDObject(0) + + + import CaloD3PDMaker + import D3PDMakerCoreComps + from D3PDMakerCoreComps.D3PDObject import D3PDObject + + TileRawChannelSGKey='TileRawChannelCnt' + + def makeTileRawChannelD3PDObject (name, prefix, object_name='TileRawChannelD3PDObject', getter = None, + sgkey = None, + label = None): + if sgkey == None: sgkey = TileRawChannelSGKey + if label == None: label = prefix + + if not getter: + getter = CaloD3PDMaker.SGTileRawChannelGetterTool \ + (name + '_Getter', + TypeName = 'TileRawChannelContainer', + SGKey = sgkey, + Label = label) + + # create the selected cells + from D3PDMakerConfig.D3PDMakerFlags import D3PDMakerFlags + return D3PDMakerCoreComps.VectorFillerTool (name, + Prefix = prefix, + Getter = getter, + ObjectName = object_name, + SaveMetadata = \ + D3PDMakerFlags.SaveObjectMetadata()) + + TileRawChannelD3PDObject = D3PDObject( makeTileRawChannelD3PDObject, 'tileraw_', + 'TileRawChannelD3PDObject' ) + + TileRawChannelD3PDObject.defineBlock (0, 'RawChannel', + CaloD3PDMaker.TileRawChannelFillerTool, + SaveHardwareInfo=True, + SaveRawChannel= True, + ) + + + TileRawChannelD3PDObject.defineBlock (1, 'Hardware', + CaloD3PDMaker.TileRawChannelFillerTool, + SaveHardwareInfo=True, + SaveRawChannel= False, + ) + + alg+=TileRawChannelD3PDObject(0) + + + +#from CaloD3PDMaker.LArRawChannelD3PDObject import * +#alg+=LArRawChannelD3PDObject(0) + +from TrigCaloD3PDMaker.TrigCaloD3PDMakerConf \ + import D3PD__TriggerTowerFillerTool, D3PD__TriggerTowerFillerTool + +TriggerTowerCollection_sgkey = 'TriggerTowers' +TriggerTowerD3PDObject = make_SGDataVector_D3PDObject( + 'DataVector<LVL1::TriggerTower>', + TriggerTowerCollection_sgkey, + 'tt_', + 'TriggerTowerD3PDObject') + +# Level 0 +TriggerTowerD3PDObject.defineBlock(0, 'Basics', + D3PD__TriggerTowerFillerTool) +alg+=TriggerTowerD3PDObject(0) + + + +import CaloD3PDMaker +import D3PDMakerCoreComps +from D3PDMakerCoreComps.D3PDObject import D3PDObject + + + +#if readRaw : +# # put OF iteration results on SG +# ToolSvc.LArRawChannelBuilderToolOFCIter.StoreTiming=True + +print 'CHECK : svcMgr.DetDescrCnvSvc : ' +print svcMgr.DetDescrCnvSvc +#svcMgr.StoreGateSvc.Dump=True +svcMgr.MessageSvc.OutputLevel = INFO + diff --git a/LArCalorimeter/LArL1Sim/src/LArSCL1Maker.cxx b/LArCalorimeter/LArL1Sim/src/LArSCL1Maker.cxx new file mode 100755 index 00000000000..5d00ba9355f --- /dev/null +++ b/LArCalorimeter/LArL1Sim/src/LArSCL1Maker.cxx @@ -0,0 +1,864 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// +======================================================================+ +// + + +// + Author ........: Denis O. Damazio + +// + Institut ......: BNL + +// + Creation date .: 18/11/2013 + +// + + +// +======================================================================+ +// +// ........ includes +// +#include "LArL1Sim/LArSCL1Maker.h" +// .......... utilities +// +#include <math.h> +#include <fstream> +#include "AtlasCLHEP_RandomGenerators/RandGaussZiggurat.h" +#include "CLHEP/Random/RandomEngine.h" +#include "AthenaKernel/IAtRndmGenSvc.h" +#include "GaudiKernel/ServiceHandle.h" + +#include "LArRawEvent/LArDigit.h" +#include "LArRawEvent/LArDigitContainer.h" +#include "LArSimEvent/LArHitContainer.h" +#include "LArDigitization/LArHitList.h" + +#include "CaloIdentifier/CaloIdManager.h" +#include "CaloIdentifier/LArID.h" +#include "CaloIdentifier/CaloID_Exception.h" +#include "CaloIdentifier/CaloLVL1_ID.h" +#include "LArIdentifier/LArIdManager.h" +#include "LArIdentifier/LArOnline_SuperCellID.h" +#include "LArTools/LArCablingService.h" +#include "LArTools/LArSuperCellCablingTool.h" +#include "CaloIdentifier/CaloCell_SuperCell_ID.h" +#include "CaloTriggerTool/ICaloSuperCellIDTool.h" +// +// ........ Event Header Files: +// +//#include "EventInfo/EventID.h" +//#include "EventInfo/EventInfo.h" +#include "EventInfo/PileUpTimeEventIndex.h" +// +// ........ Gaudi needed includes +// +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/Property.h" +#include "GaudiKernel/IService.h" +#include "GaudiKernel/Algorithm.h" +#include "GaudiKernel/IChronoStatSvc.h" +#include "GaudiKernel/IToolSvc.h" +#include "GaudiKernel/ListItem.h" +#include "GaudiKernel/IIncidentSvc.h" +#include "StoreGate/StoreGateSvc.h" +#include "PathResolver/PathResolver.h" + +#include "LArElecCalib/ILArfSampl.h" +#include "LArElecCalib/ILArShape.h" +#include "LArElecCalib/ILArPedestal.h" +#include "LArElecCalib/ILArNoise.h" +#include "LArElecCalib/ILArAutoCorrNoiseTool.h" +#include "LArElecCalib/ILArADC2MeVTool.h" + +// Pile up +#include "PileUpTools/PileUpMergeSvc.h" +// trigger time +#include "AthenaKernel/ITriggerTime.h" + +// + + +using CLHEP::RandGaussZiggurat; + + +LArSCL1Maker::LArSCL1Maker(const std::string& name, ISvcLocator* pSvcLocator) : + Algorithm(name, pSvcLocator) + , m_storeGateSvc("StoreGateSvc",name) + , m_atRndmGenSvc("AtRndmGenSvc",name) + , m_rndmEngineName("LArSCL1Maker") + , m_rndmEngine(0) + , m_hitmap(0) + , m_autoCorrNoiseTool("LArAutoCorrNoiseTool") + , m_adc2mevTool("LArADC2MeVTool") + , m_fSamplKey("LARfSamplSC") + , m_shapesKey("LArShapeSC") + , m_noiseKey("LArNoiseSC") + , m_pedestalKey("LArPedestalSC") +// + -------------------------------------------------------------------- + +// + Author ........: Denis O. Damazio + +// + Creation date .: 18/11/2013 + +// + Subject: SCL1 Maker constructor + +// + -------------------------------------------------------------------- + +{ + m_first=true; +// +// ........ default values of private data +// + m_detectorStore = 0; + //m_storeGateSvc = 0; + m_chronSvc = 0; + m_mergeSvc = 0; + m_useTriggerTime = false; + m_triggerTimeToolName = "CosmicTriggerTimeTool"; + p_triggerTimeTool = 0; + + m_BeginRunPriority = 100; + + m_cablingSCSvc = 0; + //m_ttSvc = 0; + m_scHelper = 0; + + m_SubDetectors = "LAr_All"; + m_SCL1ContainerName = "LArDigitSCL1"; + + + m_NoiseOnOff = true; + m_PileUp = false; + m_noEmCalibMode = false; + m_noHadCalibMode = false; + m_chronoTest = false; + m_debugThresh = 5000.; + + // + // ........ declare the private data as properties + // + + declareProperty("EventStore",m_storeGateSvc,"StoreGate Service"); + declareProperty("SubDetectors",m_SubDetectors); + declareProperty("RndmSvc", m_atRndmGenSvc); + declareProperty("AutoCorrNoiseTool",m_autoCorrNoiseTool,"Tool handle for electronic noise covariance"); + declareProperty("ADC2MeVTool",m_adc2mevTool,"Tool handle for ADC2MeV Conversion"); + + + declareProperty("SCL1ContainerName",m_SCL1ContainerName); + + declareProperty("NoiseOnOff",m_NoiseOnOff); + + declareProperty("PileUp",m_PileUp); + declareProperty("UseTriggerTime",m_useTriggerTime); + declareProperty("TriggerTimeToolName",m_triggerTimeToolName); + + declareProperty("FirstSample",m_firstSample=-1); + declareProperty("NSamples",m_nSamples=7); + + declareProperty("ChronoTest",m_chronoTest); + declareProperty("DebugThreshold",m_debugThresh); + +// +return; +} + + +LArSCL1Maker::~LArSCL1Maker() +{ + /** SCL1 Maker destructor */ + +return; +} + + +StatusCode LArSCL1Maker::initialize() +{ +// +======================================================================+ +// + + +// + Author ........: Denis O. Damazio + +// + Creation date .: 18/11/2013 + +// + + +// +======================================================================+ +// +// ......... declaration +// + MsgStream msglog(messageService(),name()); + int outputLevel = msgSvc()->outputLevel( name() ); + + m_chronSvc = chronoSvc(); + + msglog << MSG::INFO + << "***********************************************" + << endreq; + msglog << MSG::INFO + << "* Steering options for LArSCL1Maker algorithm *" + << endreq; + msglog << MSG::INFO + << "***********************************************" + << endreq; + // + // ......... print the noise flag + // + if ( m_NoiseOnOff ) + { + msglog << MSG::INFO + << "Electronic noise will be added in each SC for selected sub-detectors." + << endreq; + } + else + { + msglog << MSG::INFO + << "No electronic noise added." + << endreq; + } + +// +// ......... print the pile-up flag +// + if (m_PileUp) + { + msglog << MSG::INFO + << "take events from PileUp service" + << endreq; + } + else + { + msglog << MSG::INFO + << "no pile up" + << endreq; + } + +// +// ......... print the trigger time flag +// + if (m_useTriggerTime) + { + msglog << MSG::INFO + << "use Trigger Time service " << m_triggerTimeToolName + << endreq; + } + else + { + msglog << MSG::INFO + << "no Trigger Time used" + << endreq; + } + +// +// ....... Access Event StoreGate +// + //StatusCode sc = service ( "StoreGateSvc" , m_storeGateSvc ) ; + StatusCode sc = m_storeGateSvc.retrieve(); + if (sc.isFailure()) + { + msglog << MSG::ERROR + << "Unable to access pointer to StoreGate" + << endreq; + return StatusCode::FAILURE; + } + +// +// locate the PileUpMergeSvc and initialize our local ptr +// + if (m_PileUp) { + if (!(service("PileUpMergeSvc", m_mergeSvc)).isSuccess() || 0 == m_mergeSvc) { + msglog << MSG::ERROR + << "Could not find PileUpMergeSvc" + << endreq; + return StatusCode::FAILURE; + } + else + { + msglog << MSG::DEBUG << "PileUpMergeSvc successfully initialized" << endreq; + } + } + +// +// .........retrieve tool computing trigger time if requested +// + if (m_useTriggerTime && m_PileUp) { + msglog << MSG::INFO + << " In case of pileup, the trigger time subtraction is done in PileUpSvc " + << endreq; + msglog << MSG::INFO + << " => LArSCL1Maker will not apply Trigger Time " << endreq; + m_useTriggerTime = false; + } + + if (m_useTriggerTime) { + IToolSvc* p_toolSvc = 0; + if (!(service("ToolSvc", p_toolSvc)).isSuccess() || 0 == p_toolSvc) { + msglog << MSG::ERROR + << "Could not find ToolSvc" + << endreq; + return StatusCode::FAILURE; + } + + IAlgTool* algtool(0); + ListItem theTool(m_triggerTimeToolName.value()); + sc = p_toolSvc->retrieveTool(theTool.type(), theTool.name(),algtool); + if (sc.isFailure()) { + msglog << MSG::ERROR + << "Unable to find tool for " << m_triggerTimeToolName.value() << endreq; + p_triggerTimeTool = 0; + } + else { + p_triggerTimeTool=dynamic_cast<ITriggerTime*>(algtool); + msglog << MSG::DEBUG << "retrieved TriggerTime tool: " + << m_triggerTimeToolName.value() << endreq; + } + } + + // + // ..... need LAr and CaloIdManager to retrieve all needed helpers + // + + const CaloIdManager* caloMgr; + const LArIdManager* larMgr; + // + // ....... Access Detector StoreGate + // + sc = service( "DetectorStore", m_detectorStore ); + if ( sc.isSuccess( ) ) { + sc = m_detectorStore->retrieve(caloMgr); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve CaloIdManager from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved CaloIdManager from DetectorStore" << endreq; + } + } + sc = m_detectorStore->retrieve(larMgr); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve LArIdManager from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved LArIdManager from DetectorStore" << endreq; + } + } + + sc = m_detectorStore->retrieve(m_OnlSCHelper); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve LArOnline_SuperCellID from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved LArOnline_SuperCellID from DetectorStore" << endreq; + } + } + + } else { + msglog << MSG::ERROR << "Could not locate DetectorStore" << endreq; + return StatusCode::FAILURE; + } + + // + //..... need of course the LVL1 helper + // + m_scHelper = caloMgr->getCaloCell_SuperCell_ID(); + if (!m_scHelper) { + msglog << MSG::ERROR << "Could not access CaloCell_SuperCell_ID helper" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully accessed CaloCell_SuperCell_ID helper" << endreq; + } + } + + // ..... need cabling services, to get channels associated to each SC + // + IToolSvc* toolSvc; + StatusCode status = service( "ToolSvc",toolSvc ); + if(status.isSuccess()) { + sc = toolSvc->retrieveTool("LArSuperCellCablingTool",m_cablingSCSvc); + if(sc.isFailure()) { + msglog << MSG::ERROR << "Could not retrieve LArSuperCellCablingTool"<< endreq; + return(StatusCode::FAILURE); + } + sc = toolSvc->retrieveTool("CaloSuperCellIDTool",m_scidtool); + if(sc.isFailure()) { + msglog << MSG::ERROR << "Could not retrieve CaloTriggerTowerService"<< endreq; + return(StatusCode::FAILURE); + } + } else { + msglog << MSG::ERROR << "Could not get ToolSvc"<< endreq; + return(StatusCode::FAILURE); + } + + // Incident Service: + IIncidentSvc* incSvc; + sc = service("IncidentSvc", incSvc); + if (sc.isFailure()) + { + msglog << MSG::ERROR + << "Unable to retrieve pointer to DetectorStore " + << endreq; + return StatusCode::FAILURE; + } + //start listening to "BeginRun" + incSvc->addListener(this, "BeginRun", m_BeginRunPriority); + + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Initialization completed successfully" + << endreq; + } + + sc = m_atRndmGenSvc.retrieve(); + if(sc.isFailure()) { + msglog << MSG::ERROR << "Could not initialize random number service." << endreq; + return sc; + } + m_rndmEngine = m_atRndmGenSvc->GetEngine(m_rndmEngineName); + if(!m_rndmEngine) { + msglog << MSG::ERROR << "Could not find RndmEngine : " << m_rndmEngineName << endreq; + return StatusCode::FAILURE ; + } + + if ( m_autoCorrNoiseTool.retrieve().isFailure() ) { + msglog << MSG::ERROR << "Unable to find LArAutoCorrNoiseTool" << endreq; + return StatusCode::FAILURE; + } + + if ( m_adc2mevTool.retrieve().isFailure() ) { + msglog << MSG::ERROR << "Unable to find LArADC2MeVTool" << endreq; + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; + +} + + + +void LArSCL1Maker::handle(const Incident& /* inc*/ ) +{ + MsgStream msglog( msgSvc(), name() ); + msglog << MSG::DEBUG << "LArSCL1Maker handle()" << endreq; + + StatusCode sc = this->updateConditions(); + if (sc.isFailure()) { + msglog << MSG::ERROR << " Error from updateConditions " << endreq; + } + + + return; +} + +StatusCode LArSCL1Maker::execute() +{ +// +======================================================================+ +// + + +// + Author: Denis O. Damazio + +// + Creation date: 2003/01/13 + +// + Subject: Make the SCL1s and put them into the container + +// + + +// +======================================================================+ + // + // ......... declarations + // + MsgStream msglog(messageService(),name()); + int outputLevel = msgSvc()->outputLevel( name() ); + + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Begining of execution" + << endreq; + } + + // + // ....... fill the LArHitEMap + // + if(m_chronoTest) { + m_chronSvc->chronoStart( "fill LArHitEMap " ); + } + + int totHit=0; + if ( this->fillEMap(totHit) == StatusCode::FAILURE ) return StatusCode::FAILURE; + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "total number of hits with |E|> 1.e-06 found = " << totHit << endreq; + } + + + if(m_chronoTest) { + m_chronSvc->chronoStop ( "fill LArHitEMap " ); + m_chronSvc->chronoPrint( "fill LArHitEMap " ); + } + + if ( totHit==0) { + msglog << MSG::WARNING << " No LAr hit in the event " << endreq; + msglog << MSG::WARNING << "cannot process this event" << endreq; + return StatusCode::SUCCESS; + } + + // + // .....get the trigger time if requested + // + double trigtime=0; + if (m_useTriggerTime && p_triggerTimeTool) { + trigtime = p_triggerTimeTool->time(); + } + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Trigger time used : " << trigtime << endreq; + } + + + LArDigitContainer *scContainer = new LArDigitContainer(); + if ( scContainer == 0 ) { + msglog << MSG::ERROR + << "Could not allocate a new LArSCDigitContainer" + << endreq; + return StatusCode::FAILURE; + } + + + // ...... register the TTL1 containers into the TES + // + StatusCode sc = m_storeGateSvc->record( scContainer , m_SCL1ContainerName) ; + if( sc.isFailure() ) + { + msglog << MSG::ERROR + << "Could not record new LArDigitContainer in TES : " + << m_SCL1ContainerName + << endreq; + return StatusCode::FAILURE; + } + + // + // ... initialise vectors for sums of energy in each TT + // + unsigned int nbSC = (unsigned int)m_scHelper->calo_cell_hash_max() ; + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Max number of LAr Super-Cell= " << nbSC + << endreq; + } + std::vector<std::vector<float> > sumEnergy ; // inner index = time slot (from 0 to visEvecSize-1) + std::vector<std::vector<float> > sumTime ; // inner index = time slot (from 0 to visEvecSize-1) + sumEnergy.resize(nbSC); + sumTime.resize(nbSC); + std::vector<float> scSumE; + int scSumEvecSize=5; + + scSumE.reserve(scSumEvecSize); + std::vector< std::vector < float> > scFloatContainerTmp; + + int it = 0; + int it_end = m_hitmap->GetNbCells(); + scContainer->reserve( nbSC ); //container ordered by hash + scFloatContainerTmp.assign( nbSC, std::vector<float>(0) ); //container ordered by hash + std::vector<bool> alreadyThere; + alreadyThere.resize( nbSC ); + alreadyThere.assign( nbSC, false ); + + std::vector<HWIdentifier> hwid; + hwid.resize( nbSC ); + hwid.assign( nbSC, HWIdentifier(0) ); + + //m_nSamples=5; + std::vector<float> samples; + samples.resize( m_nSamples ); + std::vector<short> samplesInt; + samplesInt.resize( m_nSamples ); + CaloGain::CaloGain scGain = CaloGain::LARHIGHGAIN; + + for( ; it!=it_end;++it) { + LArHitList * hitlist = m_hitmap->GetCell(it); + if ( hitlist != (LArHitList*) NULL ) { + + const std::vector<std::pair<float,float> >* timeE = hitlist->getData(); + if (timeE->size() > 0 ) { + Identifier cellId = hitlist->getIdentifier(); + + + Identifier scId = m_scidtool->offlineToSuperCellID(cellId); + IdentifierHash scHash = m_scHelper->calo_cell_hash(scId) ; + if ( scHash.value() == 999999 ) continue; + HWIdentifier hwSC = m_cablingSCSvc->createSignalChannelID(scId); + IdentifierHash scHWHash = m_OnlSCHelper->channel_Hash(hwSC); + hwid[scHWHash] = hwSC; + + float factor = 1.0; + if ( m_first ) { + printConditions(hwSC); + } + factor = (m_adc2mevTool->ADC2MEV(hwSC,(CaloGain::CaloGain)1))[1] ; + factor = 1.0/m_fracS->FSAMPL(hwSC)/factor; + + ConvertHits2Samples(hwSC, scGain, timeE, samples); + std::vector< float >& vec = scFloatContainerTmp.at(scHWHash); + if ( !alreadyThere[scHWHash] ){ + alreadyThere[scHWHash]=true; + for(unsigned int i=0; i< samples.size(); i++) + vec.push_back( factor * samples[i] ); + } else { + for(unsigned int i=0; i< samples.size(); i++) + vec[i]+= ( factor * samples[i] ); + } + } + + + } // if hit list is not empty + } // it end + + it=0; + it_end=nbSC; + int cc=0; + short MAXADC=4096; // 12 bits ADC? + std::vector<float> noise(m_nSamples); + double Rndm[32]; + for( ; it != it_end; ++it){ + + if ( alreadyThere[it] ){ + const std::vector< float >& vec= scFloatContainerTmp.at(it); + HWIdentifier id = hwid[it]; + // reset noise + noise.assign(m_nSamples,0); + + // noise definition + if ( m_NoiseOnOff ) { + float SigmaNoise = (m_NoiseSC->noise(id,0)); + int index; + const std::vector<float>* CorrGen = &(m_autoCorrNoiseTool->autoCorrSqrt(id,0,m_nSamples)); + + RandGaussZiggurat::shootArray(m_rndmEngine,m_nSamples,Rndm,0.,1.); + + for(int i=0;i<(int)m_nSamples;i++){ + noise[i]=0.; + for(int j=0;j<=i;j++){ + index = i* m_nSamples + j; + noise[i] += Rndm[j] * (*CorrGen)[index]; + } + noise[i]=noise[i]*SigmaNoise; + } + } + + int ped = m_PedestalSC->pedestal(id,0); + samplesInt.assign( m_nSamples, 0 ); + for(unsigned int i=0; i< vec.size(); i++) { + samplesInt[i]=rint(vec[i]+ped+noise[i]); + if ( samplesInt[i] >= MAXADC ) samplesInt[i]=MAXADC-1; + if ( samplesInt[i] < 0 ) samplesInt[i]=0; + } + LArDigit* dig = new LArDigit(id, scGain, samplesInt ); + scContainer->push_back(dig); + } else { + cc++; + } + } + + if(m_chronoTest) { + m_chronSvc->chronoStop( "LArSCL1Mk hit loop " ); + m_chronSvc->chronoPrint( "LArSCL1Mk hit loop " ); + m_chronSvc->chronoStart( "LArSCL1Mk SC loop " ); + } + + m_first = false; + return StatusCode::SUCCESS; +} + + +StatusCode LArSCL1Maker::finalize() +{ +// +======================================================================+ +// + + +// + Author: Denis O. Damazio + +// + Creation date: 18/11/2013 + +// + + +// +======================================================================+ +// +// ......... declaration +// + MsgStream msglog(messageService(),name()); + + msglog << MSG::INFO << " LArSCL1Maker finalize completed successfully" + << endreq; + + m_chronSvc->chronoPrint( "LArSCL1Mk hit loop " ); + m_chronSvc->chronoPrint( "LArSCL1Mk TT loop " ); + + return StatusCode::SUCCESS; + +} + + +StatusCode LArSCL1Maker::fillEMap(int& totHit) + +{ +// +======================================================================+ +// + + +// + Author: Denis O. Damazio + +// + Creation date: 18/11/2013 + +// + + +// +======================================================================+ + + + MsgStream msglog(messageService(),name()); +// +// ........ reset the map Energies +// + StatusCode sc = m_detectorStore->retrieve(m_hitmap,"LArHitEMap"); + if (sc.isFailure()) { + msglog << " cannot retrieve LArHitEMap from detector Store " << endreq; + return sc; + } + totHit = m_hitmap->GetNbCells(); + + return StatusCode::SUCCESS; +} + + +//=========================================================================== +int LArSCL1Maker::decodeInverse(int region, int eta) +{ +//=========================================================================== +// Maps [ region , eta ] to [ Ieta] +// ========================================================================== +// Problem: this is NOT a bijection, because of the "barrel end" zone. +// Convention: only the barrel part of the identifying fields is returned +//=========================================================================== + + int Ieta=0; + if(region == 0){ // Barrel + EC-OW + if(eta <= 14) { + Ieta=eta+1; + } + else { + Ieta=eta-13; + } + } + else if( region == 1 || region == 2) { // EC-IW + if(region == 1) { + Ieta=eta+12; + } + else { + Ieta=15; + } + } + else if(region == 3) { // FCAL + Ieta=eta+1; + } + return Ieta ; +} + + +StatusCode LArSCL1Maker::updateConditions(){ + + MsgStream msglog(messageService(),name()); + msglog << MSG::DEBUG << "Updating conditions" << endreq; + + if (m_detectorStore == 0 ) { + msglog << "No valid DetStore, forget this" << endreq; + return StatusCode::FAILURE; + } + + if ( (m_detectorStore->retrieve(m_fracS,m_fSamplKey)).isFailure() && !m_fracS ) { + msglog << "could not get fSampl conditions with key " << m_fSamplKey << endreq; + return StatusCode::FAILURE; + } + + if ( (m_detectorStore->retrieve(m_shapes,m_shapesKey)).isFailure() && !m_shapes ) { + msglog << "could not get shape conditions with key " << m_shapesKey << endreq; + return StatusCode::FAILURE; + } + + if ( (m_detectorStore->retrieve(m_NoiseSC,m_noiseKey)).isFailure() && !m_NoiseSC ) { + msglog << "could not get noise conditions with key " << m_noiseKey << endreq; + return StatusCode::FAILURE; + } + + if ( (m_detectorStore->retrieve(m_PedestalSC,m_pedestalKey)).isFailure() && !m_PedestalSC ) { + msglog << "could not get pedestal conditions with key " << m_pedestalKey << endreq; + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; + +} + +void LArSCL1Maker::printConditions(const HWIdentifier& hwSC){ + + MsgStream msglog(messageService(),name()); + msglog << MSG::INFO << "HW Identifier : " << hwSC.get_identifier32().get_compact() << endreq; + if ( m_shapes ) { + ILArShape::ShapeRef_t shape = m_shapes->Shape(hwSC,0); + ILArShape::ShapeRef_t shapeder = m_shapes->ShapeDer(hwSC,0); + msglog << MSG::INFO << "shape0.size() : " << shape.size() << endreq; + for(unsigned int i=0;i<shape.size();i++) + msglog << "shape[" << i << "]=" << shape[i] << " - " << shapeder[i] << "; "; + msglog << endreq; + } + if ( m_fracS ) { + msglog << MSG::INFO << "fSample : " << m_fracS->FSAMPL(hwSC) << endreq; + } + if ( m_adc2mevTool ) { + msglog << MSG::INFO << "Ramp (gain0) : " << (m_adc2mevTool->ADC2MEV(hwSC,(CaloGain::CaloGain)0))[1] << endreq; + msglog << MSG::INFO << "Ramp (gain1) : " << (m_adc2mevTool->ADC2MEV(hwSC,(CaloGain::CaloGain)1))[1] << endreq; + msglog << MSG::INFO << "Ramp (gain2) : " << (m_adc2mevTool->ADC2MEV(hwSC,(CaloGain::CaloGain)2))[1] << endreq; + } + if ( m_PedestalSC ) { + msglog << MSG::INFO << "Pedestal : " << m_PedestalSC->pedestal(hwSC,0) << endreq; + } + if ( m_NoiseSC ) { + msglog << "Noise : " << m_NoiseSC->noise(hwSC,0) << endreq; + } + if ( m_autoCorrNoiseTool ) { + const std::vector<float>* CorrGen = &(m_autoCorrNoiseTool->autoCorrSqrt(hwSC,0,m_nSamples)); + msglog << MSG::INFO << "Auto : " ; + for(size_t ii=0;ii<m_nSamples;++ii){ + msglog << CorrGen->at(ii) << " "; + } + msglog << endreq; + } +} + +void LArSCL1Maker::ConvertHits2Samples(const HWIdentifier & hwSC, CaloGain::CaloGain igain, + const std::vector<std::pair<float,float> > *TimeE, std::vector<float>& samples) +{ +// Converts hits of a particular LAr cell into energy samples +// declarations + + int nsamples ; + int nsamples_der ; + int i ; + int j ; + float energy ; + float time ; + +// ........ retrieve data (1/2) ................................ +// + ILArShape::ShapeRef_t Shape = m_shapes->Shape(hwSC,igain); + ILArShape::ShapeRef_t ShapeDer = m_shapes->ShapeDer(hwSC,igain); + + nsamples = Shape.size(); + nsamples_der = ShapeDer.size(); + + if (nsamples==0) { + std::cout << " No samples for cell = " << hwSC << std::endl; + return; + } + + std::vector<std::pair<float,float> >::const_iterator first = TimeE->begin(); + std::vector<std::pair<float,float> >::const_iterator last = TimeE->end(); + samples.clear(); + samples.assign(m_nSamples,0); + //m_firstSample=0; + + while (first != last) + { + energy = (*first).first; + time = (*first).second; + + // Atlas like mode where we use 25ns binned pulse shape and derivative to deal with time offsets +// shift between reference shape and this time + int ishift=(int)(rint(time/25.)); + double dtime=time-25.*((double)(ishift)); + + for (i=0;i<(int)m_nSamples;i++) + { + j = i - ishift + m_firstSample; + if (j >=0 && j < nsamples ) { + if (j<nsamples_der && std::fabs(ShapeDer[j])<10. ) + samples[i] += (Shape[j]- ShapeDer[j]*dtime)*energy ; + else samples[i] += Shape[j]*energy ; + } + } + + ++first; + } // loop over hits + + return; + +} diff --git a/LArCalorimeter/LArL1Sim/src/LArSCSimpleMaker.cxx b/LArCalorimeter/LArL1Sim/src/LArSCSimpleMaker.cxx new file mode 100755 index 00000000000..5a07b0fa77e --- /dev/null +++ b/LArCalorimeter/LArL1Sim/src/LArSCSimpleMaker.cxx @@ -0,0 +1,210 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// +======================================================================+ +// + + +// + Author ........: Denis Oliveira Damazio + +// + Institute .....: BNL + +// + Creation date .: 17/06/2012 + +// + + +// +======================================================================+ +// +// ........ includes +// +#include "LArL1Sim/LArSCSimpleMaker.h" +//#include "LArRawEvent/LArRawChannelContainer.h" +#include "CaloTriggerTool/ICaloSuperCellIDTool.h" +#include "CaloEvent/CaloCellContainer.h" +#include "CaloDetDescr/CaloDetDescrManager.h" +#include "CaloIdentifier/CaloIdManager.h" +#include "CaloIdentifier/CaloCell_SuperCell_ID.h" +#include "StoreGate/StoreGateSvc.h" +#include "boost/foreach.hpp" +#include <cmath> +#include "GeoModelInterfaces/IGeoModelSvc.h" + +#include <iostream> + +/** + * @brief Standard Gaudi algorithm constructor. + */ +LArSCSimpleMaker::LArSCSimpleMaker(const std::string& name, + ISvcLocator* pSvcLocator) + : AthAlgorithm(name, pSvcLocator), + m_scidtool ("CaloSuperCellIDTool"), + m_sem_mgr(0), + m_calo_id_manager(0) +{ + declareProperty ("SCIDTool", m_scidtool, + "Offline / supercell mapping tool."); + declareProperty ("CellContainer", m_cellContainer = "AllCalo", + "SG key for the input calorimeter cell container."); + declareProperty ("SCellContainer", m_sCellContainer = "SCellContainer", + "SG key for the output supercell LAr channel container."); +} + + +/** + * @brief Standard Gaudi initialize method. + */ +StatusCode LArSCSimpleMaker::initialize() +{ + + return StatusCode::SUCCESS; +} + + +/** + * @brief Algorithm execute method. + */ +StatusCode LArSCSimpleMaker::execute() +{ + + CHECK( m_scidtool.retrieve() ); + CHECK( detStore()->retrieve (m_sem_mgr, "CaloSuperCellMgr") ); + CHECK( detStore()->retrieve (m_calo_id_manager, "CaloIdManager") ); + + const CaloCell_SuperCell_ID * calo_sc_id = m_calo_id_manager->getCaloCell_SuperCell_ID(); + const CaloCell_ID * calo_cell_id = m_calo_id_manager->getCaloCell_ID(); + + const Tile_SuperCell_ID * tile_sc_id = m_calo_id_manager->getTile_SuperCell_ID(); + const TileID * tile_cell_id = m_calo_id_manager->getTileID(); + + MsgStream msglog(messageService(),name()); + + if (m_dataPool.allocated()==0) + m_dataPool.reserve (calo_sc_id->calo_cell_hash_max()); + + const DataHandle<CaloCellContainer> cells; + CHECK( evtStore()->retrieve(cells, m_cellContainer) ); + if ( msgLvl (MSG::DEBUG) ) + msg(MSG::DEBUG) << "Got container Size : " << cells->size() << endreq; + + int hash_max = calo_sc_id->calo_cell_hash_max(); + std::vector<float> energies (hash_max,0); + std::vector<uint16_t> provenances (hash_max,0); + std::vector<uint16_t> gains (hash_max,0); + std::vector<uint16_t> qualities (hash_max,0); + + + BOOST_FOREACH (const CaloCell* cell, *cells) { + Identifier cell_id = cell->ID(); + Identifier sCellID = m_scidtool->offlineToSuperCellID (cell_id); + + if (!sCellID.is_valid()) { + // msg(MSG::WARNING) << " SC ID not valid " << sCellID.get_identifier32().get_compact() << " offline id = " << cell->ID().get_identifier32().get_compact()<< endreq; + // calo_sc_id->print(cell->ID()); + continue; + } + IdentifierHash hash; + hash = calo_sc_id->calo_cell_hash (sCellID); + assert (hash < energies.size() ); + energies[hash] += cell->energy(); + provenances[hash] |= cell->provenance(); + gains[hash] = std::max(gains[hash],(uint16_t)cell->gain()); + if ( qualities[hash] + (int) cell->quality() > 65535 ){ + qualities[hash] = 65535 ; + }else + { + qualities[hash] += cell->quality(); + } + + if ( calo_cell_id->is_tile(cell_id) && tile_cell_id->sampling(cell_id)==TileID::SAMP_D){ + // Special case for SAMP_D in tile. Signal is split into two SCs + int section = tile_cell_id->section (cell_id); + int side = tile_cell_id->side (cell_id); + int module = tile_cell_id->module (cell_id); + int tower = tile_cell_id->tower (cell_id); + + int section1 = section; + int section2 = section; + int side1 = side; + int side2 = side; + int tower1= tower; + int tower2= tower-1; + + if (tower ==0){ + side1 = -1; + side2 = 1; + tower1=0; + tower2=0; + } + if (tower==10){ + section2 = TileID::BARREL; + section1 = TileID::EXTBAR; + } + + + Identifier sc_id1 = tile_sc_id->cell_id(section1,side1,module,tower1,0); + Identifier sc_id2 = tile_sc_id->cell_id(section2,side2,module,tower2,0); + + + int hash1 = calo_sc_id->calo_cell_hash (sc_id1); + int hash2 = calo_sc_id->calo_cell_hash (sc_id2); + + energies[hash1] += cell->energy()/2.; + energies[hash2] += cell->energy()/2.; + //ignore the rest for now + + /* + if (module==0){ + + std::cout<<" cell, " <<tile_cell_id->print_to_string (cell_id)<< " " <<cell->eta()<< std::endl; + std::cout<<" sc1, " <<tile_sc_id->print_to_string(sc_id1)<<std::endl; + std::cout<<" sc2, " <<tile_sc_id->print_to_string(sc_id2)<<std::endl; + + //std::cout<<" section, side, module, tower, sampling "<<section<<" " << side<< " "<< module<< " " << tower << " " << sampling<< " " << cell->eta() << std::endl; + std::cout<<" cell, " <<tile_cell_id->print_to_string (cell_id)<< " " <<cell->eta()<< std::endl; + std::cout<<" sc, " <<tile_sc_id->print_to_string(sCellID)<<std::endl; + std::vector<Identifier> ids_in_SC = m_scidtool->superCellToOfflineID (sCellID); + BOOST_FOREACH (Identifier id_in_SC, ids_in_SC) + { + std::cout<<" cells in sc " <<tile_cell_id->print_to_string(id_in_SC)<< std::endl; + } + } + */ + } + + + } + + CaloCellContainer* superCellContainer= new CaloCellContainer(SG::VIEW_ELEMENTS); + superCellContainer->reserve(energies.size()); + if ( evtStore()->record(superCellContainer, m_sCellContainer).isFailure() ){ + REPORT_ERROR (StatusCode::FAILURE) << "Could not store container"; + delete superCellContainer; + return StatusCode::FAILURE; + } + + for (unsigned int i=0; i < energies.size(); i++) { + + const CaloDetDescrElement* dde = m_sem_mgr->get_element (i); + if (!dde) { + // msg(MSG::WARNING) << " Not valid DDE, hash index = "<< i << endreq; + continue; + } + + CaloCell* ss = m_dataPool.nextElementPtr(); + ss->setCaloDDE( m_sem_mgr->get_element (i)); + ss->setEnergy( energies[i] ); + ss->setTime( 0. ); + + ss->setQuality( qualities[i] ); + if (calo_sc_id->is_tile(ss->ID())) + { + ss->setProvenance( 0 ); + ss->setGain( (CaloGain::CaloGain) 0 ); + } + else + { + ss->setProvenance( provenances[i] ); + ss->setGain( (CaloGain::CaloGain) gains[i] ); + } + superCellContainer->push_back(ss); + + } + + return StatusCode::SUCCESS; +} + diff --git a/LArCalorimeter/LArL1Sim/src/LArTTL1Calib.cxx b/LArCalorimeter/LArL1Sim/src/LArTTL1Calib.cxx new file mode 100755 index 00000000000..93779305f89 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/src/LArTTL1Calib.cxx @@ -0,0 +1,452 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// +======================================================================+ +// + + +// + Author ........: F. Ledroit + +// + Institut ......: ISN Grenoble + +// + Creation date .: 09/01/2003 + +// + + +// +======================================================================+ +// +// ........ includes +// +#include "LArL1Sim/LArTTL1Calib.h" +// .......... utilities +// + +#include "LArRawEvent/LArTTL1.h" +#include "LArRawEvent/LArTTL1Container.h" + +#include "CaloIdentifier/CaloLVL1_ID.h" +// +// ........ Event Header Files: +// +#include "EventInfo/EventID.h" +#include "EventInfo/EventInfo.h" +// +// ........ Gaudi needed includes +// +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/Property.h" +#include "GaudiKernel/Algorithm.h" +#include "GaudiKernel/ISvcLocator.h" + +#include "GaudiKernel/INTupleSvc.h" +//#include "GaudiKernel/NTupleDirectory.h" +#include "GaudiKernel/SmartDataPtr.h" + +#include "StoreGate/StoreGateSvc.h" + +// + +LArTTL1Calib::LArTTL1Calib(const std::string& name, ISvcLocator* pSvcLocator) : + Algorithm(name, pSvcLocator) +// + -------------------------------------------------------------------- + +// + Author ........: F. Ledroit + +// + Creation date .: 20/04/2004 + +// + Subject: TTL1 Calib constructor + +// + -------------------------------------------------------------------- + +{ +// +// ........ default values of private data +// +// m_thresholdGeVTTs = 0.1; // pb here: threshold depends on whether calib already applied or not... +// m_thresholdGeVTTs = 0.001; // and whether want to study calib or noise + m_thresholdGeVTTs = 0.; + m_maxNtt = 8000; + // m_maxNtt = 500; + + m_pnt = 0; + + m_storeGateSvc = 0; + m_lvl1Helper = 0; + +// + return; +} + + +LArTTL1Calib::~LArTTL1Calib() +{ + /** TTL1 Calib destructor */ + +return; +} + + +StatusCode LArTTL1Calib::initialize() +{ +// +======================================================================+ +// + + +// + Author ........: F. Ledroit + +// + Creation date .: 09/01/2003 + +// + + +// +======================================================================+ +// +// ......... declaration +// + MsgStream msglog(messageService(),name()); + msglog << MSG::DEBUG << "in initialize()" << endreq; + + NTuple::Directory* dir = 0; + NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1"); + + if ( file1 != 0) { + if ( !(dir=ntupleSvc()->createDirectory("/NTUPLES/FILE1/DIR")) ) { + msglog << MSG::ERROR << "Cannot create directory /NTUPLES/FILE1/DIR" << endreq; + } + } + else { + msglog << MSG::ERROR << "Cannot access file /NTUPLES/FILE1" << endreq; + return StatusCode::FAILURE; + } + + msglog << MSG::DEBUG << "created directory /NTUPLES/FILE1/DIR" << endreq; + + NTuplePtr pnt(ntupleSvc(),"/NTUPLES/FILE1/DIR/10"); + if (!pnt) { + pnt = ntupleSvc()->book(dir,10 ,CLID_ColumnWiseTuple,"LArLvl1TT"); + + if(pnt){ + msglog << MSG::DEBUG << "booked ntuple LArLvl1TT" << endreq; + m_pnt = pnt; + StatusCode sc = pnt->addItem ("nTotEm",m_ntotem); + if ( sc == StatusCode::FAILURE ) { + msglog<<MSG::ERROR << " could not add item to col wise ntuple" << endreq; + return StatusCode::FAILURE; + } + + sc = pnt->addItem ("EtTotEm",m_etotem); + sc = pnt->addItem ("n3x3Em",m_n3x3em); + sc = pnt->addItem ("E3x3Em",m_e3x3em); + sc = pnt->addItem ("EmostEm",m_emostem); + sc = pnt->addItem ("BarrelECEm",m_bec_mostEem); + sc = pnt->addItem ("EmFcalEm",m_emfcal_mostEem); + sc = pnt->addItem ("IetaEm",m_Ieta_mostEem); + sc = pnt->addItem ("LphiEm",m_Lphi_mostEem); + sc = pnt->addItem ("nTTaboveThrEm",m_nttem,0,m_maxNtt); + sc = pnt->addItem ("EtTTEm",m_nttem,m_ettem); + sc = pnt->addItem ("BECTTEm",m_nttem,m_becttem); + sc = pnt->addItem ("EmFcalTTEm",m_nttem,m_emfcalttem); + sc = pnt->addItem ("IetaTTEm",m_nttem,m_Ietattem); + sc = pnt->addItem ("LphiTTEm",m_nttem,m_Lphittem); + + sc = pnt->addItem ("nTotHad",m_ntothad); + sc = pnt->addItem ("EtTotHad",m_etothad); + sc = pnt->addItem ("n3x3Had",m_n3x3had); + sc = pnt->addItem ("E3x3Had",m_e3x3had); + sc = pnt->addItem ("EmostHad",m_emosthad); + sc = pnt->addItem ("HecFcalHad",m_hecfcal_mostEhad); + sc = pnt->addItem ("IetaHad",m_Ieta_mostEhad); + sc = pnt->addItem ("LphiHad",m_Lphi_mostEhad); + sc = pnt->addItem ("nTTaboveThrHad",m_ntthad,0,m_maxNtt); + sc = pnt->addItem ("EtTTHad",m_ntthad,m_etthad); + sc = pnt->addItem ("HecFcalTTHad",m_ntthad,m_hecfcaltthad); + sc = pnt->addItem ("IetaTTHad",m_ntthad,m_Ietatthad); + sc = pnt->addItem ("LphiTTHad",m_ntthad,m_Lphitthad); + + } else { + msglog << MSG::ERROR << "could not book Ntuple" << endreq; + } + } else { + msglog << MSG::ERROR << "Ntuple is already booked" << endreq; + } +// +// ....... Access Event StoreGate +// + StatusCode sc = service ( "StoreGateSvc" , m_storeGateSvc ) ; + if (sc.isFailure()) + { + msglog << MSG::ERROR + << "Unable to access pointer to StoreGate" + << endreq; + return StatusCode::FAILURE; + } + + StoreGateSvc* detStore = 0; + const CaloLVL1_ID* lvl1Helper; + sc = service( "DetectorStore", detStore ); + if ( sc.isSuccess( ) ) { + sc = detStore->retrieve(lvl1Helper); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve CaloLVL1_ID from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + msglog << MSG::DEBUG << "Successfully retrieved CaloLVL1_ID from DetectorStore" << endreq; + } + } + else { + msglog << MSG::ERROR << "Could not locate DetectorStore" << endreq; + return StatusCode::FAILURE; + } + m_lvl1Helper = lvl1Helper; + + msglog << MSG::DEBUG + << "Initialization completed successfully" + << endreq; + + + return StatusCode::SUCCESS; + +} + + + + + +StatusCode LArTTL1Calib::execute() +{ +// +======================================================================+ +// + + +// + Author: F. Ledroit + +// + Creation date: 2003/01/13 + +// + + +// +======================================================================+ + // + // ......... declarations + // + MsgStream msglog(messageService(),name()); + + msglog << MSG::DEBUG + << "Begining of execution" + << endreq; + + LArTTL1Container::const_iterator ttl1iter; + Identifier ttOffId; + + std::vector<std::string> containerVec; // FIX ME USING PROPERTIES !!! + containerVec.push_back("LArTTL1EM"); + containerVec.push_back("LArTTL1HAD"); + + +// +// ............ loop over the wanted ttl1 containers (1 em, 1 had) +// + + for (unsigned int iTTL1Container=0;iTTL1Container<containerVec.size();iTTL1Container++) + { + // + // ..... Get the pointer to the TTL1 Container from StoreGate + // + msglog << MSG::DEBUG << " asking for: " + << containerVec[iTTL1Container] + << endreq; + + const LArTTL1Container* ttl1_container ; + StatusCode sc = m_storeGateSvc->retrieve( ttl1_container, containerVec[iTTL1Container] ) ; + + if ( sc.isFailure() || !ttl1_container) + { + msglog << MSG::ERROR << "Could not retrieve a LArTTL1Container " + << containerVec[iTTL1Container] + << endreq; + return StatusCode::FAILURE; + } + + // + msglog << MSG::DEBUG + << "number of ttl1s: " << ttl1_container->size() + << endreq; + + int samp = iTTL1Container; + int indexem=0; + int indexhad=0; + + float e_mostE=0.; + int bec_mostE=99; + int hecfcal_mostE=99; + int Ieta_mostE=99; + int Lphi_mostE=99; + int phiMax_mostE = -99; + + // + // ....... loop over ttl1s and get informations + // + + for(ttl1iter=ttl1_container->begin(); + ttl1iter != ttl1_container->end();ttl1iter++) { + // ttOnlId = (*ttl1iter)->ttOnlineID(); + ttOffId = (*ttl1iter)->ttOfflineID(); + std::vector<float> sampleV = (*ttl1iter)->samples(); + float e = sampleV[3] / 1000. ; // go to GeV + int reg = m_lvl1Helper->region(ttOffId); + int eta = m_lvl1Helper->eta(ttOffId); + int bec=99; + int hecfcal=99; + int Ieta=99; + decodeInverseTTChannel(reg, samp, eta, bec, hecfcal, Ieta); + int Lphi= m_lvl1Helper->phi(ttOffId) + 1; + int phiMax = m_lvl1Helper->phi_max(ttOffId); + if(e>e_mostE) { + e_mostE=e; + bec_mostE = bec; + hecfcal_mostE = hecfcal; + Ieta_mostE = Ieta; + Lphi_mostE = Lphi; + phiMax_mostE = phiMax; + } + + if(samp == 0) { + // EM towers + m_etotem += e; + if(indexem < m_maxNtt) { + if(e > m_thresholdGeVTTs) { + m_ettem[indexem] = e; + m_becttem[indexem] = bec; + m_emfcalttem[indexem] = hecfcal; + m_Ietattem[indexem] = Ieta; + m_Lphittem[indexem] = Lphi; + indexem++; + } + } + else { + if (m_maxNtt == indexem) { + msglog << MSG::WARNING + << "Too many em TT. Save only " << m_maxNtt + << endreq; + break; + } + } + } + else { + // HAD towers + m_etothad += e; + if(indexhad < m_maxNtt) { + if(e > m_thresholdGeVTTs) { + m_etthad[indexhad] = e; + m_hecfcaltthad[indexhad] = hecfcal; + m_Ietatthad[indexhad] = Ieta; + m_Lphitthad[indexhad] = Lphi; + indexhad++; + } + } + else { + if (m_maxNtt == indexhad) { + msglog << MSG::WARNING + << "Too many had TT. Save only " << m_maxNtt + << endreq; + break; + } + } + } + + } // .... end of loop over ittl1s + msglog << MSG::DEBUG + << "number of em/had ttl1s: " << indexem << " / " << indexhad + << endreq; + + + int etaLim=13; // barrel-ec transition + if(hecfcal_mostE >0 ) {etaLim=14;} + if(samp == 0) { + // em: + m_ntotem = ttl1_container->size() ; + m_nttem = indexem; + for (int i=0; i<indexem; ++i) { + int deltaEta = abs(Ieta_mostE - m_Ietattem[i]); + int deltaPhi = abs(Lphi_mostE - m_Lphittem[i]); // not quite 3x3 + if( (deltaEta <= 1 || deltaEta >= etaLim) // because phi granularity changes + && (deltaPhi <= 1 || deltaPhi == phiMax_mostE) ) { // at eta=2.5 and 3.2 + m_n3x3em +=1; + m_e3x3em += m_ettem[i]; + } + } + m_emostem = e_mostE; + m_bec_mostEem = bec_mostE; + m_emfcal_mostEem = hecfcal_mostE; + m_Ieta_mostEem = Ieta_mostE; + m_Lphi_mostEem = Lphi_mostE; + } + else { + // had: + m_ntothad = ttl1_container->size() ; + m_ntthad = indexhad; + for (int i=0; i<indexhad; ++i) { + int deltaEta = abs( Ieta_mostE- m_Ietatthad[i]); + int deltaPhi = abs( Lphi_mostE- m_Lphitthad[i]); + if( (deltaEta <= 1 || deltaEta >=etaLim) + && (deltaPhi <= 1 || deltaPhi == phiMax_mostE) ) { + m_n3x3had +=1; + m_e3x3had += m_etthad[i]; + } + } + m_emosthad = e_mostE; + m_hecfcal_mostEhad = hecfcal_mostE; + m_Ieta_mostEhad = Ieta_mostE; + m_Lphi_mostEhad = Lphi_mostE; + } + + } // .... end of loop over containers + msglog << MSG::DEBUG + << "end of loop over containers, write record" + << endreq; + + StatusCode sc = ntupleSvc()->writeRecord(m_pnt); + if (sc!=StatusCode::SUCCESS) { + msglog << MSG::ERROR << "writeRecord failed" << endreq; + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; + +} + + +StatusCode LArTTL1Calib::finalize() +{ +// +======================================================================+ +// + + +// + Author: F. Ledroit + +// + Creation date: 2003/01/13 + +// + + +// +======================================================================+ +// +// ......... declaration +// + MsgStream msglog(messageService(),name()); + +// + msglog << MSG::DEBUG << " LArTTL1Calib finalize completed successfully" + << endreq; + +// + return StatusCode::SUCCESS; + +} + +void +LArTTL1Calib::decodeInverseTTChannel(int region, int layer, int eta, int & bec, int & emhf, int & Ieta) const +{ + // fix me !!! this method is duplicated from LArCablingService + if(region == 0){ + if(eta <= 14) { + bec=0; // EM barrel + emhf=0; + Ieta=eta+1; + } + else { + bec=1; // EMEC or HEC + emhf=layer; + Ieta=eta-13; + } + } + else if( region == 1 || region == 2) { + bec=1; + emhf = layer ; + if(region == 1) { + Ieta=eta+12; + } + else { + Ieta=15; + } + } + else if(region == 3) { // FCAL only + bec=1; + emhf=2+layer; + Ieta=eta+1; + } +} + + diff --git a/LArCalorimeter/LArL1Sim/src/LArTTL1Maker.cxx b/LArCalorimeter/LArL1Sim/src/LArTTL1Maker.cxx new file mode 100755 index 00000000000..3af0e11e734 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/src/LArTTL1Maker.cxx @@ -0,0 +1,2649 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// +======================================================================+ +// + + +// + Author ........: F. Ledroit + +// + Institut ......: ISN Grenoble + +// + Creation date .: 09/01/2003 + +// + + +// +======================================================================+ +// +// ........ includes +// +#include "LArL1Sim/LArTTL1Maker.h" +// .......... utilities +// +#include <math.h> +#include <fstream> +#include "AtlasCLHEP_RandomGenerators/RandGaussZiggurat.h" +#include "CLHEP/Random/RandomEngine.h" +#include "AthenaKernel/IAtRndmGenSvc.h" +#include "GaudiKernel/ServiceHandle.h" + +#include "LArRawEvent/LArTTL1.h" +#include "LArRawEvent/LArTTL1Container.h" +#include "LArSimEvent/LArHitContainer.h" +#include "LArDigitization/LArHitList.h" + +#include "CaloIdentifier/CaloIdManager.h" +#include "CaloIdentifier/LArID.h" +#include "CaloIdentifier/CaloID_Exception.h" +#include "CaloIdentifier/CaloLVL1_ID.h" +#include "LArIdentifier/LArIdManager.h" +#include "LArTools/LArCablingService.h" +#include "CaloTriggerTool/CaloTriggerTowerService.h" +// +// ........ Event Header Files: +// +//#include "EventInfo/EventID.h" +//#include "EventInfo/EventInfo.h" +#include "EventInfo/PileUpTimeEventIndex.h" +// +// ........ Gaudi needed includes +// +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/Property.h" +#include "GaudiKernel/IService.h" +#include "GaudiKernel/Algorithm.h" +#include "GaudiKernel/IChronoStatSvc.h" +#include "GaudiKernel/IToolSvc.h" +#include "GaudiKernel/ListItem.h" +#include "GaudiKernel/IIncidentSvc.h" +#include "StoreGate/StoreGateSvc.h" +#include "PathResolver/PathResolver.h" + +// Pile up +#include "PileUpTools/PileUpMergeSvc.h" +// trigger time +#include "AthenaKernel/ITriggerTime.h" + +// + + +using CLHEP::RandGaussZiggurat; + + +LArTTL1Maker::LArTTL1Maker(const std::string& name, ISvcLocator* pSvcLocator) : + Algorithm(name, pSvcLocator) + , m_storeGateSvc("StoreGateSvc",name) + , m_atRndmGenSvc("AtRndmGenSvc",name) + , m_rndmEngineName("LArTTL1Maker") + , m_rndmEngine(0) + , m_hitmap(0) + , m_fSamplKey("LArfSampl") +// + -------------------------------------------------------------------- + +// + Author ........: F. Ledroit + +// + Creation date .: 09/01/2003 + +// + Subject: TTL1 Maker constructor + +// + -------------------------------------------------------------------- + +{ +// +// ........ default values of private data +// + m_detectorStore = 0; + //m_storeGateSvc = 0; + m_chronSvc = 0; + m_mergeSvc = 0; + m_useTriggerTime = false; + m_triggerTimeToolName = "CosmicTriggerTimeTool"; + p_triggerTimeTool = 0; + m_useMapfromStore = true; + + m_BeginRunPriority = 100; + + m_cablingSvc = 0; + m_ttSvc = 0; + m_lvl1Helper = 0; + m_emHelper = 0; + m_hecHelper = 0; + m_fcalHelper = 0; + + m_SubDetectors = "LAr_All"; + m_EmTTL1ContainerName = "LArTTL1EM"; + m_HadTTL1ContainerName = "LArTTL1HAD"; + + m_EmBarrelHitContainerName = "LArHitEMB"; + m_EmEndCapHitContainerName = "LArHitEMEC"; + m_HecHitContainerName = "LArHitHEC"; + m_ForWardHitContainerName = "LArHitFCAL"; + + m_NoiseOnOff = true; + m_PileUp = false; + m_noEmCalibMode = false; + m_noHadCalibMode = false; + m_chronoTest = false; + m_debugThresh = 5000.; + + m_calibCoeffEmb.resize(s_NBETABINS); + m_calibCoeffEmec.resize(s_NBETABINS); + m_calibCoeffHec.resize(s_NBETABINS); + for(int ieta=0 ; ieta<s_NBETABINS ; ieta++) { + m_calibCoeffEmb[ieta]=1.; + m_calibCoeffEmec[ieta]=1.; + m_calibCoeffHec[ieta]=1.; + } + const int nEta = 4; + m_calibCoeffFcalEm.resize(nEta); + m_calibCoeffFcalHad.resize(nEta); + for(int ieta=0 ; ieta<nEta ; ieta++) { + m_calibCoeffFcalEm[ieta]=.03; + m_calibCoeffFcalHad[ieta]=.03; + } + + // + // ........ declare the private data as properties + // + + declareProperty("EventStore",m_storeGateSvc,"StoreGate Service"); + declareProperty("SubDetectors",m_SubDetectors); + declareProperty("RndmSvc", m_atRndmGenSvc); + + declareProperty("EmBarrelHitContainerName",m_EmBarrelHitContainerName); + declareProperty("EmEndCapHitContainerName",m_EmEndCapHitContainerName); + declareProperty("HecHitContainerName",m_HecHitContainerName); + declareProperty("ForWardHitContainerName",m_ForWardHitContainerName); + + declareProperty("EmTTL1ContainerName",m_EmTTL1ContainerName); + declareProperty("HadTTL1ContainerName",m_HadTTL1ContainerName); + + declareProperty("NoiseOnOff",m_NoiseOnOff); + + declareProperty("PileUp",m_PileUp); + declareProperty("UseTriggerTime",m_useTriggerTime); + declareProperty("TriggerTimeToolName",m_triggerTimeToolName); + + declareProperty("EmBarrelCalibrationCoeffs",m_calibCoeffEmb); + declareProperty("EmEndCapCalibrationCoeffs",m_calibCoeffEmec); + declareProperty("HECCalibrationCoeffs",m_calibCoeffHec); + declareProperty("EmFcalCalibrationCoeffs",m_calibCoeffFcalEm); + declareProperty("HadFcalCalibrationCoeffs",m_calibCoeffFcalHad); + + declareProperty("NoEmCalibrationMode",m_noEmCalibMode); + declareProperty("NoHadCalibrationMode",m_noHadCalibMode); + declareProperty("ChronoTest",m_chronoTest); + declareProperty("DebugThreshold",m_debugThresh); + + declareProperty("useMapFromStore",m_useMapfromStore,"Use LArHitEMap already filled from detector store"); + +// +return; +} + + +LArTTL1Maker::~LArTTL1Maker() +{ + /** TTL1 Maker destructor */ + +return; +} + + +StatusCode LArTTL1Maker::initialize() +{ +// +======================================================================+ +// + + +// + Author ........: F. Ledroit + +// + Creation date .: 09/01/2003 + +// + + +// +======================================================================+ +// +// ......... declaration +// + MsgStream msglog(messageService(),name()); + int outputLevel = msgSvc()->outputLevel( name() ); + + m_chronSvc = chronoSvc(); + + msglog << MSG::INFO + << "***********************************************" + << endreq; + msglog << MSG::INFO + << "* Steering options for LArTTL1Maker algorithm *" + << endreq; + msglog << MSG::INFO + << "***********************************************" + << endreq; + // + // ......... print the noise flag + // + if ( m_NoiseOnOff ) + { + msglog << MSG::INFO + << "Electronic noise will be added in each TT for selected sub-detectors." + << endreq; + } + else + { + msglog << MSG::INFO + << "No electronic noise added." + << endreq; + } + +// +// ......... print the pile-up flag +// + if (m_PileUp) + { + msglog << MSG::INFO + << "take events from PileUp service" + << endreq; + } + else + { + msglog << MSG::INFO + << "no pile up" + << endreq; + } + +// +// ......... print the trigger time flag +// + if (m_useTriggerTime) + { + msglog << MSG::INFO + << "use Trigger Time service " << m_triggerTimeToolName + << endreq; + } + else + { + msglog << MSG::INFO + << "no Trigger Time used" + << endreq; + } + +// +// ......... print the calibration flags +// +// currently the calib modes are not used anymore -> turning INFO logs into DEBUG + if (m_noEmCalibMode) + { + msglog << MSG::DEBUG + << "NO calibration mode chosen for EM towers " + << " == technical option. Should not be used for physics !!! " + << endreq; + } + else + { + msglog << MSG::DEBUG + << "standard calibration chosen for EM towers" + << endreq; + } + + if (m_noHadCalibMode) + { + msglog << MSG::DEBUG + << "NO calibration mode chosen for HEC towers " + << " == technical option. Should not be used for physics !!! " + << endreq; + } + else { + msglog << MSG::DEBUG + << "standard calibration mode chosen for HEC towers " + << endreq; + } + +// +// ....... Access Event StoreGate +// + //StatusCode sc = service ( "StoreGateSvc" , m_storeGateSvc ) ; + StatusCode sc = m_storeGateSvc.retrieve(); + if (sc.isFailure()) + { + msglog << MSG::ERROR + << "Unable to access pointer to StoreGate" + << endreq; + return StatusCode::FAILURE; + } + +// +// locate the PileUpMergeSvc and initialize our local ptr +// + if (m_PileUp) { + if (!(service("PileUpMergeSvc", m_mergeSvc)).isSuccess() || 0 == m_mergeSvc) { + msglog << MSG::ERROR + << "Could not find PileUpMergeSvc" + << endreq; + return StatusCode::FAILURE; + } + else + { + msglog << MSG::DEBUG << "PileUpMergeSvc successfully initialized" << endreq; + } + } + +// +// .........retrieve tool computing trigger time if requested +// + if (m_useTriggerTime && m_PileUp) { + msglog << MSG::INFO + << " In case of pileup, the trigger time subtraction is done in PileUpSvc " + << endreq; + msglog << MSG::INFO + << " => LArTTL1Maker will not apply Trigger Time " << endreq; + m_useTriggerTime = false; + } + + if (m_useTriggerTime) { + IToolSvc* p_toolSvc = 0; + if (!(service("ToolSvc", p_toolSvc)).isSuccess() || 0 == p_toolSvc) { + msglog << MSG::ERROR + << "Could not find ToolSvc" + << endreq; + return StatusCode::FAILURE; + } + + IAlgTool* algtool(0); + ListItem theTool(m_triggerTimeToolName.value()); + sc = p_toolSvc->retrieveTool(theTool.type(), theTool.name(),algtool); + if (sc.isFailure()) { + msglog << MSG::ERROR + << "Unable to find tool for " << m_triggerTimeToolName.value() << endreq; + p_triggerTimeTool = 0; + } + else { + p_triggerTimeTool=dynamic_cast<ITriggerTime*>(algtool); + msglog << MSG::DEBUG << "retrieved TriggerTime tool: " + << m_triggerTimeToolName.value() << endreq; + } + } + + // + // ..... need LAr and CaloIdManager to retrieve all needed helpers + // + + const CaloIdManager* caloMgr; + const LArIdManager* larMgr; + // + // ....... Access Detector StoreGate + // + sc = service( "DetectorStore", m_detectorStore ); + if ( sc.isSuccess( ) ) { + sc = m_detectorStore->retrieve(caloMgr); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve CaloIdManager from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved CaloIdManager from DetectorStore" << endreq; + } + } + sc = m_detectorStore->retrieve(larMgr); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve LArIdManager from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved LArIdManager from DetectorStore" << endreq; + } + } + } else { + msglog << MSG::ERROR << "Could not locate DetectorStore" << endreq; + return StatusCode::FAILURE; + } + + // + //..... need of course the LVL1 helper + // + m_lvl1Helper = caloMgr->getLVL1_ID(); + if (!m_lvl1Helper) { + msglog << MSG::ERROR << "Could not access CaloLVL1_ID helper" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully accessed CaloLVL1_ID helper" << endreq; + } + } + + // + //..... also need the LArEM helper (e.g. to deal with the barrel end part) + // + sc = m_detectorStore->retrieve(m_emHelper); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve LArEM helper from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved LArEM helper from DetectorStore" << endreq; + } + } + + // + //..... also need the LArHEC helper to avoid adding up energy from 4th compartment + // + sc = m_detectorStore->retrieve(m_hecHelper); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve LArHEC helper from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved LArHEC helper from DetectorStore" << endreq; + } + } + // + //..... also need the LArFCAL helper to use hash ids to store all gains + // + sc = m_detectorStore->retrieve(m_fcalHelper); + if (sc.isFailure()) { + msglog << MSG::ERROR << "Unable to retrieve LArFCAL helper from DetectorStore" << endreq; + return StatusCode::FAILURE; + } else { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Successfully retrieved LArFCAL helper from DetectorStore" << endreq; + } + } + + // + // ..... need cabling services, to get channels associated to each TT + // + IToolSvc* toolSvc; + StatusCode status = service( "ToolSvc",toolSvc ); + if(status.isSuccess()) { + sc = toolSvc->retrieveTool("LArCablingService",m_cablingSvc); + if(sc.isFailure()) { + msglog << MSG::ERROR << "Could not retrieve LArCablingService"<< endreq; + return(StatusCode::FAILURE); + } + sc = toolSvc->retrieveTool("CaloTriggerTowerService",m_ttSvc); + if(sc.isFailure()) { + msglog << MSG::ERROR << "Could not retrieve CaloTriggerTowerService"<< endreq; + return(StatusCode::FAILURE); + } + } else { + msglog << MSG::ERROR << "Could not get ToolSvc"<< endreq; + return(StatusCode::FAILURE); + } + + // Incident Service: + IIncidentSvc* incSvc; + sc = service("IncidentSvc", incSvc); + if (sc.isFailure()) + { + msglog << MSG::ERROR + << "Unable to retrieve pointer to DetectorStore " + << endreq; + return StatusCode::FAILURE; + } + //start listening to "BeginRun" + incSvc->addListener(this, "BeginRun", m_BeginRunPriority); + + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Initialization completed successfully" + << endreq; + } + + sc = m_atRndmGenSvc.retrieve(); + if(sc.isFailure()) { + msglog << MSG::ERROR << "Could not initialize random number service." << endreq; + return sc; + } + m_rndmEngine = m_atRndmGenSvc->GetEngine(m_rndmEngineName); + if(!m_rndmEngine) { + msglog << MSG::ERROR << "Could not find RndmEngine : " << m_rndmEngineName << endreq; + return StatusCode::FAILURE ; + } + + return StatusCode::SUCCESS; + +} + + + +void LArTTL1Maker::handle(const Incident& /* inc*/ ) +{ + MsgStream msglog( msgSvc(), name() ); + msglog << MSG::DEBUG << "LArTTL1Maker handle()" << endreq; + + StatusCode sc = this->retrieveDatabase(); + if (sc.isFailure()) { + msglog << MSG::ERROR << " Error from retrieveDatabase " << endreq; + } + + + // + // ...... init hit map + // + if ( this->initHitMap() == StatusCode::FAILURE ) { + msglog << MSG::ERROR << " Error from initHitMap() " << endreq; + } + + // + // ...... Read auxiliary data files + // + if ( this->readAuxiliary() == StatusCode::FAILURE ) { + + msglog << MSG::ERROR << " Error from readAuxiliary() " << endreq; + + } + + + return; +} + +StatusCode LArTTL1Maker::retrieveDatabase() +{ +// can not do all that in initialize with the new calibration classes/Iov service + + MsgStream msglog( msgSvc(), name() ); + msglog << MSG::DEBUG << "LArTTL1Maker retrieveDatabase()" << endreq; + StatusCode sc; + + sc = m_detectorStore->retrieve(m_dd_fSampl,m_fSamplKey); + if (sc.isFailure() || !m_dd_fSampl) { + msglog << MSG::ERROR + << "Unable to retrieve ILArfSampl from DetectorStore" + << endreq; + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; +} + + + +StatusCode LArTTL1Maker::execute() +{ +// +======================================================================+ +// + + +// + Author: F. Ledroit + +// + Creation date: 2003/01/13 + +// + Subject: Make the TTL1s and put them into the container + +// + + +// +======================================================================+ + // + // ......... declarations + // + MsgStream msglog(messageService(),name()); + int outputLevel = msgSvc()->outputLevel( name() ); + + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Begining of execution" + << endreq; + } + + // + // ....... fill the LArHitEMap + // + if(m_chronoTest) { + m_chronSvc->chronoStart( "fill LArHitEMap " ); + } + + int totHit=0; + if ( this->fillEMap(totHit) == StatusCode::FAILURE ) return StatusCode::FAILURE; + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "total number of hits with |E|> 1.e-06 found = " << totHit << endreq; + } + + + if(m_chronoTest) { + m_chronSvc->chronoStop ( "fill LArHitEMap " ); + m_chronSvc->chronoPrint( "fill LArHitEMap " ); + } + + if (!m_useMapfromStore && totHit==0) { + msglog << MSG::WARNING << " No LAr hit in the event " << endreq; + } + + // + // .....get the trigger time if requested + // + double trigtime=0; + if (m_useTriggerTime && p_triggerTimeTool) { + trigtime = p_triggerTimeTool->time(); + } + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG << "Trigger time used : " << trigtime << endreq; + } + + // + // ....... create the LAr TTL1 Containers + // + LArTTL1Container *ttL1ContainerEm = new LArTTL1Container(); + if ( ttL1ContainerEm == 0 ){ + msglog << MSG::ERROR + << "Could not allocate a new LArTTL1Container" + << endreq; + return StatusCode::FAILURE; + } + LArTTL1Container *ttL1ContainerHad = new LArTTL1Container(); + if ( ttL1ContainerHad == 0 ){ + msglog << MSG::ERROR + << "Could not allocate a second new LArTTL1Container" + << endreq; + return StatusCode::FAILURE; + } + + // + // ...... register the TTL1 containers into the TES + // + StatusCode sc = m_storeGateSvc->record( ttL1ContainerEm , m_EmTTL1ContainerName) ; + if( sc.isFailure() ) + { + msglog << MSG::ERROR + << "Could not record new LArTTL1Container in TES : " + << m_EmTTL1ContainerName + << endreq; + return StatusCode::FAILURE; + } + sc = m_storeGateSvc->record( ttL1ContainerHad , m_HadTTL1ContainerName) ; + if( sc.isFailure() ) + { + msglog << MSG::ERROR + << "Could not record second new LArTTL1Container in TES : " + << m_HadTTL1ContainerName + << endreq; + return StatusCode::FAILURE; + } + + // + // ... initialise vectors for sums of energy in each TT + // + unsigned int nbTT = (unsigned int)m_lvl1Helper->tower_hash_max() ; + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Max number of LAr Trigger Towers= " << nbTT + << endreq; + } + std::vector<std::vector<float> > sumEnergy ; // inner index = time slot (from 0 to visEvecSize-1) + std::vector<std::vector<float> > sumEnergy2 ; // to allow barrel/endcap separation in 15th TT + FCAL2/3 + sumEnergy.resize(nbTT); + sumEnergy2.resize(nbTT); + std::vector<float> ttSumE; + int ttSumEvecSize=0; + int refTime = 0; + if(!m_PileUp) { + ttSumEvecSize=2*25-1; + refTime = 25-1; + } else { + ttSumEvecSize=s_MAXSAMPLES+s_NBSAMPLES-1; + refTime = s_MAXSAMPLES-1; + } + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Number of time slots considered= " << ttSumEvecSize + << " reference time= " << refTime + << endreq; + } + ttSumE.resize(ttSumEvecSize); + for(unsigned int iTT=0;iTT<nbTT;iTT++){ + sumEnergy[iTT]=ttSumE; + sumEnergy2[iTT]=ttSumE; + } + + m_chronSvc->chronoStart( "LArTTL1Mk hit loop " ); + + int it = 0; + int it_end = m_hitmap->GetNbCells(); + msglog << MSG::DEBUG + << "Size of the hit map= " << it_end + << endreq; + + // + // .... loop on hit lists in the map + // .... and fill sumEnergy vector according to hit membership to TT + // .... (one sum per time slot for each TT) + // + float outOfTimeE=0.; + int nOutOfTime=0; + float inTimeE=0.; + int nInTime=0; + float printEthresh=20.; + int nMissingGain=0; + for( ; it!=it_end;++it) { + LArHitList * hitlist = m_hitmap->GetCell(it); + if (hitlist != 0 ) { + + const std::vector<std::pair<float,float> >* timeE = hitlist->getData(); + if (timeE->size() > 0 ) { + Identifier cellId = hitlist->getIdentifier(); + + int specialCase=0; + bool skipCell=false; + // + // .... skip disconnected cells (normally not in the TT list) + // + if(m_lvl1Helper->is_lar_em(cellId)) { + if(!m_emHelper->is_connected(cellId)) skipCell=true; + } + else if(m_lvl1Helper->is_lar_hec(cellId)) { + if(!m_hecHelper->is_connected(cellId)) skipCell=true; + } + else if(m_lvl1Helper->is_lar_fcal(cellId)) { + if(!m_fcalHelper->is_connected(cellId)) skipCell=true; + } + // + // ... skip cells not summed up in LVL1 (end of barrel PS and 4th compartment of HEC) + // + if(!m_ttSvc->is_in_lvl1(cellId)) skipCell=true; + + if(!skipCell) { + // + // ...determine to which TT this channel belongs + // + Identifier ttId = m_ttSvc->whichTTID(cellId); + + if(m_chronoTest) { + m_chronSvc->chronoStart( "retrieve RG " ); + } + + // + // ....... determine the sampling fraction + //........ and the relative (layer) gains for this cellId + // + float cellSampFraction = m_dd_fSampl->FSAMPL(cellId); + // std::cout << "cellid, SF= " << m_lvl1Helper->show_to_string(cellId) << " " << cellSampFraction << std::endl; + float relGain = 0.; + float sinTheta = 0.; + + int region = m_lvl1Helper->region(ttId); + int eta = m_lvl1Helper->eta(ttId); + int Ieta = decodeInverse(region,eta); + + // case EM cell + if(m_lvl1Helper->is_lar_em(cellId)) { + bool barrelEnd = m_lvl1Helper->is_barrel_end(ttId); + std::vector<float> vecRG ; + if(m_emHelper->is_em_barrel(cellId)) { + // Barrel + sinTheta = m_sinThetaEmb[Ieta-1]; // Ieta starts at 1 + } else { + // End-Cap + if(!barrelEnd) { + sinTheta = m_sinThetaEmec[Ieta-1]; + } else { + // patching the fact that TT offline ID is ambiguous + // decodeInverse(eta,region) returns 15 (ok barrel) instead of 1 (ec value). + sinTheta = m_sinThetaEmec[0]; + specialCase = 1; + msglog << MSG::VERBOSE << " special case " + << m_emHelper->show_to_string(ttId) << endreq ; + } + } + relGain = 1.; // no relative gain for EMB and EMEC + } + + // case HEC cell + else if(m_lvl1Helper->is_lar_hec(cellId)) { + sinTheta = m_sinThetaHec[Ieta-1]; // Ieta starts at 1 + relGain = 1.; + } + + // case FCAL cell + else if(m_lvl1Helper->is_lar_fcal(cellId)) { + IdentifierHash fcalHash = m_fcalHelper->channel_hash(cellId); + relGain = m_cellRelGainFcal [fcalHash] ; + if(relGain < 0.001) { + nMissingGain++; + msglog << MSG::WARNING << " No relative gain value found for FCAL cell " + << m_emHelper->show_to_string(cellId) << " (index " << fcalHash + << "), setting default value 1. " << endreq ; + relGain = 1. ; + } + sinTheta = 1.; // this factor is included in the relative gain + } + if(m_chronoTest) { + m_chronSvc->chronoStop( "retrieve RG " ); + m_chronSvc->chronoPrint( "retrieve RG " ); + } + + IdentifierHash ttHash = m_lvl1Helper->tower_hash(ttId) ; + + // + // .... loop on hits in hit list + // + std::vector<std::pair<float,float> >::const_iterator first = timeE->begin(); + std::vector<std::pair<float,float> >::const_iterator last = timeE->end(); + while (first != last) { + float hitEnergy = (*first).first; + float hitTime = (*first).second - trigtime; + if(hitTime>99000.) { + if(hitEnergy > printEthresh) { + msglog << MSG::WARNING + << " Found pathological hit time, cellId= " << m_emHelper->show_to_string(cellId) + << " time= " << hitTime + << " energy= " << hitEnergy + << endreq; + } + // provisionnaly fix a bug in PileUpEvent. Normally not needed starting from 10.5.0 + hitTime-=100000.; + } + // remove pathological energies (found in some Grid simulated DC2/Rome samples) + if (fabs(hitEnergy)>1e+9) { + msglog << MSG::WARNING << " Pathological energy ignored cellId= " + << m_emHelper->show_to_string(cellId) << " energy= " << hitEnergy << endreq; + hitEnergy = 0.; + hitTime = 0.; + } + // + // ....determine time with respect to ref shape + // + int iShift = 0; + if(!m_PileUp) { + // keep 1ns granularity + //iShift = (int)(hitTime+0.5); + iShift = static_cast<int>(floor(hitTime+0.5)); + } else { + // round to 25ns + //iShift = (int)(hitTime/25.+0.5); + iShift = static_cast<int>(floor(hitTime/25.+0.5)); + } + // + // ....make time positive to allow using it as an index + // + int iTime = iShift + refTime ; + // + // .... keep only hits in the timing window + // + if(iTime>=0 && iTime < ttSumEvecSize ) { + + if(!specialCase) { + // standard case + // ....... make the energy sum + ttSumE = sumEnergy[ttHash]; + ttSumE[iTime] += hitEnergy / cellSampFraction * sinTheta * relGain; + sumEnergy[ttHash] = ttSumE; + } + else { + // ec part of barrel-end or FCAL3 + // ....... make the energy sum + ttSumE = sumEnergy2[ttHash]; + ttSumE[iTime] += hitEnergy / cellSampFraction * sinTheta * relGain; + sumEnergy2[ttHash] = ttSumE; + } + // msglog << MSG::VERBOSE << "applied relative layer gain " << relGain + // << " to a hit of cell " << m_emHelper->show_to_string(cellId) << endreq; + inTimeE+=hitEnergy; + nInTime++; + } // only hits in timing window + else { + outOfTimeE+=hitEnergy; + nOutOfTime++; + if(hitEnergy > printEthresh) { + if(!m_PileUp) { + msglog << MSG::DEBUG + << "Found a hit out of the timing window, hitTime= " << hitTime + << " with more than " << printEthresh << " MeV: hitEnergy= " << hitEnergy << " MeV" + << endreq; + } + else { + if (outputLevel <= MSG::VERBOSE) { + msglog << MSG::VERBOSE + << "Found a hit out of the timing window, hitTime= " << hitTime + << " with more than " << printEthresh << " MeV: hitEnergy= " << hitEnergy << " MeV" + << endreq; + } + } + } + } + ++first; + } // end of loop on hits in the list + } // skip cell condition + } // check timeE->size() > 0 + } // check hitlist > 0 + } // end of loop on hit lists + + msglog << MSG::DEBUG + << "Number of missing relative FCAL gains for this event = " + << nMissingGain + << endreq; + if (outputLevel <= MSG::VERBOSE) { + if (inTimeE == 0 || nInTime == 0) + msglog << MSG::VERBOSE << "No in time energy" << endreq; + else + msglog << MSG::VERBOSE + << "Out of time energy = " << outOfTimeE << " MeV" + << " represents " << 100.* outOfTimeE/inTimeE + << " % of in time energy for " << nOutOfTime << " (" << 100.*nOutOfTime/nInTime << " %) hits" + << endreq; + } + if(outOfTimeE > 0.02*inTimeE) { + msglog << MSG::WARNING + << "Out of time energy = " << outOfTimeE << " MeV" + << " larger than 2% of in time energy = " << inTimeE << " MeV; nb of out of time hits = " + << nInTime << " (" + << (nInTime > 0 ? 100.*nOutOfTime/nInTime : 0) + << " %)" + << endreq; + } + + if(m_chronoTest) { + m_chronSvc->chronoStop( "LArTTL1Mk hit loop " ); + m_chronSvc->chronoPrint( "LArTTL1Mk hit loop " ); + m_chronSvc->chronoStart( "LArTTL1Mk TT loop " ); + } + + std::vector<std::string> emHadString(2); + emHadString[0]="ElectroMagnetic"; + emHadString[1]="Hadronic"; + + std::vector<Identifier>::const_iterator itTt = m_lvl1Helper->tower_begin() ; + std::vector<Identifier>::const_iterator itEnd = m_lvl1Helper->tower_end() ; + + // + // ....... loop on Trigger Towers + // ....... and build signals using pulse shape and noise for each TT + // + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << " Starting loop on Trigger Towers " + << endreq; + } + for(; itTt!=itEnd;++itTt){ + + Identifier towerId = (*itTt) ; + + // + // ........ skip Tile cal + // + if(!m_lvl1Helper->is_tile(towerId)){ + + int emHad = m_lvl1Helper->sampling(towerId); + // convert offline id to hash id + IdentifierHash ttHash = m_lvl1Helper->tower_hash(towerId) ; + int region = m_lvl1Helper->region(towerId); + int eta = m_lvl1Helper->eta(towerId); + int Ieta = decodeInverse(region,eta); + + // + // .... compute the signal for current trigger tower + // + if(m_chronoTest) { + m_chronSvc->chronoStart( "compute signal " ); + } + std::vector<float> analogSum (s_NBSAMPLES) ; + std::vector<float> analogSum2(s_NBSAMPLES) ; + std::vector<float> sumTTE = sumEnergy[ttHash]; + std::vector<float> sumTTE2 = sumEnergy2[ttHash]; + int nSpecialCase=0; + + + bool hasE = false; + for ( unsigned int i=0; i<sumTTE.size();++i){ + if (fabs(sumTTE[i]) > 0.){ + hasE=true; + break; + } + } + + bool hasE2 = false; + for ( unsigned int i=0; i<sumTTE2.size();++i){ + if (fabs(sumTTE2[i]) > 0.){ + hasE2=true; + break; + } + } + + + //if ( fabs(sumTTE[refTime]) > 0.|| fabs(sumTTE2[refTime]) > 0. ) { + if ( hasE || hasE2 ) { + + //if(fabs(sumTTE[refTime]) > 0.) { + if( hasE ) { + analogSum = computeSignal(towerId,Ieta,0,sumTTE,refTime) ; + } + + // treat special case of ec part of the "barrel end" and special case of FCAL2/FCAL3 + // fix me !! this way, we double the saturation energy in this tower...(because it is cut into 2 pieces) + //if(fabs(sumTTE2[refTime]) > 0.) { + if( hasE2) { + nSpecialCase+=1; + if(nSpecialCase>1) { + msglog << MSG::WARNING + << " more than 1 special case, current Trigger Tower is " << emHadString[emHad] << ": " + << m_emHelper->show_to_string(towerId) << " " + << endreq; + } + analogSum2 = computeSignal(towerId,Ieta,1,sumTTE2,refTime) ; + for( int isamp=0 ; isamp < s_NBSAMPLES ; isamp++) { + analogSum[isamp] += analogSum2[isamp] ; + } + } + + if ( sumTTE[refTime] > m_debugThresh) { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << " current Trigger Tower is " << emHadString[emHad] << ": " + << m_emHelper->show_to_string(towerId) << " " + << endreq; + msglog << MSG::DEBUG + << " transverse E (i.e. sum E / sampling fraction * sin_theta * rel gain)= (at ref. time, before calib)" + << sumTTE[refTime] << " + " << sumTTE2[refTime] << " (special cases) " + << endreq; + } + } + else if ( sumTTE[refTime] > 0.) { + if (outputLevel <= MSG::VERBOSE) { + msglog << MSG::VERBOSE + << " current Trigger Tower is " << emHadString[emHad] << ": " + << m_emHelper->show_to_string(towerId) << " " + << endreq; + msglog << MSG::VERBOSE + << " [very low] transverse E (i.e. sum E / sampling fraction * sin_theta * rel gain)= (at ref. time, before calib)" + << sumTTE[refTime] << " + " << sumTTE2[refTime] << " (special cases) " + << endreq; + } + } + + } + if(m_chronoTest) { + m_chronSvc->chronoStop ( "compute signal " ); + m_chronSvc->chronoPrint( "compute signal " ); + } + + // + // ........ add the noise + // + if(m_chronoTest) { + m_chronSvc->chronoStart( "adding noise " ); + } + std::vector<float> fullSignal(s_NBSAMPLES) ; + if(m_NoiseOnOff) { + fullSignal = computeNoise(towerId,Ieta,analogSum); + } else { + fullSignal = analogSum; + } + + if(m_chronoTest) { + m_chronSvc->chronoStop ( "adding noise " ); + m_chronSvc->chronoPrint ( "adding noise " ); + } + + if(sumTTE[refTime] > m_debugThresh) { + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << " uncalibrated amplitudes around peak (+-3 time slots): " + << sumTTE[refTime-3] << ", " + << sumTTE[refTime-2] << ", " + << sumTTE[refTime-1] << ", " + << sumTTE[refTime] << ", " + << sumTTE[refTime+1] << ", " + << sumTTE[refTime+2] << ", " + << sumTTE[refTime+3] + << endreq; + msglog << MSG::DEBUG + << " calibrated signal is " + << analogSum[0] << ", " + << analogSum[1] << ", " + << analogSum[2] << ", " + << analogSum[3] << ", " + << analogSum[4] << ", " + << analogSum[5] << ", " + << analogSum[6] + << endreq; + msglog << MSG::DEBUG + << " shape of calibrated signal is " + << analogSum[0]/analogSum[3] << ", " + << analogSum[1]/analogSum[3] << ", " + << analogSum[2]/analogSum[3] << ", " + << analogSum[3]/analogSum[3] << ", " + << analogSum[4]/analogSum[3] << ", " + << analogSum[5]/analogSum[3] << ", " + << analogSum[6]/analogSum[3] + << endreq; + msglog << MSG::DEBUG + << " after adding noise, full signal is " + << fullSignal[0] << ", " + << fullSignal[1] << ", " + << fullSignal[2] << ", " + << fullSignal[3] << ", " + << fullSignal[4] << ", " + << fullSignal[5] << ", " + << fullSignal[6] + << endreq; + } + if (outputLevel <= MSG::VERBOSE) { + for (unsigned int iTime=0; iTime<sumTTE.size(); iTime++) { + msglog << MSG::VERBOSE + << " iTime [range=0-28 or 0-299] = " << iTime + << " hit energy = " << sumTTE[iTime] + << endreq; + } + } + } else if(sumTTE[refTime] > 0.) { + if (outputLevel <= MSG::VERBOSE) { + msglog << MSG::VERBOSE + << " uncalibrated amplitudes around peak (+-3 time slots): " + << sumTTE[refTime-3] << ", " + << sumTTE[refTime-2] << ", " + << sumTTE[refTime-1] << ", " + << sumTTE[refTime] << ", " + << sumTTE[refTime+1] << ", " + << sumTTE[refTime+2] << ", " + << sumTTE[refTime+3] + << endreq; + msglog << MSG::VERBOSE + << " calibrated signal is " + << analogSum[0] << ", " + << analogSum[1] << ", " + << analogSum[2] << ", " + << analogSum[3] << ", " + << analogSum[4] << ", " + << analogSum[5] << ", " + << analogSum[6] + << endreq; + msglog << MSG::VERBOSE + << " shape of calibrated signal is " + << analogSum[0]/analogSum[3] << ", " + << analogSum[1]/analogSum[3] << ", " + << analogSum[2]/analogSum[3] << ", " + << analogSum[3]/analogSum[3] << ", " + << analogSum[4]/analogSum[3] << ", " + << analogSum[5]/analogSum[3] << ", " + << analogSum[6]/analogSum[3] + << endreq; + msglog << MSG::VERBOSE + << " after adding noise, full signal is " + << fullSignal[0] << ", " + << fullSignal[1] << ", " + << fullSignal[2] << ", " + << fullSignal[3] << ", " + << fullSignal[4] << ", " + << fullSignal[5] << ", " + << fullSignal[6] + << endreq; + } + } + + if(fabs(fullSignal[s_PEAKPOS]) > 0.) { + // + // ...... create the LArTTL1 and push it into the appropriate TTL1 container + // + LArTTL1 * ttL1 ; + HWIdentifier ttChannel; + ttL1 = new LArTTL1(ttChannel,towerId,fullSignal); + if (emHad) { ttL1ContainerHad->push_back(ttL1); } + else { ttL1ContainerEm->push_back(ttL1); } + } + + } // end excluding Tile cal + } // end of for loop on TT + if(m_chronoTest) { + m_chronSvc->chronoStop ( "LArTTL1Mk TT loop " ); + m_chronSvc->chronoPrint ( "LArTTL1Mk TT loop " ); + } + + msglog << MSG::DEBUG << "number of created TTL1s (Em, Had) = " + << ttL1ContainerEm->size() + << " , " + << ttL1ContainerHad->size() + << endreq; + + // + // ...... lock the TTL1 containers + // + sc = m_storeGateSvc->setConst( ttL1ContainerEm ) ; + if( sc.isFailure() ) + { + msglog << MSG::ERROR + << "Could not lock em LArTTL1Container : " + << m_EmTTL1ContainerName + << endreq; + return StatusCode::FAILURE; + } + sc = m_storeGateSvc->setConst( ttL1ContainerHad ) ; + if( sc.isFailure() ) + { + msglog << MSG::ERROR + << "Could not lock had LArTTL1Container : " + << m_HadTTL1ContainerName + << endreq; + return StatusCode::FAILURE; + } + + + return StatusCode::SUCCESS; + +} + + +StatusCode LArTTL1Maker::finalize() +{ +// +======================================================================+ +// + + +// + Author: F. Ledroit + +// + Creation date: 2003/01/13 + +// + + +// +======================================================================+ +// +// ......... declaration +// + MsgStream msglog(messageService(),name()); + + msglog << MSG::INFO << " LArTTL1Maker finalize completed successfully" + << endreq; + + m_chronSvc->chronoPrint( "LArTTL1Mk hit loop " ); + m_chronSvc->chronoPrint( "LArTTL1Mk TT loop " ); + + return StatusCode::SUCCESS; + +} + + +std::vector<float> LArTTL1Maker::computeSignal(const Identifier towerId, const int Ieta, const int specialCase, + std::vector<float> ttSumEnergy, const int refTime) const + +{ +// +======================================================================+ +// + + +// + Author: F. Ledroit + +// + Creation date: 2003/01/13 + +// + + +// +======================================================================+ +// + + MsgStream msglog(messageService(),name()); + int outputLevel = msgSvc()->outputLevel( name() ); + + std::vector<float> bareSignal(s_NBSAMPLES) ; + std::vector<float> pulseShape(s_MAXSAMPLES) ; + std::vector<float> pulseShapeDer(s_MAXSAMPLES) ; + + bool emb = m_lvl1Helper->is_emb(towerId); + bool barrelEnd = m_lvl1Helper->is_barrel_end(towerId); + bool emec= m_lvl1Helper->is_emec(towerId); + + unsigned int visEvecSize=ttSumEnergy.size(); + if (outputLevel <= MSG::VERBOSE) { + msglog << MSG::VERBOSE << "computeSignal: special case = " + << specialCase << endreq ; + } + + // + // ..... case EM + // + if (emb || barrelEnd || emec) { + // + // ..... retrieve calib coeffs for the tower + // + float calibCoeff=0.; + if (emb) { + calibCoeff = m_calibCoeffEmb[Ieta-1]; + } + else if (emec) { + calibCoeff = m_calibCoeffEmec[Ieta-1]; + } + else { // barrelEnd + if(!specialCase) { + calibCoeff = m_calibCoeffEmb[14]; + } + else { + calibCoeff = m_calibCoeffEmec[0]; + } + } + if(calibCoeff < 0.001) { + msglog << MSG::WARNING << " No calibration coefficient value found for tower " + << m_emHelper->show_to_string(towerId) + << " setting default value 6. " + << endreq ; + calibCoeff = 6.; + } + + // + // ... loop on time samples + // + for (unsigned int iTime=0; iTime<visEvecSize; iTime++) { + if(fabs(ttSumEnergy[iTime]) > 0.) { + if(!m_noEmCalibMode) { + // apply calibration coefficient + ttSumEnergy[iTime] *= calibCoeff; + if (outputLevel <= MSG::VERBOSE) { + msglog << MSG::VERBOSE << " ComputeSignal: applied EM calibration coefficient (iTime) " + << calibCoeff << " (" << iTime << ") " << endreq; + } + } + + // with respect to saturation + // to compute the amplitude + float theEnergy = ttSumEnergy[iTime] ; + if( ttSumEnergy[iTime] > m_refEnergyEm [0]) { + theEnergy = m_refEnergyEm [0] ; + } + // to determine the shape + if( ttSumEnergy[iTime] > m_refEnergyEm [s_NBENERGIES-1]) { + // don't know the shape after highest energy point -> stick to 'last' shape + ttSumEnergy[iTime] = m_refEnergyEm [s_NBENERGIES-1] ; + } + + // + // ... determine regime: linear(iene=0) or saturation + // + for (int iene = 0 ; iene < s_NBENERGIES ; iene++) { + if ( ttSumEnergy[iTime] <= m_refEnergyEm[iene] ) { + pulseShape = m_pulseShapeEm[iene] ; + pulseShapeDer = m_pulseShapeDerEm[iene] ; + + // make samples + for( int iSamp=0 ; iSamp < s_NBSAMPLES ; iSamp++) { + if(!m_PileUp) { + // go back to signed value of time + int hitTime = iTime - refTime; + // go to 25 ns sampling + //int time = (int)(hitTime/25.+0.5) ; + int time = static_cast<int>(floor(hitTime/25.+0.5)); + // ....determine fine shift within 25ns + float dTime = (float)(hitTime - 25*time); + int j= iSamp - time ; + if(j >=0 && j < s_MAXSAMPLES) { + bareSignal[iSamp] += (pulseShape[j] - pulseShapeDer[j]*dTime) * theEnergy ; + } + } + else { + int j= iSamp - iTime + refTime ; + if(j >=0 && j < s_MAXSAMPLES) { + bareSignal[iSamp] += pulseShape[j] * theEnergy ; + } + } + } // end loop on samples + + break ; + } + } // end of loop on reference energies + } // end condition visE>0 + } // end loop on times + } // end case EM + // + // ..... case HEC + // + else if (m_lvl1Helper->is_hec(towerId)) { + pulseShape = m_pulseShapeHec ; + pulseShapeDer = m_pulseShapeDerHec ; + // + // ... loop on time samples + // + for (unsigned int iTime=0; iTime<visEvecSize; iTime++) { + float theEnergy = ttSumEnergy[iTime]; + if(!m_noHadCalibMode) { + // apply calibration coefficient + theEnergy *= m_calibCoeffHec[Ieta-1]; + } + + float satEt = m_satEnergyHec[Ieta-1]; + for( int iSamp=0 ; iSamp < s_NBSAMPLES ; iSamp++) { + if(!m_PileUp) { + // go back to signed value of time + int hitTime = iTime - refTime; + // go to 25 ns sampling + //int time = (int)(hitTime/25.+0.5) ; + int time = static_cast<int>(floor(hitTime/25.+0.5)); + // ....determine fine shift within 25ns + float dTime = (float)(hitTime - 25*time); + int j= iSamp - time ; + if(j >=0 && j < s_MAXSAMPLES) { + bareSignal[iSamp] += (pulseShape[j] - pulseShapeDer[j]*dTime) * theEnergy ; + } + } + else { + int j= iSamp - iTime + refTime ; + if(j >=0 && j < s_MAXSAMPLES) { + bareSignal[iSamp] += pulseShape[j] * theEnergy ; + } + } + // saturation + if ( bareSignal[iSamp] > satEt) { + bareSignal[iSamp] = satEt; + } + } + } + } + // + // ..... case FCAL + // + else if (m_lvl1Helper->is_fcal(towerId)) { + + int emOrHad = m_lvl1Helper->sampling(towerId) ; + int module = emOrHad + 1 ; + std::vector<float> calibCoeff2= m_calibCoeffFcal[module-1]; + if( emOrHad && (Ieta==3 || Ieta==4) ) { + module+=1; // FCAL3 + } + pulseShape = m_pulseShapeFcal[module-1] ; + pulseShapeDer = m_pulseShapeDerFcal[module-1] ; + + // + // ... loop on time samples (only one if no pile-up) + // + for (unsigned int iTime=0; iTime<visEvecSize; iTime++) { + float theEnergy = ttSumEnergy[iTime]; + if((!m_noEmCalibMode&&module==1)||(!m_noHadCalibMode&&module>1)) { + // apply calibration coefficient + theEnergy *= calibCoeff2[Ieta-1]; + } + + // make samples + for( int iSamp=0 ; iSamp < s_NBSAMPLES ; iSamp++) { + if(!m_PileUp) { + // go back to signed value of time + int hitTime = iTime - refTime; + // go to 25 ns sampling + //int time = (int)(hitTime/25.+0.5) ; + int time = static_cast<int>(floor(hitTime/25.+0.5)); + // ....determine fine shift within 25ns + float dTime = (float)(hitTime - 25*time); + int j= iSamp - time ; + if(j >=0 && j < s_MAXSAMPLES) { + bareSignal[iSamp] += (pulseShape[j] - pulseShapeDer[j]*dTime) * theEnergy ; + } + } + else { + int j= iSamp - iTime + refTime ; + if(j >=0 && j < s_MAXSAMPLES) { + bareSignal[iSamp] += pulseShape[j] * theEnergy ; + } + } + } + } + } + else { + msglog << MSG::WARNING << " LArTTL1Maker: computeSignal unknown towerId " + << m_lvl1Helper->show_to_string(towerId) + << endreq; + return bareSignal ; + } + + return bareSignal ; + +} + + +std::vector<float> LArTTL1Maker::computeNoise(const Identifier towerId, const int Ieta, + std::vector<float>& inputV) +{ + +// +======================================================================+ +// + + +// + Author: F. Ledroit + +// + Creation date: 2003/01/13 + +// + + +// +======================================================================+ + + MsgStream msglog(messageService(),name()); + + std::vector<float> outputV (s_NBSAMPLES) ; + + // + // .... first retrieve auto-correlation matrix and noise rms + // + + float sigmaNoise = 0; + std::vector<float>* autoC = 0; + std::vector<float> noiseRms(4) ; + + bool emb = m_lvl1Helper->is_emb(towerId); + bool barrelEnd = m_lvl1Helper->is_barrel_end(towerId); + bool emec= m_lvl1Helper->is_emec(towerId); + // + // ..... case EM + // + if (emb || barrelEnd || emec) { + // + // ..... retrieve noise info for the tower + // + autoC = &(m_autoCorrEm); + if (emb) { + sigmaNoise = m_noiseRmsEmb[Ieta-1]; + } + else if (emec) { + sigmaNoise = m_noiseRmsEmec[Ieta-1]; + } + else { // barrelEnd + sigmaNoise = sqrt(m_noiseRmsEmb[14]*m_noiseRmsEmb[14] + m_noiseRmsEmec[0]*m_noiseRmsEmec[0]); + msglog << MSG::VERBOSE << " LArTTL1Maker: barrelEnd noise rms" + << sigmaNoise + << endreq; + } + } // end case EM + // + // ..... case HEC + // + else if (m_lvl1Helper->is_hec(towerId)) { + + autoC = &(m_autoCorrHec[Ieta-1]); + sigmaNoise = m_noiseRmsHec[Ieta-1]; + } + // + // ..... case FCAL + // + else if (m_lvl1Helper->is_fcal(towerId)) { + + int emOrHad = m_lvl1Helper->sampling(towerId) ; + int module = emOrHad + 1 ; + if( emOrHad && (Ieta==3 || Ieta==4) ) { + module+=1; // FCAL3 + } + autoC = &m_autoCorrFcal; + // sigmaNoise = m_noiseRmsFcal[module-1]; + noiseRms = m_noiseRmsFcal[module-1]; + sigmaNoise = noiseRms[Ieta-1]; + msglog << MSG::VERBOSE << " LArTTL1Maker: noise FCAL= " << sigmaNoise + << " module= " << module << "Ieta= " << Ieta << endreq; + } + else { + msglog << MSG::WARNING << " LArTTL1Maker: computeNoise unknown towerId " + << m_lvl1Helper->show_to_string(towerId) + << endreq; + return inputV ; + } + + if(fabs((*autoC)[0]) < 0.001) { + msglog << MSG::WARNING << " No autocorrelation matrix found for tower " + << m_emHelper->show_to_string(towerId) << " " + << "setting default values 1.00 0.10 -0.30 -0.20 -0.05 -0.01 -0.01 " + << endreq ; + (*autoC)[0] = 1.00; + (*autoC)[1] = 0.10; + (*autoC)[2] =-0.30; + (*autoC)[3] =-0.20; + (*autoC)[4] =-0.05; + (*autoC)[5] =-0.01; + (*autoC)[6] =-0.01; + } + if(sigmaNoise < 0.001) { + msglog << MSG::WARNING << " No noise rms value found for tower " + << m_emHelper->show_to_string(towerId) << " " + << "setting default value 300 MeV " + << endreq ; + sigmaNoise = 300.; + } + + // + // .... now compute the noise 'signal' + // + + float c11,c21,c31,c41,c51,c22,c32,c33,c42,c43,c44,c52,c53,c54,c55; + float c61,c62,c63,c64,c65,c66,c71,c72,c73,c74,c75,c76,c77; + + c11 = sigmaNoise*(*autoC)[0]; + c21 = sigmaNoise*(*autoC)[1]; + c31 = sigmaNoise*(*autoC)[2]; + c41 = sigmaNoise*(*autoC)[3]; + c51 = sigmaNoise*(*autoC)[4]; + c61 = sigmaNoise*(*autoC)[5]; + c71 = sigmaNoise*(*autoC)[6]; +// + c22 = sqrt(c11*c11-c21*c21); + c32 = (c21*c11-c21*c31)/c22; + c33 = sqrt(c11*c11-c31*c31-c32*c32); + c42 = (c31*c11-c21*c41)/c22; + c43 = (c21*c11-c31*c41-c32*c42)/c33; + c44 = sqrt(c11*c11-c41*c41-c42*c42-c43*c43); + c52 = (c41*c11-c21*c51)/c22; + c53 = (c31*c11-c31*c51-c32*c52)/c33; + c54 = (c21*c11-c41*c51-c42*c52-c43*c53)/c44; + c55 = sqrt(c11*c11-c51*c51-c52*c52-c53*c53-c54*c54); + c62 = (c51*c11-c21*c61)/c22; + c63 = (c41*c11-c31*c61-c32*c62)/c33; + c64 = (c31*c11-c41*c61-c42*c62-c43*c63)/c44; + c65 = (c21*c11-c51*c61-c52*c62-c53*c63-c54*c64)/c55; + c66 = sqrt(c11*c11-c61*c61-c62*c62-c63*c63-c64*c64-c65*c65); + c72 = (c61*c11-c21*c71)/c22; + c73 = (c51*c11-c31*c71-c32*c72)/c33; + c74 = (c41*c11-c41*c71-c42*c72-c43*c73)/c44; + c75 = (c31*c11-c51*c71-c52*c72-c53*c73-c54*c74)/c55; + c76 = (c21*c11-c61*c71-c62*c72-c63*c73-c64*c74-c65*c75)/c66; + c77 = sqrt(c11*c11-c71*c71-c72*c72-c73*c73-c74*c74-c75*c75-c76*c76); + + double rndm[s_NBSAMPLES]; + RandGaussZiggurat::shootArray(m_rndmEngine,(int)s_NBSAMPLES,rndm,0.,1.); + outputV[0] = inputV[0] + c11*rndm[0]; + outputV[1] = inputV[1] + c21*rndm[0] + c22*rndm[1]; + outputV[2] = inputV[2] + c31*rndm[0] + c32*rndm[1] + c33*rndm[2]; + outputV[3] = inputV[3] + c41*rndm[0] + c42*rndm[1] + c43*rndm[2] + c44*rndm[3]; + outputV[4] = inputV[4] + c51*rndm[0] + c52*rndm[1] + c53*rndm[2] + c54*rndm[3] + c55*rndm[4]; + outputV[5] = inputV[5] + c61*rndm[0] + c62*rndm[1] + c63*rndm[2] + c64*rndm[3] + c65*rndm[4] + c66*rndm[5]; + outputV[6] = inputV[6] + c71*rndm[0] + c72*rndm[1] + c73*rndm[2] + c74*rndm[3] + c75*rndm[4] + c76*rndm[5] + c77*rndm[6]; + + return outputV ; + +} + + +StatusCode LArTTL1Maker::fillEMap(int& totHit) + +{ +// +======================================================================+ +// + + +// + Author: F. Ledroit + +// + Creation date: 2003/01/14 + +// + + +// +======================================================================+ + + MsgStream msglog(messageService(),name()); + int outputLevel = msgSvc()->outputLevel( name() ); + +// +// ........ reset the map Energies +// + if (m_useMapfromStore) { + StatusCode sc = m_detectorStore->retrieve(m_hitmap,"LArHitEMap"); + if (sc.isFailure()) { + msglog << " cannot retrieve LArHitEMap from detector Store " << endreq; + return sc; + } + } + + else { + + m_hitmap->EnergyReset(); + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "Execute: energy reset done" + << endreq; + } + + + Identifier cellId; + float hitEnergy; + float hitTime; + int skipHit=0; + + // + // ............ loop over the wanted hit containers (one per sub-detector) + // + + for (unsigned int iHitContainer=0;iHitContainer<m_HitContainer.size();iHitContainer++) + { + // + // ..... Get the pointer to the Hit Container from StoreGate + // + if (outputLevel <= MSG::VERBOSE) { + msglog << MSG::VERBOSE << " asking for: " + << m_HitContainer[iHitContainer] + << endreq; + } + + if (!m_PileUp) { + const LArHitContainer* hit_container ; + StatusCode sc = m_storeGateSvc->retrieve( hit_container, + m_HitContainer[iHitContainer] ) ; + + if ( sc.isFailure() || !hit_container) { + msglog << MSG::ERROR << "Could not retrieve a LArHitContainer " + << m_HitContainer[iHitContainer] + << endreq; + return StatusCode::FAILURE; + } + + if (outputLevel <= MSG::VERBOSE) { + msglog << MSG::VERBOSE + << "number of hits: " << hit_container->size() + << endreq; + } + + // + // ....... loop over hits and get informations + // + + LArHitContainer::const_iterator hititer; + for(hititer=hit_container->begin();hititer != hit_container->end();hititer++) { + cellId = (*hititer)->cellID(); + hitEnergy = (float)(*hititer)->energy(); + hitTime = (float)(*hititer)->time(); + + // + // ....... fill the Map ; don't keep hits with ridiculously small energy + // + if (fabs(hitEnergy) > 1.e-06) { + if ( !m_hitmap->AddEnergy(cellId,hitEnergy,hitTime) ) { + msglog << MSG::FATAL << " Cell " << cellId.getString() + << " could not add the energy= " << hitEnergy << " (MeV)" + << endreq; + return(StatusCode::FAILURE); + } + ++totHit; + } + else { + ++skipHit; + } + + } // .... end of loop over hits + } // ... end of NO pile-up condition + else { + + typedef PileUpMergeSvc::TimedList<LArHitContainer>::type TimedHitContList; + TimedHitContList hitContList; + // + // ...retrieve list of pairs (time,container) from PileUp service + // + + if (!(m_mergeSvc->retrieveSubEvtsData(m_HitContainer[iHitContainer] + ,hitContList).isSuccess()) && hitContList.size()==0) { + msglog << MSG::ERROR + << "Could not fill TimedHitContList" + << endreq; + return StatusCode::FAILURE; + } + + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "number of containers in the list: " << hitContList.size() + << endreq; + } + + // + // ...loop over this list + // + TimedHitContList::iterator iFirstCont(hitContList.begin()); + TimedHitContList::iterator iEndCont(hitContList.end()); + double SubEvtTimOffset; + + while (iFirstCont != iEndCont) { + // get time for this subevent + // new coding of time information (January 05) + const PileUpTimeEventIndex* time_evt = &(iFirstCont->first); + SubEvtTimOffset = time_evt->time(); + // SubEvtTimOffset = iFirstCont->first; + + // get LArHitContainer for this subevent + const LArHitContainer& firstCont = *(iFirstCont->second); + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "number of hits in container: " << firstCont.size() + << ", first five are:" + << endreq; + } + + int numHit=0; + // Loop over cells in this LArHitContainer + LArHitContainer::const_iterator f_cell=firstCont.begin(); + LArHitContainer::const_iterator l_cell=firstCont.end(); + while (f_cell != l_cell) { + hitEnergy = (float) (*f_cell)->energy(); + hitTime = (float) SubEvtTimOffset; + cellId = (*f_cell)->cellID(); + ++f_cell; + + if(numHit<5) { + msglog << MSG::DEBUG + << "cellId " << m_lvl1Helper->show_to_string(cellId) + << ", energy= " << hitEnergy + << ", time= " << hitTime + << endreq; + } + ++numHit; + ++totHit; + + // + // ....... fill the Map + // + if ( !m_hitmap->AddEnergy(cellId,hitEnergy,hitTime) ) + { + msglog << MSG::FATAL << " Cell " << cellId.getString() + << " could not add the energy= " << hitEnergy << " (?eV)" + << endreq; + return(StatusCode::FAILURE); + } + + } // ... end of while loop on hits + + ++iFirstCont; + } // ... end of while loop on containers + + } // ... end of pile-up condition + } // .... end of loop over containers + + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::DEBUG + << "skipped " << skipHit << " hits with abs(energy) less than 1.e-06 " << endreq; + } + + } + + return StatusCode::SUCCESS; +} + + +StatusCode LArTTL1Maker::readAuxiliary() +{ + // + // ...... Read auxiliary data file for EM (barrel and EC) + // + + MsgStream msglog( msgSvc(), name() ); + msglog << MSG::DEBUG << "executing readAuxiliary()" << endreq; + + float refEnergy; + std::vector<float> pulseShape(s_MAXSAMPLES) ; + std::vector<float> pulseShapeDer(s_MAXSAMPLES) ; + std::string barrel_endcap ; + int Ieta ; + float sinTheta=0. ; + std::vector<float> layerRelGain(s_NBDEPTHS) ; + std::vector<float> noiseRms(3) ; + std::vector<float> noiseRms4(4) ; + std::vector<float> autoCorr(s_NBSAMPLES) ; + + + std::string pulsedataname=PathResolver::find_file ("LArEmLvl1.data", "DATAPATH"); + if (pulsedataname == "") { + msglog << MSG::ERROR << "Could not locate LArEmLvl1.data file" << endreq; + return StatusCode::FAILURE; + } + const char * pulsedatafile= pulsedataname.c_str() ; + std::ifstream infile (pulsedatafile) ; + + if(!infile) { + msglog << MSG::ERROR + << " cannot open EM file " + << endreq; + return StatusCode::FAILURE; + } else { + msglog << MSG::DEBUG + << " EM file opened " + << endreq; + } + + // first read the pulse shape for all energies + // valid for both barrel and endcap (from Xavier de la Broise) + for(int iene=0 ; iene<s_NBENERGIES ; iene++) { + infile >> refEnergy + >> pulseShape[0] + >> pulseShape[1] + >> pulseShape[2] + >> pulseShape[3] + >> pulseShape[4] + >> pulseShape[5] + >> pulseShape[6] + >> pulseShape[7] + >> pulseShape[8] + >> pulseShape[9] + >> pulseShape[10] + >> pulseShape[11] + >> pulseShape[12] + >> pulseShape[13] + >> pulseShape[14] + >> pulseShape[15] + >> pulseShape[16] + >> pulseShape[17] + >> pulseShape[18] + >> pulseShape[19] + >> pulseShape[20] + >> pulseShape[21] + >> pulseShape[22] + >> pulseShape[23]; + m_refEnergyEm.push_back(refEnergy); + m_pulseShapeEm.push_back(pulseShape); + infile >> pulseShapeDer[0] + >> pulseShapeDer[1] + >> pulseShapeDer[2] + >> pulseShapeDer[3] + >> pulseShapeDer[4] + >> pulseShapeDer[5] + >> pulseShapeDer[6] + >> pulseShapeDer[7] + >> pulseShapeDer[8] + >> pulseShapeDer[9] + >> pulseShapeDer[10] + >> pulseShapeDer[11] + >> pulseShapeDer[12] + >> pulseShapeDer[13] + >> pulseShapeDer[14] + >> pulseShapeDer[15] + >> pulseShapeDer[16] + >> pulseShapeDer[17] + >> pulseShapeDer[18] + >> pulseShapeDer[19] + >> pulseShapeDer[20] + >> pulseShapeDer[21] + >> pulseShapeDer[22] + >> pulseShapeDer[23]; + m_pulseShapeDerEm.push_back(pulseShapeDer); + } + + // then read autocorrelation function valid for barrel and endcap + // ("default" auto-correlation function) + m_autoCorrEm.resize(s_NBSAMPLES); + infile >> m_autoCorrEm[0] + >> m_autoCorrEm[1] + >> m_autoCorrEm[2] + >> m_autoCorrEm[3] + >> m_autoCorrEm[4] + >> m_autoCorrEm[5] + >> m_autoCorrEm[6] ; + + // then read the separator for barrel + infile >> barrel_endcap; + + for(int ieta=0 ; ieta<s_NBETABINS ; ieta++) { + // then read the calibration factor (from MC) + // and the noise rms (from Bill Cleland) for each eta bin, in barrel + infile >> Ieta + >> sinTheta; + m_sinThetaEmb.push_back(sinTheta); + infile >> noiseRms[0] + >> noiseRms[1] + >> noiseRms[2] ; + m_noiseRmsEmb.push_back(noiseRms[2]); + // if later we want to disentangle between pre-sum and summing electronic noise... + // m_noiseRmsEmb.push_back(noiseRms); + } + + // now read the separator for endcap + infile >> barrel_endcap; + + for(int ieta=0 ; ieta<s_NBETABINS ; ieta++) { + // then read the calibration coefficient (from MC) + // and the noise rms (from Bill Cleland) for each eta bin, in barrel + infile >> Ieta + >> sinTheta; + m_sinThetaEmec.push_back(sinTheta); + infile >> noiseRms[0] + >> noiseRms[1] + >> noiseRms[2] ; + m_noiseRmsEmec.push_back(noiseRms[2]); + // if later we want to disentangle between pre-sum and summing electronic noise... + // m_noiseRmsEmec.push_back(noiseRms); + + } + + infile.close() ; + msglog << MSG::DEBUG + << " EM file closed " + << endreq; + + msglog << MSG::INFO + << " 1 pulse shape per energy for EMB+EMEC : " + << endreq; + for(int iene=0;iene<s_NBENERGIES;iene++) { + msglog << MSG::INFO + << m_refEnergyEm[iene] << " MeV: " + << (m_pulseShapeEm[iene])[0] << ", " + << (m_pulseShapeEm[iene])[1] << ", " + << (m_pulseShapeEm[iene])[2] << ", " + << (m_pulseShapeEm[iene])[3] << ", " + << (m_pulseShapeEm[iene])[4] << ", " + << (m_pulseShapeEm[iene])[5] << ", " + << (m_pulseShapeEm[iene])[6] << ", " + << (m_pulseShapeEm[iene])[7] << ", " + << (m_pulseShapeEm[iene])[8] << ", " + << (m_pulseShapeEm[iene])[9] << ", " + << (m_pulseShapeEm[iene])[10] << ", " + << (m_pulseShapeEm[iene])[11] << ", " + << (m_pulseShapeEm[iene])[12] << ", " + << (m_pulseShapeEm[iene])[13] << ", " + << (m_pulseShapeEm[iene])[14] << ", " + << (m_pulseShapeEm[iene])[15] << ", " + << (m_pulseShapeEm[iene])[16] << ", " + << (m_pulseShapeEm[iene])[17] << ", " + << (m_pulseShapeEm[iene])[18] << ", " + << (m_pulseShapeEm[iene])[19] << ", " + << (m_pulseShapeEm[iene])[20] << ", " + << (m_pulseShapeEm[iene])[21] << ", " + << (m_pulseShapeEm[iene])[22] << ", " + << (m_pulseShapeEm[iene])[23] << ", " + << endreq; + } + msglog << MSG::INFO + << " 1 auto-corr matrix for EMB+EMEC : " + << m_autoCorrEm[0] << ", " + << m_autoCorrEm[1] << ", " + << m_autoCorrEm[2] << ", " + << m_autoCorrEm[3] << ", " + << m_autoCorrEm[4] << ", " + << m_autoCorrEm[5] << ", " + << m_autoCorrEm[6] + << endreq; + + msglog << MSG::DEBUG + << " Finished reading 1 calib coeff + 1 noise value per eta bin for EM: " + << endreq; + for(int ieta=0 ; ieta<s_NBETABINS ; ieta++) { +// currently the calib coeffs are all set to 1 -> turning INFO logs into DEBUG + msglog << MSG::DEBUG + << "Ieta= " << ieta+1 + << ", calib coeff EMB: " << m_calibCoeffEmb[ieta] + << ", calib coeff EMEC: " << m_calibCoeffEmec[ieta] + << endreq; + msglog << MSG::INFO + << "Ieta= " << ieta+1 + << ", noise rms EMB: " << m_noiseRmsEmb[ieta] << " MeV" + << ", noise rms EMEC: " << m_noiseRmsEmec[ieta] << " MeV" + << endreq; + } + + // + // ...... Read auxiliary data file for HEC + // + + pulsedataname=PathResolver::find_file ("LArHecLvl1.data", "DATAPATH"); + if (pulsedataname == "") { + msglog << MSG::ERROR << "Could not locate LArHecLvl1.data file" << endreq; + return StatusCode::FAILURE; + } + pulsedatafile= pulsedataname.c_str() ; + infile.open(pulsedatafile) ; + + if(!infile) { + msglog << MSG::ERROR + << " cannot open HEC file " + << endreq; + return StatusCode::FAILURE; + } else { + msglog << MSG::DEBUG + << " HEC file opened " + << endreq; + } + + // first read the pulse shape for each eta bin (from Leonid Kurchaninov) + m_pulseShapeHec.resize(s_MAXSAMPLES); + infile >> m_pulseShapeHec [0] + >> m_pulseShapeHec [1] + >> m_pulseShapeHec [2] + >> m_pulseShapeHec [3] + >> m_pulseShapeHec [4] + >> m_pulseShapeHec [5] + >> m_pulseShapeHec [6] + >> m_pulseShapeHec [7] + >> m_pulseShapeHec [8] + >> m_pulseShapeHec [9] + >> m_pulseShapeHec [10] + >> m_pulseShapeHec [11] + >> m_pulseShapeHec [12] + >> m_pulseShapeHec [13] + >> m_pulseShapeHec [14] + >> m_pulseShapeHec [15] + >> m_pulseShapeHec [16] + >> m_pulseShapeHec [17] + >> m_pulseShapeHec [18] + >> m_pulseShapeHec [19] + >> m_pulseShapeHec [20] + >> m_pulseShapeHec [21] + >> m_pulseShapeHec [22] + >> m_pulseShapeHec [23] ; + + m_pulseShapeDerHec.resize(s_MAXSAMPLES); + infile >> m_pulseShapeDerHec [0] + >> m_pulseShapeDerHec [1] + >> m_pulseShapeDerHec [2] + >> m_pulseShapeDerHec [3] + >> m_pulseShapeDerHec [4] + >> m_pulseShapeDerHec [5] + >> m_pulseShapeDerHec [6] + >> m_pulseShapeDerHec [7] + >> m_pulseShapeDerHec [8] + >> m_pulseShapeDerHec [9] + >> m_pulseShapeDerHec [10] + >> m_pulseShapeDerHec [11] + >> m_pulseShapeDerHec [12] + >> m_pulseShapeDerHec [13] + >> m_pulseShapeDerHec [14] + >> m_pulseShapeDerHec [15] + >> m_pulseShapeDerHec [16] + >> m_pulseShapeDerHec [17] + >> m_pulseShapeDerHec [18] + >> m_pulseShapeDerHec [19] + >> m_pulseShapeDerHec [20] + >> m_pulseShapeDerHec [21] + >> m_pulseShapeDerHec [22] + >> m_pulseShapeDerHec [23] ; + + m_satEnergyHec.resize(s_NBETABINS); + m_sinThetaHec.resize(s_NBETABINS); + m_noiseRmsHec.resize(s_NBETABINS); + m_autoCorrHec.push_back(autoCorr); + for(int ieta=1 ; ieta<s_NBETABINS ; ieta++) { + // now read for each eta bin: + // - the calibration factor (determined from MC) + // - the transverse saturating energies (same in all eta bins but Ieta=15) + // - the value of sin_theta + // - the noise rms + infile >> Ieta + >> m_satEnergyHec[ieta] + >> m_sinThetaHec[ieta] + >> m_noiseRmsHec[ieta]; + + // and the autocorrelation function for each eta bin (from Leonid Kurchaninov) + infile >> autoCorr[0] + >> autoCorr[1] + >> autoCorr[2] + >> autoCorr[3] + >> autoCorr[4] + >> autoCorr[5] + >> autoCorr[6] ; + m_autoCorrHec.push_back(autoCorr); + } + infile.close() ; + msglog << MSG::DEBUG + << " HEC file closed " + << endreq; + msglog << MSG::INFO + << " 1 pulse shape for HEC : " + << m_pulseShapeHec[0] << ", " + << m_pulseShapeHec[1] << ", " + << m_pulseShapeHec[2] << ", " + << m_pulseShapeHec[3] << ", " + << m_pulseShapeHec[4] << ", " + << m_pulseShapeHec[5] << ", " + << m_pulseShapeHec[5] << ", " + << m_pulseShapeHec[6] << ", " + << m_pulseShapeHec[7] << ", " + << m_pulseShapeHec[8] << ", " + << m_pulseShapeHec[9] << ", " + << m_pulseShapeHec[10] << ", " + << m_pulseShapeHec[11] << ", " + << m_pulseShapeHec[12] << ", " + << m_pulseShapeHec[13] << ", " + << m_pulseShapeHec[14] << ", " + << m_pulseShapeHec[15] << ", " + << m_pulseShapeHec[16] << ", " + << m_pulseShapeHec[17] << ", " + << m_pulseShapeHec[18] << ", " + << m_pulseShapeHec[19] << ", " + << m_pulseShapeHec[20] << ", " + << m_pulseShapeHec[21] << ", " + << m_pulseShapeHec[22] << ", " + << m_pulseShapeHec[23] + << endreq; + msglog << MSG::DEBUG + << "Finished reading calib coeff, noise rms, sat ene, auto corr for each eta bin for HEC: " + << endreq; + for(int ieta=1 ; ieta<s_NBETABINS ; ieta++) { +// currently the calib coeffs are all set to 1 -> turning INFO logs into DEBUG + msglog << MSG::DEBUG + << " Ieta= " << ieta+1 + << ", calib coeff HEC: " << m_calibCoeffHec[ieta] + << endreq; + msglog << MSG::INFO + << " Ieta= " << ieta+1 + << ", noise rms HEC: " << m_noiseRmsHec[ieta] << " MeV" + << ", sat ene HEC: " << m_satEnergyHec[ieta] << " MeV" + << endreq; + msglog << MSG::INFO + << " Ieta= " << ieta+1 + << ", auto corr HEC: " + << (m_autoCorrHec[ieta])[0] << ", " + << (m_autoCorrHec[ieta])[1] << ", " + << (m_autoCorrHec[ieta])[2] << ", " + << (m_autoCorrHec[ieta])[3] << ", " + << (m_autoCorrHec[ieta])[4] << ", " + << (m_autoCorrHec[ieta])[5] << ", " + << (m_autoCorrHec[ieta])[6] << ", " + << endreq; + } + + // + // ...... Read auxiliary data files for FCAL + // + + pulsedataname=PathResolver::find_file ("LArFcalLvl1.data", "DATAPATH"); + if (pulsedataname == "") { + msglog << MSG::ERROR << "Could not locate LArFcalLvl1.data file" << endreq; + return StatusCode::FAILURE; + } + pulsedatafile= pulsedataname.c_str() ; + infile.open(pulsedatafile) ; + + if(!infile) { + msglog << MSG::ERROR + << " cannot open FCAL file " + << endreq; + return StatusCode::FAILURE; + } else { + msglog << MSG::DEBUG + << " FCAL file opened " + << endreq; + } + + const int nMod = 3; + // m_noiseRmsFcal.resize(nMod); + int imod=0; + + for(int iMod=0; iMod < nMod; iMod++) { + + // first read the module number + noise rms (from J. Rutherfoord) + // infile >> imod >> m_noiseRmsFcal[iMod] ; + infile >> imod + >> noiseRms4[0] + >> noiseRms4[1] + >> noiseRms4[2] + >> noiseRms4[3] ; + m_noiseRmsFcal.push_back(noiseRms4); + + // read the pulse shape for this module (from John Rutherfoord) + infile >> pulseShape [0] + >> pulseShape [1] + >> pulseShape [2] + >> pulseShape [3] + >> pulseShape [4] + >> pulseShape [5] + >> pulseShape [6] + >> pulseShape [7] + >> pulseShape [8] + >> pulseShape [9] + >> pulseShape [10] + >> pulseShape [11] + >> pulseShape [12] + >> pulseShape [13] + >> pulseShape [14] + >> pulseShape [15] + >> pulseShape [16] + >> pulseShape [17] + >> pulseShape [18] + >> pulseShape [19] + >> pulseShape [20] + >> pulseShape [21] + >> pulseShape [22] + >> pulseShape [23] ; + m_pulseShapeFcal.push_back(pulseShape); + + // read the pulse shape derivative + infile >> pulseShapeDer [0] + >> pulseShapeDer [1] + >> pulseShapeDer [2] + >> pulseShapeDer [3] + >> pulseShapeDer [4] + >> pulseShapeDer [5] + >> pulseShapeDer [6] + >> pulseShapeDer [7] + >> pulseShapeDer [8] + >> pulseShapeDer [9] + >> pulseShapeDer [10] + >> pulseShapeDer [11] + >> pulseShapeDer [12] + >> pulseShapeDer [13] + >> pulseShapeDer [14] + >> pulseShapeDer [15] + >> pulseShapeDer [16] + >> pulseShapeDer [17] + >> pulseShapeDer [18] + >> pulseShapeDer [19] + >> pulseShapeDer [20] + >> pulseShapeDer [21] + >> pulseShapeDer [22] + >> pulseShapeDer [23] ; + m_pulseShapeDerFcal.push_back(pulseShapeDer); + } + + m_autoCorrFcal.resize(s_NBSAMPLES); + // finally read standard autocorrelation matrix + infile >> m_autoCorrFcal[0] + >> m_autoCorrFcal[1] + >> m_autoCorrFcal[2] + >> m_autoCorrFcal[3] + >> m_autoCorrFcal[4] + >> m_autoCorrFcal[5] + >> m_autoCorrFcal[6] ; + + infile.close() ; + msglog << MSG::DEBUG + << " FCAL file closed " + << endreq; + + std::vector<float> auxV(3); + m_calibCoeffFcal.push_back(m_calibCoeffFcalEm); + m_calibCoeffFcal.push_back(m_calibCoeffFcalHad); + + msglog << MSG::DEBUG + << "Finished reading noise, calib coeff and pulse shape for each module for FCAL: " + << endreq; + for(int iMod=0; iMod < nMod; iMod++) { + + msglog << MSG::INFO + << " iMod= " << iMod + << ", noise rms FCAL (eta bin 1,2,3,4): " + << m_noiseRmsFcal[iMod] << " (transverse) MeV " + << endreq; + if(iMod < 2){ + msglog << MSG::INFO + << " iMod= " << iMod + << ", calib coeff FCAL (eta bin 1,2,3,4): " + << (m_calibCoeffFcal[iMod])[0] << ", " + << (m_calibCoeffFcal[iMod])[1] << ", " + << (m_calibCoeffFcal[iMod])[2] << ", " + << (m_calibCoeffFcal[iMod])[3] + << endreq; + } + msglog << MSG::INFO + << " iMod= " << iMod + << ", pulse shape FCAL: " + << (m_pulseShapeFcal[iMod])[0] << ", " + << (m_pulseShapeFcal[iMod])[1] << ", " + << (m_pulseShapeFcal[iMod])[2] << ", " + << (m_pulseShapeFcal[iMod])[3] << ", " + << (m_pulseShapeFcal[iMod])[4] << ", " + << (m_pulseShapeFcal[iMod])[5] << ", " + << (m_pulseShapeFcal[iMod])[6] << ", " + << (m_pulseShapeFcal[iMod])[7] << ", " + << (m_pulseShapeFcal[iMod])[8] << ", " + << (m_pulseShapeFcal[iMod])[9] << ", " + << (m_pulseShapeFcal[iMod])[10] << ", " + << (m_pulseShapeFcal[iMod])[11] << ", " + << (m_pulseShapeFcal[iMod])[12] << ", " + << (m_pulseShapeFcal[iMod])[13] << ", " + << (m_pulseShapeFcal[iMod])[14] << ", " + << (m_pulseShapeFcal[iMod])[15] << ", " + << (m_pulseShapeFcal[iMod])[16] << ", " + << (m_pulseShapeFcal[iMod])[17] << ", " + << (m_pulseShapeFcal[iMod])[18] << ", " + << (m_pulseShapeFcal[iMod])[19] << ", " + << (m_pulseShapeFcal[iMod])[20] << ", " + << (m_pulseShapeFcal[iMod])[21] << ", " + << (m_pulseShapeFcal[iMod])[22] << ", " + << (m_pulseShapeFcal[iMod])[23] + << endreq; + } + msglog << MSG::INFO + << "auto corr FCAL: " + << m_autoCorrFcal[0] << ", " + << m_autoCorrFcal[1] << ", " + << m_autoCorrFcal[2] << ", " + << m_autoCorrFcal[3] << ", " + << m_autoCorrFcal[4] << ", " + << m_autoCorrFcal[5] << ", " + << m_autoCorrFcal[6] + << endreq; + + // now the relative gains + pulsedataname=PathResolver::find_file ("Fcal_ptweights_table7.data", "DATAPATH"); + if (pulsedataname == "") { + msglog << MSG::ERROR << "Could not locate Fcal_ptweights_table7.data file" << endreq; + return StatusCode::FAILURE; + } + pulsedatafile= pulsedataname.c_str() ; + infile.open(pulsedatafile) ; + + if(!infile) { + msglog << MSG::ERROR + << " cannot open FCAL gains file " + << endreq; + return StatusCode::FAILURE; + } else { + msglog << MSG::DEBUG + << " FCAL gains file opened " + << endreq; + } + + // first read the gain file for all cells (from John Rutherfoord) + + m_cellRelGainFcal.resize(m_fcalHelper->channel_hash_max()); + // file structure: + const unsigned int colNum=14; + // number of cells per 'group' (sharing the same gain) + const unsigned int maxCell = 4; + std::vector<std::string> TTlabel; + TTlabel.resize(colNum); + std::vector<float> gain; + gain.resize(colNum); + int iline=0; + int ngain=0; + + while ( infile >> TTlabel [0] >> gain [0] + >> TTlabel [1] >> gain [1] + >> TTlabel [2] >> gain [2] + >> TTlabel [3] >> gain [3] + >> TTlabel [4] >> gain [4] + >> TTlabel [5] >> gain [5] + >> TTlabel [6] >> gain [6] + >> TTlabel [7] >> gain [7] + >> TTlabel [8] >> gain [8] + >> TTlabel [9] >> gain [9] + >> TTlabel [10] >> gain [10] + >> TTlabel [11] >> gain [11] + >> TTlabel [12] >> gain [12] + >> TTlabel [13] >> gain [13]) { + + msglog << MSG::DEBUG + << " TTlabel[0], gain[0]= " << TTlabel[0] << ", " << gain[0] + << endreq; + msglog << MSG::DEBUG + << " [...] " + << endreq; + msglog << MSG::DEBUG + << " TTlabel[13], gain[13]= " << TTlabel[13] << ", " << gain[13] + << endreq; + + //int barrel_ec_fcal = 2; + //int pos_neg = -2; // C side (neg z) + int detZside = -1; + if(iline < 32) {detZside = 1;} // A side (pos z) + //if(iline < 32) {pos_neg = 2;} // A side (pos z) + for(unsigned int icol=0; icol<colNum; icol++){ + + int em_had = 0; // FCAL1 + int layer=0; + if(icol > 7) { + em_had = 1; // FCAL2+FCAL3 + if(icol > 11) {layer=1;} + } + int region = 3 ; // FCAL + + std::string TTlab = TTlabel[icol]; + //int Ieta = int(TTlab[0])-48; + //int Lphi = int(TTlab[1])-64; + int eta = int(TTlab[0])-49; + int phi = int(TTlab[1])-65; + int nPhi = 16; + if(detZside<0) { + phi = ( phi<nPhi/2 ? nPhi/2 -1 -phi : 3*nPhi/2 -1 -phi ); + } + + int group = 1; + if(TTlab.size() > 3) {group = int(TTlab[3])-48;} + + // added Nov 2009 following introduction of TTCell map "FcalFix" + // modified again Jun 2010 following introduction of TTCell map "FcalFix2" (the C side is not symmetric wrt the A side) + // modified again Jan 2011 in an attempt to correctly load the gains to the correct cells (not all cells were being assigned gains) + if(em_had) { // FCAL-had only + if( (layer==1 && detZside>0) || (layer==0 && detZside<0) ) { + if(eta==0||eta==2) { + eta+=1; + //group+=1; //- OLD CODE (pre Jan2011 change) + + } else { + if(layer==1) {group+=1;} else { group+=2; } + } + } + else { + if(eta==1||eta==3) { + eta-=1; + //group+=2; //- OLD CODE (pre Jan2011 change) + if(layer==1) {group+=1;} else { group+=2; } + } + } + } + // + // ... create an offline TT+layer identifier from the labels read in. + // + Identifier ttId = m_lvl1Helper->layer_id(detZside, em_had, region, eta, phi, layer) ; + msglog << MSG::DEBUG + << "ttId= " << m_lvl1Helper->show_to_string(ttId) + << endreq; + // + //... fill a vector with offline identifier of cells belonging to this tower (with layer info) + // + std::vector<Identifier> cellIdVec = m_ttSvc->createCellIDvecLayer(ttId); + // + // .... loop on all LAr offline channels belonging to the trigger tower (with layer info) + // + std::vector<unsigned int> hashVec; + // number of connected cells in the current tower + unsigned int nCell=0; + for (unsigned int ichan=0;ichan<cellIdVec.size();ichan++){ + Identifier cellId=cellIdVec[ichan]; + + int cellPhi = m_fcalHelper->phi(cellId); + if(cellPhi >= 0) { // protection skipping unconnected channels (they have eta=phi=-999) - normally, they should not be in the list + nCell++; + // + // .... convert to FCAL hash + // + IdentifierHash fcalHash = m_fcalHelper->channel_hash(cellId); + // + //... save hash indices + // + hashVec.push_back(fcalHash); + } // end condition cellPhi>=0 + + } // end of loop on channels in the tower + msglog << MSG::DEBUG + << "nCell= " << nCell + << endreq; + + // + //... loop on the 4 cells of the current group and put their gain in the gain vector (indexed by hashes) + // + for(unsigned int iCell=0; iCell<maxCell; iCell++) { + unsigned int index0 = (group-1)*maxCell + iCell ; + if(index0 < nCell) { // unconnected channels have highest hash ids (because highest eta values) + unsigned int index = hashVec[index0] ; + m_cellRelGainFcal [index] = gain[icol] ; + ngain++; + msglog << MSG::DEBUG + << " index, gain= " << index << ", " << gain[icol] + << endreq; + } + } + + } // end of loop on columns + iline++; + } // end of loop on lines + + + msglog << MSG::DEBUG + << " number of lines found= " + << iline + << endreq; + infile.close() ; + msglog << MSG::INFO + << " FCAL gains file closed, extracted " << ngain << " gains" + << endreq; + + return StatusCode::SUCCESS; + +} +//=========================================================================== +int LArTTL1Maker::decodeInverse(int region, int eta) +{ +//=========================================================================== +// Maps [ region , eta ] to [ Ieta] +// ========================================================================== +// Problem: this is NOT a bijection, because of the "barrel end" zone. +// Convention: only the barrel part of the identifying fields is returned +//=========================================================================== + + int Ieta=0; + if(region == 0){ // Barrel + EC-OW + if(eta <= 14) { + Ieta=eta+1; + } + else { + Ieta=eta-13; + } + } + else if( region == 1 || region == 2) { // EC-IW + if(region == 1) { + Ieta=eta+12; + } + else { + Ieta=15; + } + } + else if(region == 3) { // FCAL + Ieta=eta+1; + } + return Ieta ; +} + + +StatusCode LArTTL1Maker::initHitMap() +{ + + MsgStream msglog(messageService(),name()); + int outputLevel = msgSvc()->outputLevel( name() ); + + if (m_useMapfromStore) { + msglog << MSG::INFO << " Use LArHitEMap from detector store. should be filled by LArDigitMaker (digitmaker1)" << endreq; + return StatusCode::SUCCESS; + } + + m_hitmap = new LArHitEMap(); + +// +// ......... make the Sub-detector flag vector +// + static std::vector<bool> SubDetFlag; + for (int i=0; i < LArHitEMap::NUMDET ; i++) + { + SubDetFlag.push_back(false); + } + +// +// ......... make the LArHit container name list +// + if ( m_SubDetectors == "LAr_All" ) + { + m_HitContainer.push_back(m_EmBarrelHitContainerName); + SubDetFlag[LArHitEMap::EMBARREL_INDEX] = true; + m_HitContainer.push_back(m_EmEndCapHitContainerName); + SubDetFlag[LArHitEMap::EMENDCAP_INDEX] = true; + m_HitContainer.push_back(m_HecHitContainerName); + SubDetFlag[LArHitEMap::HADENDCAP_INDEX] = true; + m_HitContainer.push_back(m_ForWardHitContainerName); + SubDetFlag[LArHitEMap::FORWARD_INDEX] = true; + } + else if ( m_SubDetectors == "LAr_Em" ) + { + m_HitContainer.push_back(m_EmBarrelHitContainerName); + SubDetFlag[LArHitEMap::EMBARREL_INDEX] = true; + m_HitContainer.push_back(m_EmEndCapHitContainerName); + SubDetFlag[LArHitEMap::EMENDCAP_INDEX] = true; + } + else if ( m_SubDetectors == "LAr_EmBarrel" ) + { + m_HitContainer.push_back(m_EmBarrelHitContainerName); + SubDetFlag[LArHitEMap::EMBARREL_INDEX] = true; + } + else if ( m_SubDetectors == "LAr_EmEndCap" ) + { + m_HitContainer.push_back(m_EmEndCapHitContainerName); + SubDetFlag[LArHitEMap::EMENDCAP_INDEX] = true; + } + else if ( m_SubDetectors == "LAr_HEC" ) + { + m_HitContainer.push_back(m_HecHitContainerName); + SubDetFlag[LArHitEMap::HADENDCAP_INDEX] = true; + } + else if ( m_SubDetectors == "LAr_Fcal" ) + { + m_HitContainer.push_back(m_ForWardHitContainerName); + SubDetFlag[LArHitEMap::FORWARD_INDEX] = true; + } + else + { + // + //........ Unknown case + // + msglog << MSG::FATAL + << "Invalid SubDetector property : " << m_SubDetectors + << "Valid ones are: LAr_All, LAr_Em, LAr_EmBarrel, " + << "LAr_EmEndCap, LAr_HEC and LAr_Fcal " + << endreq; + + return(StatusCode::FAILURE); + } + +// +// put all cell energies to 0. +// + + if ( ! m_hitmap->Initialize(SubDetFlag) ) + { + msglog << MSG::FATAL + << "Making of the cell table failed" + << endreq; + return StatusCode::FAILURE; + } + + if (outputLevel <= MSG::DEBUG) { + msglog << MSG::INFO + << "Number of created cells in LArHitEMap " + << m_hitmap->GetNbCells() + << endreq; + } + + + return StatusCode::SUCCESS; +} diff --git a/LArCalorimeter/LArL1Sim/src/components/LArL1Sim_entries.cxx b/LArCalorimeter/LArL1Sim/src/components/LArL1Sim_entries.cxx new file mode 100755 index 00000000000..12d72f169a2 --- /dev/null +++ b/LArCalorimeter/LArL1Sim/src/components/LArL1Sim_entries.cxx @@ -0,0 +1,23 @@ +// +// Entry file for ATHENA component libraries +// +#include "GaudiKernel/DeclareFactoryEntries.h" + +#include "LArL1Sim/LArTTL1Maker.h" +#include "LArL1Sim/LArSCL1Maker.h" +#include "LArL1Sim/LArTTL1Calib.h" +#include "LArL1Sim/LArSCSimpleMaker.h" + +DECLARE_ALGORITHM_FACTORY( LArTTL1Maker ) +DECLARE_ALGORITHM_FACTORY( LArSCL1Maker ) +DECLARE_ALGORITHM_FACTORY( LArTTL1Calib ) +DECLARE_ALGORITHM_FACTORY( LArSCSimpleMaker ) + +DECLARE_FACTORY_ENTRIES(LArL1Sim) { + + DECLARE_ALGORITHM(LArTTL1Maker) + DECLARE_ALGORITHM(LArSCL1Maker) + DECLARE_ALGORITHM(LArTTL1Calib) + DECLARE_ALGORITHM(LArSCSimpleMaker) + +} diff --git a/LArCalorimeter/LArL1Sim/src/components/LArL1Sim_load.cxx b/LArCalorimeter/LArL1Sim/src/components/LArL1Sim_load.cxx new file mode 100755 index 00000000000..f9afc8aceed --- /dev/null +++ b/LArCalorimeter/LArL1Sim/src/components/LArL1Sim_load.cxx @@ -0,0 +1,8 @@ +// +// Load file for ATHENA component library +// + + +#include "GaudiKernel/LoadFactoryEntries.h" + +LOAD_FACTORY_ENTRIES(LArL1Sim) -- GitLab