diff --git a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArCalibPatchingAlg.icc b/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArCalibPatchingAlg.icc index 8322373a5411a669dc7173b86527400d41e11191..2e650436131d7444168a80fdf4c0fb93bf5de7af 100644 --- a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArCalibPatchingAlg.icc +++ b/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArCalibPatchingAlg.icc @@ -647,7 +647,7 @@ bool LArCalibPatchingAlg<CONDITIONSCONTAINER>::getAverage(const HWIdentifier chi nCells++; } for(idx=0; idx<patch.size(); ++idx) { - patch[idx] = LArCaliWave( patch_wave[idx].getWave(), patch[idx].getDt(), patch[idx].getDAC(), patch[idx].getFlag()); + patch[idx] = LArCaliWave( patch_wave[idx].getWave(), patch[idx].getDt(), patch[idx].getDAC(), patch[idx].getIsPulsedInt(), patch[idx].getFlag()); } } if (nCells==0) { diff --git a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMFitterTool.h b/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMFitterTool.h deleted file mode 100644 index 7ca574e910639b06a73681dcf56447c211e1e20a..0000000000000000000000000000000000000000 --- a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMFitterTool.h +++ /dev/null @@ -1,95 +0,0 @@ -//Dear emacs, this is -*- c++ -*- - -/* - Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration -*/ - - -/** - * @author D. Prieur prieur@lapp.in2p3.fr - * 5. 4. 2004 - */ - -#ifndef LArTCMFitterTool_H -#define LArTCMFitterTool_H - -#include "AthenaBaseComps/AthAlgTool.h" - -#include "LArRawConditions/LArCaliWave.h" -#include "LArRawConditions/LArPhysWave.h" - -class TMinuit; -class TCanvas; -class TH1D; - -static const InterfaceID IID_LArTCMFitterTool("LArTCMFitterTool", 1 , 0); - -class LArTCMFitterTool : public AthAlgTool -{ - -public: - - // Retrieve interface ID - static const InterfaceID& interfaceID() { return IID_LArTCMFitterTool; } - - - // constructor - LArTCMFitterTool(const std::string& type, const std::string& name, const IInterface* parent ); - virtual ~LArTCMFitterTool(); - - virtual StatusCode initialize(); - virtual StatusCode finalize(){return StatusCode::SUCCESS;} - - LArPhysWave Fit(const LArCaliWave&, const LArPhysWave&); - void fcn(int &, double *, double &, double *, int); - - double getTaud() const { return m_taud; }; - double getTauexp() const { return m_tauexp; }; - double getTaur() const { return m_taur; }; - double getf() const { return m_f; }; - double getW0() const { return m_w0; }; - double getAlpha() const { return m_alpha; }; - - double getcaliStart() const { return m_caliStart; }; - double getphysShift() const { return m_physShift; }; - - double getEmean() const { return m_Emean; }; - double getMphyMcal() const { return m_MphyMcal; }; - double getchi2() const { return m_chi2; }; - - void setminuitoutputlevel(int level) {m_minuitoutputlevel = level; }; - -private: - - TMinuit* m_minuit; - - //int m_nparams; - int m_minuitoutputlevel; - - LArCaliWave m_gCali; - LArPhysWave m_gPhys; - - //int m_iGain, m_iLayer, m_iPhi, m_iEta, m_dac, m_ipdg; //not necessary - - double m_taud, m_f, m_tauexp, m_taur, m_w0, m_alpha; - double m_tauderr, m_ferr, m_tauexperr, m_taurerr, m_w0err, m_alphaerr; - - double m_caliStart, m_caliStarterr; - double m_physShift, m_physShifterr; - - double m_Emean, m_MphyMcal, m_chi2; - - double computeChi2(const LArPhysWave &, const LArPhysWave &); - - - unsigned int getStart(const LArCaliWave &); - unsigned int getStart(const LArWave &); - unsigned int getindMax(const LArWave &); - - TH1D LArCaliWave2TH1D(const LArCaliWave &); - TH1D LArPhysWave2TH1D(const LArPhysWave &); - TH1D LArWave2TH1D(const LArWave &); - -}; - -#endif diff --git a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMPhysWavePredictor.h b/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMPhysWavePredictor.h deleted file mode 100644 index 94f1fa1e443fe7b13c6ed03fabae58e51e605dda..0000000000000000000000000000000000000000 --- a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMPhysWavePredictor.h +++ /dev/null @@ -1,70 +0,0 @@ -//Dear emacs, this is -*- c++ -*- - -/* - Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration -*/ - - -/** - * @author D. Prieur prieur@lapp.in2p3.fr - * 5. 4. 2004 - */ - -#ifndef LARTCMPHYSWAVEPREDICTOR -#define LARTCMPHYSWAVEPREDICTOR - -#include "AthenaBaseComps/AthAlgorithm.h" -#include "GaudiKernel/ToolHandle.h" - -#include "LArRawConditions/LArWFParams.h" -#include "LArRawConditions/LArPhysWave.h" -#include "LArRawConditions/LArCaliWave.h" - -#include "LArCalibUtils/LArWFParamTool.h" -#include "LArCalibUtils/LArPhysWaveTool.h" -#include "LArCabling/LArOnOffIdMapping.h" -#include "StoreGate/ReadCondHandleKey.h" - -#include <vector> -#include <string> - -class LArTCMPhysWavePredictor : public AthAlgorithm -{ - public: - LArTCMPhysWavePredictor(const std::string & name, ISvcLocator * pSvcLocator); - - ~LArTCMPhysWavePredictor(); - - //standard algorithm methods - StatusCode initialize() ; - StatusCode execute() {return StatusCode::SUCCESS;} //empty method - StatusCode stop(); - StatusCode finalize(){ATH_CHECK(m_cablingKey.initialize()); return StatusCode::SUCCESS;} - private: - - SG::ReadCondHandleKey<LArOnOffIdMapping> m_cablingKey{this, "OnOffMap", "LArOnOffIdMap", "SG key for mapping object"}; - - bool m_testmode; - bool m_datafromfile; //switch to take data from file/db - int m_minuitoutputlevel; - bool m_rootrawdump; - std::string m_rootoutputfile; - - bool m_filter_cells; - int m_filter_layer; - - bool m_filter_phi; - int m_filter_phi_min; - int m_filter_phi_max; - - bool m_filter_eta; - int m_filter_eta_min; - int m_filter_eta_max; - - - std::vector<std::string> m_keylist; - std::string m_groupingType; - -}; - -#endif diff --git a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMPredictor.h b/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMPredictor.h deleted file mode 100755 index 32c4e2792cde541e58480485357c8211d3893a4e..0000000000000000000000000000000000000000 --- a/LArCalorimeter/LArCalibUtils/LArCalibUtils/LArTCMPredictor.h +++ /dev/null @@ -1,87 +0,0 @@ -//Dear emacs, this is -*- c++ -*- - -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - - -/** - * @author D. Prieur prieur@lapp.in2p3.fr - * 5. 4. 2004 - */ - - -#ifndef LARTCMPREDICTOR_H -#define LARTCMPREDICTOR_H - -#include "GaudiKernel/AlgTool.h" -#include "GaudiKernel/MsgStream.h" - -#include "LArRawConditions/LArCaliWave.h" -#include "LArRawConditions/LArPhysWave.h" - -//class TMinuit; - -class LArTCMPredictor -{ -public: - // constructor - LArTCMPredictor(); - virtual ~LArTCMPredictor(); - - - LArPhysWave predictLArPhysWave(const LArCaliWave &); - - double getTaud() const { return m_taud; }; - double getTauexp() const { return m_tauexp; }; - double getTaur() const { return m_taur; }; - double getf() const { return m_f; }; - double getW0() const { return m_w0; }; - double getAlpha() const { return m_alpha; }; - - double getcaliStart() const { return m_caliStart; }; - double getphysShift() const { return m_physShift; }; - - void setTaud(double taud) { m_taud=taud; }; - void setTauexp(double tauexp) { m_tauexp=tauexp; }; - void setTaur(double taur) { m_taur=taur; }; - void setf(double f) { m_f=f; }; - void setW0(double w0) { m_w0=w0; }; - void setAlpha(double alpha) { m_alpha=alpha; }; - - void setcaliStart(double caliStart) { m_caliStart=caliStart; }; - void setphysShift(double physShift) { m_physShift=physShift; }; - -// template <class T> unsigned int getMax(const T &PULSE); - double getMax(const LArCaliWave &larCaliWave); - double getMax(const LArPhysWave &larPhysWave); -private: - - LArCaliWave m_gCali; - LArPhysWave m_gPhysPred; - - double m_caliStart; //, m_caliStarterr; - double m_physShift; //, m_physShifterr; - - double m_taud, m_f, m_tauexp, m_taur, m_w0, m_alpha; - - double evalLArCaliWave(int); - double Integrand(int, int); - double qfast(int , int ); - - double linearInterpolate(const LArWave&, double); - LArWave translate(const LArWave&, int); - unsigned int getindMax(const LArWave &); - - - unsigned int getStart(const LArCaliWave &); - unsigned int getStart(const LArWave &); - - - - double Gt2(double t) const; - double Gt3(double t) const; - -}; - -#endif diff --git a/LArCalorimeter/LArCalibUtils/src/LArCaliWaveAverage.cxx b/LArCalorimeter/LArCalibUtils/src/LArCaliWaveAverage.cxx index 7d91a3e43dabac15f79745ce295aa4a454ef12ec..180a61fb807f0b049181d421173e4c529d77f6ca 100644 --- a/LArCalorimeter/LArCalibUtils/src/LArCaliWaveAverage.cxx +++ b/LArCalorimeter/LArCalibUtils/src/LArCaliWaveAverage.cxx @@ -346,7 +346,8 @@ LArCaliWave LArCaliWaveAverage::WaveAverage(std::vector<LArCaliWave> ToBeAverage ToBeAveraged[0].getErrors(), ToBeAveraged[0].getTriggers(), ToBeAveraged[0].getDt(), - (ToBeAveraged[0].getDAC() + (ToBeAveraged[0].getIsPulsedInt()<<16)), + ToBeAveraged[0].getDAC(), + ToBeAveraged[0].getIsPulsedInt(), ToBeAveraged[0].getFlag() ); return theCaliWaveAverage; diff --git a/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilder.cxx b/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilder.cxx index 10e2883c03bf0478e14f32c5b3c59f670596fc93..6c86d6eac92e100aad98a02c7a7ed899f408d3b0 100644 --- a/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilder.cxx +++ b/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilder.cxx @@ -253,8 +253,7 @@ StatusCode LArCaliWaveBuilder::executeWithAccumulatedDigits() WaveMap::iterator itm = waveMap.find(index); if ( itm == waveMap.end() ) { // A new LArCaliWave is booked - LArCaliWave wave(samplesum.size()*m_NStep, m_dt, dac, pulsed); - wave.setFlag( LArWave::meas ); + LArCaliWave wave(samplesum.size()*m_NStep, m_dt, dac, pulsed, LArWave::meas); itm = (waveMap.insert(WaveMap::value_type(index,wave))).first; ATH_MSG_DEBUG("index: "<<index<<" new wave inserted"); } @@ -317,8 +316,7 @@ StatusCode LArCaliWaveBuilder::executeWithStandardDigits() WaveMap::iterator itm = waveMap.find((*it)->DAC()); if ( itm == waveMap.end() ) { // A new LArCaliWave is booked - LArCaliWave wave(samples.size()*m_NStep, m_dt, (*it)->DAC(),0x1); - wave.setFlag( LArWave::meas ); + LArCaliWave wave(samples.size()*m_NStep, m_dt, (*it)->DAC(), 0x1, LArWave::meas ); itm = (waveMap.insert(WaveMap::value_type((*it)->DAC(),wave))).first; } @@ -460,11 +458,12 @@ StatusCode LArCaliWaveBuilder::stop() continue ; } else { dacWaves.emplace_back(((thisWave)+(-pedAve)).getWave() , - thisWave.getErrors(), - thisWave.getTriggers(), - thisWave.getDt(), - (thisWave.getDAC() + (thisWave.getIsPulsedInt()<<24)), - thisWave.getFlag() ); + thisWave.getErrors(), + thisWave.getTriggers(), + thisWave.getDt(), + thisWave.getDAC(), + thisWave.getIsPulsedInt(), + thisWave.getFlag() ); NCaliWave++; } diff --git a/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilderXtalk.cxx b/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilderXtalk.cxx index 652088b6c2c2e2bb631534254942da91508cf440..b01627bbc9046a55544a8f1e931e78a57f565f6f 100644 --- a/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilderXtalk.cxx +++ b/LArCalorimeter/LArCalibUtils/src/LArCaliWaveBuilderXtalk.cxx @@ -541,8 +541,7 @@ StatusCode LArCaliWaveBuilderXtalk::execute() WaveMap::iterator itm = waveMap.find((*it)->DAC()); if ( itm == waveMap.end() ) { // A new LArCaliWave is booked - LArCaliWave wave(samplesum.size()*m_NStep, m_dt, (*it)->DAC(), 0x1); - wave.setFlag( LArWave::meas ) ; + LArCaliWave wave(samplesum.size()*m_NStep, m_dt, (*it)->DAC(), 0x1, LArWave::meas ) ; // Be careful! Don't call addAccumulatedEvent() twice!!! // cf. LArCaliwaveBuilder ref. 1.38 @@ -664,7 +663,8 @@ StatusCode LArCaliWaveBuilderXtalk::stop() thisWave.getErrors(), thisWave.getTriggers(), thisWave.getDt(), - (thisWave.getDAC() + (thisWave.getIsPulsedInt()<<24)), + thisWave.getDAC(), + thisWave.getIsPulsedInt(), thisWave.getFlag() ); dacWaves.push_back(newWave); diff --git a/LArCalorimeter/LArCalibUtils/src/LArDeltaRespTool.cxx b/LArCalorimeter/LArCalibUtils/src/LArDeltaRespTool.cxx index a0aa3da1e5c6e99d1c2a6b5b222c00b240238a15..8a1660cae8b7242225dc32a590db8df406aa9feb 100755 --- a/LArCalorimeter/LArCalibUtils/src/LArDeltaRespTool.cxx +++ b/LArCalorimeter/LArCalibUtils/src/LArDeltaRespTool.cxx @@ -95,6 +95,7 @@ void LArDeltaRespTool::compute_deltaresp() m_gDelta = LArCaliWave( gDelta.getWave() , m_gCali.getDt() , deltaRespDAC , + m_gCali.getIsPulsedInt(), m_gCali.getFlag() ); } diff --git a/LArCalorimeter/LArCalibUtils/src/LArMasterWaveBuilder.cxx b/LArCalorimeter/LArCalibUtils/src/LArMasterWaveBuilder.cxx index 38d325c303944deffcc5cbe89628a83d58a7ddf9..a6fc81f16ddb5528df10e55f3981509ecf4dbcbb 100644 --- a/LArCalorimeter/LArCalibUtils/src/LArMasterWaveBuilder.cxx +++ b/LArCalorimeter/LArCalibUtils/src/LArMasterWaveBuilder.cxx @@ -442,9 +442,9 @@ StatusCode LArMasterWaveBuilder::stop() LArCaliWaveContainer::LArCaliWaves& dacWaves = mwContainer->get(chID, gain_it); LArCaliWave masterWave( fitWave[1].getWave(), - dt, -1, LArWave::mwf ); + dt, -1, 0x1, LArWave::mwf ); LArCaliWave dac0Wave ( fitWave[0].getWave(), - dt, -2, LArWave::dac0 ); + dt, -2, 0x1, LArWave::dac0 ); dacWaves.push_back( masterWave ); dacWaves.push_back( dac0Wave ); nMasterWaves ++ ; diff --git a/LArCalorimeter/LArCalibUtils/src/LArPhysWavePredictor.cxx b/LArCalorimeter/LArCalibUtils/src/LArPhysWavePredictor.cxx index 84a0fa2b45ec0f39a1056e6c27ca69f41f357346..9ca4b43926deb557200b66abb06e02dda9a95df6 100644 --- a/LArCalorimeter/LArCalibUtils/src/LArPhysWavePredictor.cxx +++ b/LArCalorimeter/LArCalibUtils/src/LArPhysWavePredictor.cxx @@ -540,7 +540,8 @@ StatusCode LArPhysWavePredictor::stop() theLArCaliWave.getErrors(), theLArCaliWave.getTriggers(), theLArCaliWave.getDt(), - theLArCaliWave.getDAC(), + theLArCaliWave.getDAC(), + theLArCaliWave.getIsPulsedInt(), theLArCaliWave.getFlag() ); } } diff --git a/LArCalorimeter/LArCalibUtils/src/LArTCMFitterTool.cxx b/LArCalorimeter/LArCalibUtils/src/LArTCMFitterTool.cxx deleted file mode 100644 index c0f57aefb3090900ba1ed2ebd0d56eeb05a2ca54..0000000000000000000000000000000000000000 --- a/LArCalorimeter/LArCalibUtils/src/LArTCMFitterTool.cxx +++ /dev/null @@ -1,309 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#include "LArCalibUtils/LArTCMFitterTool.h" -#include "LArCalibUtils/LArTCMPredictor.h" - -#include "TROOT.h" -#include "TApplication.h" -#include "TSystem.h" - -#include <TMinuit.h> - -#include <TH1.h> - -//must be renAMED TO LArTCMFitterToolTool.cxx - -static LArTCMFitterTool* static_LArTCMFitterTool_pointer = NULL; - -void LArTCMFitterTool_fcn_wrapper(int &npar, double *gin, double &f, double *par, int iflag) { - if (static_LArTCMFitterTool_pointer==NULL) { - std::cout << "Ugly things happening here..." << std::endl; - exit(-1); - } - static_LArTCMFitterTool_pointer->fcn(npar, gin, f, par, iflag); -} - - -LArTCMFitterTool::LArTCMFitterTool ( const std::string& type, const std::string& name,const IInterface* parent ) - : AthAlgTool(type,name,parent), - m_minuit(0), - m_minuitoutputlevel(0), - m_taud(0), - m_f(0), - m_tauexp(0), - m_taur(0), - m_w0(0), - m_alpha(0), - m_tauderr(0), - m_ferr(0), - m_tauexperr(0), - m_taurerr(0), - m_w0err(0), - m_alphaerr(0), - m_caliStart(0), - m_caliStarterr(0), - m_physShift(0), - m_physShifterr(0), - m_Emean(0), - m_MphyMcal(0), - m_chi2(0) -{ -} - -// destructor -LArTCMFitterTool::~LArTCMFitterTool() -{ -} - -///////////////////////////////////// -StatusCode LArTCMFitterTool::initialize() -{ - if( static_LArTCMFitterTool_pointer != NULL ) { - ATH_MSG_ERROR ( "Problem with static pointer to the data object used by minuits FCN-function!!!" ); - return StatusCode::FAILURE; - } - - //m_nparams = 8; -// m_minuit = new TMinuit(8); //initialize TMinuit with a maximum of nparams parameters -// static_LArTCMFitterTool_pointer = this; -// m_minuit->SetFCN(LArTCMFitterTool_fcn_wrapper); -// log<<MSG::DEBUG<<"Minuit and FCN initialized"<<endmsg; - - //double arglist[10]; - //int ierflg = 0; - //m_minuit->SetPrintLevel(1); // -1: no output, 1: std output - //m_minuit->mnexcm("SET NOW", arglist, 0, ierflg); // no warning - - m_minuitoutputlevel = -1; - gSystem->Load("libMinuit"); - - return StatusCode::SUCCESS; -} - -//////////////////////////////////////////////////////////////////////////////////////////////// -LArPhysWave LArTCMFitterTool::Fit(const LArCaliWave &larCaliWave, const LArPhysWave &larPhysWave) -{ - m_gCali = larCaliWave; - m_gPhys = larPhysWave; - - m_minuit = new TMinuit(8); //initialize TMinuit with a maximum of nparams parameters - static_LArTCMFitterTool_pointer = this; - m_minuit->SetFCN(LArTCMFitterTool_fcn_wrapper); - m_minuit->SetPrintLevel(m_minuitoutputlevel); // -1: no output, 1: std output - ATH_MSG_INFO("Minuit and FCN initialized"); - - - // Find beginning of the caliWave - m_caliStart = (double)getStart(m_gCali)-2.; // this need to be tuned - if (m_caliStart<0.) m_caliStart=0.; // - - - // Set initial starting values and step sizes for the 8 parameters - Double_t vstart[8] = {450., 340. , 0.065, 0.184, 1., m_caliStart, 27, 0.025}; - Double_t step[8] = {1., 1. , 0.001 , 0.001 , 0.1, 0.5, 1., 1e-3}; - - double arglist[10]; - int ierflg = 0; - - m_minuit->mnparm(0, "taud", vstart[0], step[0], 50, 750., ierflg); - m_minuit->mnparm(1, "tauexp", vstart[1], step[1], 300, 800., ierflg); - m_minuit->mnparm(2, "f", vstart[2], step[2], 0.0001, 0.9, ierflg); - m_minuit->mnparm(3, "w0", vstart[3], step[3], 0.1e-3, 1., ierflg); - m_minuit->mnparm(4, "taur", vstart[4], step[4], 1e-5, 50., ierflg); - m_minuit->mnparm(5, "caliStart", vstart[5], step[5], 0., 80., ierflg); - m_minuit->mnparm(6, "physShift", vstart[6], step[6], 0., 80., ierflg); - m_minuit->mnparm(7, "alpha", vstart[7], step[7], 1e-6, 1., ierflg); - - // Fixed parameters - // taud, tauexp, f, caliStart and alpha are fixed - // w0, taur and physShift are left free - arglist[0]= 1; // taud - arglist[1]= 2; // tauexp - arglist[2]= 3; // f - arglist[3]= 6; // caliStart - arglist[4]= 8; // alpha - m_minuit->mnexcm("FIX",arglist,5,ierflg); - - // Scan of parameters - arglist[0] = 7; // scan of physShift - m_minuit->mnexcm("SCAN",arglist,1,ierflg); - - // Set minuit run - arglist[0] = 1; - m_minuit->mnexcm("SET ERR", arglist ,1,ierflg); - arglist[0] = 0; // 0 fast / 1 standart / 2 improve - m_minuit->mnexcm("SET STR", arglist ,1,ierflg); - - // Start minimization method - arglist[0] = 1e11; - arglist[1] = 0.1; - //m_minuit->mnexcm("MIGRAD", arglist ,2,ierflg); - //m_minuit->mnexcm("SIMPLEX", arglist ,2,ierflg); - //m_minuit->mnexcm("MIGRAD", arglist ,2,ierflg); - m_minuit->mnexcm("MINIMIZE", arglist ,2,ierflg); // launch minimization - - // Print fit results - double amin, edm, errdef; - int nvpar, nparx, icstat; - m_minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat); - m_minuit->mnprin(3,amin); - - // Get results - m_minuit->GetParameter(0,m_taud,m_tauderr); - m_minuit->GetParameter(1,m_tauexp,m_tauexperr); - m_minuit->GetParameter(2,m_f,m_ferr); - m_minuit->GetParameter(3,m_w0,m_w0err); - m_minuit->GetParameter(4,m_taur,m_taurerr); - m_minuit->GetParameter(5,m_caliStart,m_caliStarterr); - m_minuit->GetParameter(6,m_physShift,m_physShifterr); - m_minuit->GetParameter(7,m_alpha,m_alphaerr); - m_chi2=amin; - - ATH_MSG_INFO("Final set of parameters:"); - ATH_MSG_INFO("taud: "<<m_taud<<" tauexp: "<<m_tauexp<<" f: "<<m_f<<" w0: "<<m_w0<<" taur: "<<m_taur<<" caliStart: "<<m_caliStart<<" phyShift: "<<m_physShift); - - delete m_minuit; - - // re-produce LArPhysWave from final fit parameters - LArTCMPredictor larTCMPredictor; - larTCMPredictor.setTaud(m_taud); - larTCMPredictor.setTauexp(m_tauexp); - larTCMPredictor.setf(m_f); - larTCMPredictor.setW0(m_w0); - larTCMPredictor.setTaur(m_taur); - larTCMPredictor.setcaliStart(m_caliStart); - larTCMPredictor.setphysShift(m_physShift); - larTCMPredictor.setAlpha(m_alpha); - LArPhysWave predLArPhysWave = larTCMPredictor.predictLArPhysWave(larCaliWave); - m_chi2 = computeChi2(m_gPhys,predLArPhysWave); // fill also m_Emean - - //Compute the ratio MphyMcal - double Mcal = larTCMPredictor.getMax(larCaliWave); - double Mphy = larTCMPredictor.getMax(predLArPhysWave); - - if(Mcal!=0) m_MphyMcal = Mphy/Mcal; - else m_MphyMcal=-999; - ATH_MSG_DEBUG("MphyMcal: "<<m_MphyMcal); - - //need to normalize and scale predLArPhysWave by MphyMcal - const unsigned int N = predLArPhysWave.getSize(); - for (unsigned int i=0; i<N; ++i) { - predLArPhysWave.setSample(i,predLArPhysWave.getSample(i)/Mcal); //the wave max is normalized to Mphy/Mcal value - needed for OFC - //predLArPhysWave.setSample(i,predLArPhysWave.getSample(i)); // this line must be used for debug only !!!! - } - - return predLArPhysWave; -} - - -void LArTCMFitterTool::fcn(int &/*npar*/, double */*gin*/, double &f, double *par, int /*iflag*/) { - - LArTCMPredictor larTCMPredictor; - - // Set parameters for LArPhysWave prediction - larTCMPredictor.setTaud(par[0]); - larTCMPredictor.setTauexp(par[1]); - larTCMPredictor.setf(par[2]); - larTCMPredictor.setW0(par[3]); - larTCMPredictor.setTaur(par[4]); - larTCMPredictor.setcaliStart(par[5]); - larTCMPredictor.setphysShift(par[6]); - larTCMPredictor.setAlpha(par[7]); - - // LArPhysWave prediction - LArPhysWave larpredPhysWave = larTCMPredictor.predictLArPhysWave(m_gCali); - - f = computeChi2(m_gPhys,larpredPhysWave); -} - - -double LArTCMFitterTool::computeChi2(const LArPhysWave &larPhysWave, const LArPhysWave &larpredPhysWave) { - - //const Int_t nbins = 160; //fit only on the nbins first bin - // const Int_t nbins = 175; - const int nbins = (larPhysWave.getSize()<larpredPhysWave.getSize())? larPhysWave.getSize() : larpredPhysWave.getSize(); - //std::cout<<"nbins: "<<nbins<<std::endl; - double chisq = 0.; - double S_yf = 0.; - double S_ff = 0.; - double S_yy = 0.; - double y,z; - - // Fit range - int start_fit = 0; - int end_fit = nbins; - - // Compute chi2 - for (Int_t i=start_fit;i<end_fit; ++i) { - y = larPhysWave.getSample(i); - z = larpredPhysWave.getSample(i); - S_yf += y*z; - S_ff += z*z; - S_yy += y*y; - } - if (S_ff==0) S_ff=0.00000001; - - //chi2 et Emean - chisq = S_yy - S_yf*S_yf/S_ff; - m_Emean = S_yf/S_ff; - return chisq; -} - - -unsigned int LArTCMFitterTool::getStart(const LArCaliWave &larCaliWave) -{ - LArWave larWave(larCaliWave.getWave(),larCaliWave.getDt()); - return getStart(larWave); -} - - -unsigned int LArTCMFitterTool::getStart(const LArWave &larWave) -{ - const unsigned int N = larWave.getSize(); - double max=larWave.getSample(getindMax(larWave)); - double treshold = max*0.01; //need to implement member/propertie for treshold - unsigned int k=0; - do k++; while((larWave.getSample(k)<treshold)&&(k<N)); - return k-1; -} - - -unsigned int LArTCMFitterTool::getindMax(const LArWave &larWave) -{ - const unsigned int N = larWave.getSize(); - unsigned int indmax=0; - for (unsigned int i=0; i<N; ++i) { - if(larWave.getSample(i)>larWave.getSample(indmax)) indmax=i; - //if(larWave[i]>larWave[indmax]) indmax=i; //do not work... - } - //std::cout<<"indmax"<<indmax<<std::endl; - return indmax; -} - - -TH1D LArTCMFitterTool::LArCaliWave2TH1D(const LArCaliWave &larCaliWave) -{ - LArWave larWave(larCaliWave.getWave(),larCaliWave.getDt()); - return LArWave2TH1D(larWave); -} - -TH1D LArTCMFitterTool::LArPhysWave2TH1D(const LArPhysWave &larPhysWave) -{ - LArWave larWave(larPhysWave.getWave(),larPhysWave.getDt()); - return LArWave2TH1D(larWave); -} - -TH1D LArTCMFitterTool::LArWave2TH1D(const LArWave &larWave) -{ - const unsigned int N = larWave.getSize(); - TH1D hist("hist","hist",N,0,N); - for(unsigned int i=0; i<N; ++i) { - hist.SetBinContent(i+1,larWave.getSample(i)); - } - return hist; -} - - - diff --git a/LArCalorimeter/LArCalibUtils/src/LArTCMPhysWavePredictor.cxx b/LArCalorimeter/LArCalibUtils/src/LArTCMPhysWavePredictor.cxx deleted file mode 100644 index d2a701ffd07f15da006a85e074b1be45e6c52906..0000000000000000000000000000000000000000 --- a/LArCalorimeter/LArCalibUtils/src/LArTCMPhysWavePredictor.cxx +++ /dev/null @@ -1,546 +0,0 @@ -/* - Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration -*/ - -#include "LArCalibUtils/LArTCMPhysWavePredictor.h" -#include "LArCalibUtils/LArTCMFitterTool.h" - -#include "GaudiKernel/ToolHandle.h" - -#include "LArRawConditions/LArCaliWaveContainer.h" -#include "LArRawConditions/LArPhysWaveContainer.h" -#include "CaloIdentifier/CaloCell_ID.h" - -#include <iostream> -#include <fstream> -#include "TFile.h" -#include "TH1.h" -#include "TH2.h" -#include "TProfile.h" -#include "TTree.h" - -typedef LArPhysWaveContainer::ConstConditionsMapIterator PhysWaveIt; - -LArTCMPhysWavePredictor::LArTCMPhysWavePredictor (const std::string& name, ISvcLocator* pSvcLocator) : - AthAlgorithm(name, pSvcLocator) -{ - declareProperty("TestMode",m_testmode=false); - - declareProperty("DataFromFile",m_datafromfile=false); - declareProperty("MinuitOutputLevel",m_minuitoutputlevel=1); // -1: no output, 1: std output - declareProperty("RootRawDump",m_rootrawdump=false); // to dump caliwave,physwave and fit result to root file. - declareProperty("RootOutputFile",m_rootoutputfile=std::string("dumpTCM.root")) ; - - //Cell filters - declareProperty("FilterCells",m_filter_cells=false); - - declareProperty("FilterCells_Layer",m_filter_layer=-1); - - declareProperty("FilterCells_Phi",m_filter_phi=false); - declareProperty("FilterCells_Phi_Min",m_filter_phi_min=-1); - declareProperty("FilterCells_Phi_Max",m_filter_phi_max=3000); //fancy values... - - declareProperty("FilterCells_Eta",m_filter_eta=false); - declareProperty("FilterCells_Eta_Min",m_filter_eta_min=-3000); //fancy values... - declareProperty("FilterCells_Eta_Max",m_filter_eta_max=3000); - - // list of gains to be processed - m_keylist.clear() ; - m_keylist.push_back("HIGH"); - m_keylist.push_back("MEDIUM"); - m_keylist.push_back("LOW"); - declareProperty("KeyList",m_keylist); - - declareProperty("GroupingType", m_groupingType); - -} - -/////////////////////////////////////////////////// -LArTCMPhysWavePredictor::~LArTCMPhysWavePredictor() -{} - -//////////////////////////////////////////////// -StatusCode LArTCMPhysWavePredictor::initialize() -{ - ATH_MSG_INFO ( " Initialize..."); - return StatusCode::SUCCESS ; -} - -StatusCode LArTCMPhysWavePredictor::stop() -{ - ATH_MSG_INFO ( "From stop..."); - - ToolHandle<LArWFParamTool> larWFParamTool("LArWFParamTool"); - ATH_CHECK( larWFParamTool.retrieve() ); - - ToolHandle<LArTCMFitterTool> larTCMFitterTool("LArTCMFitterTool"); - ATH_CHECK( larTCMFitterTool.retrieve() ); - - larTCMFitterTool->setminuitoutputlevel(m_minuitoutputlevel); - - - const CaloCell_ID* idHelper = nullptr; - ATH_CHECK( detStore()->retrieve (idHelper, "CaloCell_ID") ); - const LArEM_ID* emId = idHelper->em_idHelper(); - if (!emId) { - ATH_MSG_ERROR ( "Could not access lar EM ID helper" ); - return StatusCode::FAILURE; - } - - SG::ReadCondHandle<LArOnOffIdMapping> cablingHdl{m_cablingKey}; - const LArOnOffIdMapping* cabling{*cablingHdl}; - if(!cabling) { - ATH_MSG_ERROR("Do not have mapping object " << m_cablingKey.key()); - return StatusCode::FAILURE; - } - - //open ouput file for dumping waves - TFile f(m_rootoutputfile.c_str(),"recreate"); //need to implement: if (m_rootrawdump) - TTree* tree = 0; - - // TTree structure - struct tree_struct { - Int_t iGain; - Int_t iLayer; - Int_t iPhi; - Int_t iEta; - Int_t ipdg; - - Double_t taud; - Double_t tauexp; - Double_t f; - Double_t w0; - Double_t taur; - Double_t SplineShift; - Double_t top; - Double_t Emean; - Double_t MphyMcal; - Double_t chi2; - - Int_t DeadCell; - Int_t fit_result; - - Double_t peak_residus; - - Int_t CAL_nb; - Int_t CAL_chan; - }; - tree_struct ts1{}; - - const Int_t nlayer = 4; - const Int_t ngain = 3; - const Int_t nsampling = 4; - const Int_t iphilimit[2][nsampling] = {{0,0,0,0},{4,4,16,16}}; //verifier layer 0 !!! - const Int_t ietalimit[2][nsampling] = {{0,0,0,0},{57,451,57,27}}; //verifier layer 0 !!! - - Char_t temp[100]; - TH2F* CellMap[nlayer][ngain]; - //TH2F* NentriesMap[nlayer][ngain]; - - if ( m_rootrawdump ) { - - //create histogram - for(Int_t layer=0;layer<nlayer;layer++){ - for(Int_t gain=0;gain<ngain;gain++) { - - //fill CellMap - memset(temp,0,sizeof(temp)); - sprintf(temp,"CellMaplayer%.1d_gain%.1d",layer,gain); - CellMap[layer][gain] = new TH2F(temp,temp,ietalimit[1][layer],ietalimit[0][layer],ietalimit[1][layer],iphilimit[1][layer],iphilimit[0][layer],iphilimit[1][layer]); - - //fill NentriesMap - //memset(temp,0,sizeof(temp)); - //sprintf(temp,"NentriesMaplayer%.1d_gain%.1d",layer,gain); - //NentriesMap[layer][gain] = new TH2F(temp,temp,ietalimit[1][layer],ietalimit[0][layer],ietalimit[1][layer],iphilimit[1][layer],iphilimit[0][layer],iphilimit[1][layer]); - } - } - - - tree = new TTree("tree","tree essai"); - tree->Branch("iGain",&ts1.iGain,"iGain/I"); - tree->Branch("iLayer",&ts1.iLayer,"iLayer/I"); - tree->Branch("iPhi",&ts1.iPhi,"iPhi/I"); - tree->Branch("iEta",&ts1.iEta,"iEta/I"); - tree->Branch("ipdg",&ts1.ipdg,"ipdg/I"); - - tree->Branch("taud",&ts1.taud,"taud/D"); - tree->Branch("tauexp",&ts1.tauexp,"tauexp/D"); - tree->Branch("f",&ts1.f,"f/D"); - tree->Branch("w0",&ts1.w0,"w0/D"); - tree->Branch("taur",&ts1.taur,"taur/D"); - tree->Branch("SplineShift",&ts1.SplineShift,"SplineShift/D"); - tree->Branch("top",&ts1.top,"top/D"); - tree->Branch("Emean",&ts1.Emean,"Emean/D"); - tree->Branch("MphyMcal",&ts1.MphyMcal,"MphyMcal/D"); - tree->Branch("chi2",&ts1.chi2,"chi2/D"); - - tree->Branch("DeadCell",&ts1.DeadCell,"DeadCell/I"); - tree->Branch("fit_result",&ts1.fit_result,"fit_result/I"); - - tree->Branch("peak_residus",&ts1.peak_residus,"peak_residus/D"); - tree->Branch("CAL_nb",&ts1.CAL_nb,"CAL_nb/I"); - tree->Branch("CAL_chan",&ts1.CAL_chan,"CAL_chan/I"); - } - - //retrieve Physics pulses either from database or from file - // Files are self root files produced by emtb for debug purpose - LArPhysWaveContainer* physWaveContainer = nullptr; - if (!m_datafromfile) { - // retrieve Physics pulses from database - const LArPhysWaveContainer* constphysWaveContainer = nullptr; - ATH_CHECK( detStore()->retrieve(constphysWaveContainer) ); - physWaveContainer = (LArPhysWaveContainer*) constphysWaveContainer; - ATH_MSG_INFO ( "LArPhysWaveContainer found in StoreGate." ); - } else { - //from self external file; to debug and compare with emtb - // for the moment just test P13 phi=9 layer 2 - physWaveContainer = new LArPhysWaveContainer(); - ATH_CHECK( physWaveContainer->setGroupingType(m_groupingType,msg()) ); - - ATH_CHECK( physWaveContainer->initialize() ); - - TFile physicsinput("/afs/cern.ch/user/p/prieur/scratch0/OF-P13_21-PhyPredDataFile.root"); - physicsinput.cd("PhysicsPulses"); - - int ilayer=2; - int iphi = 9; - int ieta_start = 25; - int ieta_end = 35; - int igain; - if ((ilayer==1)||(ilayer==2)) igain = CaloGain::LARMEDIUMGAIN; - else igain = CaloGain::LARHIGHGAIN; - - for(int ieta = ieta_start;ieta<ieta_end;ieta++) { - char name[80]; - memset(name,0,sizeof(name)); - sprintf(name,"hphys%.1d%.1d%.2d%.3d",igain,ilayer,iphi,ieta); - TProfile* physicspulse=(TProfile*)gDirectory->Get(name); - if(physicspulse==NULL) { - std::cout<<"Could not get physics pulse "<<name<<" from file."<<std::endl; - continue; - } - - Identifier id=emId->channel_id(+1,ilayer,+0,ieta,iphi); - HWIdentifier chID = cabling->createSignalChannelID(id); - //if ((ilayer==1)||(ilayer==2)) larPhysWave = new LArPhysWave(175,1,1); //gain 1 for debug - //else larPhysWave = new LArPhysWave(175,1,1); - LArPhysWave larPhysWave (175, 1, 1); - ATH_MSG_DEBUG("name: "<<name<<" ilayer: "<<ilayer<<" ieta: "<<ieta<<" iphi: "<<iphi<<" id "<< chID.get_compact()); - - // decode id to check if correct - //HWIdentifier testchID = larPhysWave->channelID(); - //Identifier testid=cabling->cnvToIdentifier(testchID); - //std::cout <<"physwave "<<" eta "<<emId->eta(testid)<<" phi "<< emId->phi(testid)<<" layer "<<emId->sampling(testid)<<" gain "<<larPhysWave->getGain()<<std::endl; - - //Fill LArPhysWave - for(int i=0;i<physicspulse->GetNbinsX();i++) { - larPhysWave.setSample(i,physicspulse->GetBinContent(i+1)); - //std::cout<<physicspulse->GetBinContent(i+1)<<" "<<larPhysWave->getSample(i)<<std::endl; - } - - //physWaveContainer->push_back(larPhysWave); // old container - physWaveContainer->setPdata(chID, larPhysWave, igain); // new container - - } - - //log<<MSG::INFO<<"LArPhysWaveContainer size: "<<physWaveContainer->size()<<endmsg; - - } // end else m_datafromfile - - unsigned nchannel = 0 ; - - //Create LArPhysWaveContainer for Predicted Pulses - LArPhysWaveContainer* larpredPhysWaveContainer = new LArPhysWaveContainer(); - ATH_CHECK( larpredPhysWaveContainer->setGroupingType(m_groupingType,msg()) ); - - ATH_CHECK( larpredPhysWaveContainer->initialize() ); - - //Create container for Params -// LArWFParamsContainer* larWFParamsContainer=new LArWFParamsContainer(); - - // get the LArCaliWave from the detector store - std::vector<std::string>::const_iterator key_it=m_keylist.begin(); - std::vector<std::string>::const_iterator key_it_e=m_keylist.end(); - - for (;key_it!=key_it_e;++key_it) { // Loop over all LArCaliWave containers that are to be processed from m_keylist - - std::string keyCali = *key_it ; - - const LArCaliWaveContainer* caliWaveContainer = nullptr; - ATH_CHECK( detStore()->retrieve(caliWaveContainer,keyCali) ); - - ATH_MSG_INFO ( "Processing LArCaliWaveContainer from StoreGate, key = " << keyCali ); - - for ( unsigned gain = CaloGain::LARHIGHGAIN ; gain < CaloGain::LARNGAIN ; ++gain ) { // Loop over gains in current container - - ATH_MSG_VERBOSE ( "Now processing gain = " << gain << " in LArCaliWaveContainer with key = " << keyCali ); - - // loop over cali wave container - typedef LArCaliWaveContainer::ConstConditionsMapIterator const_iterator; - const_iterator itVec = caliWaveContainer->begin(gain); - const_iterator itVec_e = caliWaveContainer->end(gain); - - if ( itVec == itVec_e ) { - ATH_MSG_INFO ( "LArCaliWaveContainer (key=" << keyCali << ") has no wave with gain = " << gain ); - continue ; - } - - for (; itVec != itVec_e; ++itVec) { // loop over cells for the current gain - - LArCaliWaveContainer::LArCaliWaves::const_iterator cont_it = (*itVec).begin(); - LArCaliWaveContainer::LArCaliWaves::const_iterator cont_it_e = (*itVec).end(); - - if ( cont_it == cont_it_e ) { - ATH_MSG_DEBUG ( "Empty channel found in LArCaliWave container: skipping..." ); - continue ; - } else { - ATH_MSG_DEBUG ( (*itVec).size() << " LArCaliWaves found for channel 0x" << MSG::hex << itVec.channelId() << MSG::dec ); - } - - const HWIdentifier chID = itVec.channelId(); - Identifier id ; - - try { - id = cabling->cnvToIdentifier(chID); - } catch ( const LArID_Exception& ) { - ATH_MSG_ERROR ( "Cabling exception caught for channel 0x" << MSG::hex << chID << MSG::dec ); - } - int layer = emId->sampling(id); - int eta = emId->eta(id); - int phi = emId->phi(id); - - if ( m_filter_cells ) { //Filter Cells - - if ( m_filter_layer!=-1 && layer!=m_filter_layer ) { - ATH_MSG_INFO ( "Skipping channel 0x" << MSG::hex << chID << MSG::dec << " because of FilterCells_Layer cut."); - continue; - } - - if ( m_filter_eta ) - if ( (eta<m_filter_eta_min) || (eta>m_filter_eta_max) ) { - ATH_MSG_INFO ( "Skipping channel 0x" << MSG::hex << chID << MSG::dec << " because of FilterCells_Eta cut."); - continue; - } - - if ( m_filter_phi ) { - - if ( (layer==1 && ( phi<m_filter_phi_min/4 || phi>m_filter_phi_max/4) ) || - (layer!=1 && ( phi<m_filter_phi_min || phi>m_filter_phi_max ) ) - ) { - ATH_MSG_INFO ( "Skipping channel 0x" << MSG::hex << chID << MSG::dec << " because of FilterCells_Phi cut."); - continue; - } - - } - - } // end filter cells - - - ATH_MSG_VERBOSE( "Channel 0x" << MSG::hex << chID << MSG::dec - << " (Layer = " << layer - << " - Eta = " << eta - << " - Phi = " << phi - << " - Gain = " << gain - << ") passed filtering selection. Now processing..."); - - // Always use the caliwave with higher dac value for each gain. - // Predicted waves are fitted to gain 0 physics waves for sampling 0 - // and to gain 1 physics waves for sampling 1,2,3 - int idac=0, physgain=0; - if ((layer==0)&&(gain==0)) { idac=5000; physgain=0; } //use intermediate dac for avoiding saturation - was 10000 in emtb - if ((layer==1)&&(gain==0)) { idac=500; physgain=1; } - if ((layer==1)&&(gain==1)) { idac=8000; physgain=1; } - if ((layer==2)&&(gain==0)) { idac=500; physgain=1; } //fit seems netter with MG phys data - if ((layer==2)&&(gain==1)) { idac=5000; physgain=1; } - if ((layer==3)&&(gain==0)) { idac=500; physgain=0; } - if ((layer==3)&&(gain==1)) { idac=3000; physgain=0; } //no phys data for back in MG - - ATH_MSG_VERBOSE ( "... will use LArCaliWave with DAC = " << idac << " and LArPhysWave with Gain = " << physgain ); - - unsigned ndac = 0 ; - - for (;cont_it!=cont_it_e;++cont_it) { // Loop over available DAC wavwsvalues for the given cell - - const LArCaliWave& larCaliWave = (*cont_it); - int dac = larCaliWave.getDAC() ; - - ATH_MSG_VERBOSE ( ++ndac << ". DAC = " << dac ); - - if ( dac!=idac ) continue; // skip unwanted DAC waves - - ATH_MSG_INFO ( " Processing channel 0x" << MSG::hex << chID << MSG::dec - << " -> Layer = " << layer - << " - Eta = " << eta - << " - Phi = " << phi - << " - Gain = " << gain - << " - LArCaliWave DAC = " << dac ); - - // get corresponding physcs wave - const LArPhysWave& larPhysWave = physWaveContainer->get(chID,physgain); - - if ( larPhysWave.isEmpty() ) { - - ATH_MSG_ERROR ( "No LArPhysWave found for channel 0x" << MSG::hex << itVec.channelId() << MSG::dec ); - - } else { // LArPhysWave found, perform TCM fit and phys wave prediction from fit params - - //log << MSG::INFO << "Now processing channel " << ++nchannel << " with TCM fit " << endmsg ; - - ATH_MSG_DEBUG ( "larPhysWave gain = " << gain ); - ATH_MSG_DEBUG ( "larPhysWave size = " << larPhysWave.getSize() ); - ATH_MSG_DEBUG ( "larPhysWave flag = " << larPhysWave.getFlag() ); - - LArPhysWave predlarPhysWave = - larTCMFitterTool->Fit(larCaliWave,larPhysWave); // <--- - - predlarPhysWave.setFlag(LArWave::predFitPhys); - - ATH_MSG_INFO ( " End of TCM fit" ); - - //output fit values - ATH_MSG_DEBUG("Fit values from larTCMFitterTool:"); - ATH_MSG_DEBUG(larTCMFitterTool->getTaud()<<" "<<larTCMFitterTool->getTauexp()<<" "<<larTCMFitterTool->getTaur() - <<" "<<larTCMFitterTool->getf()<<" "<<larTCMFitterTool->getW0()); - ATH_MSG_DEBUG(larTCMFitterTool->getcaliStart()<<" "<<larTCMFitterTool->getphysShift()); - ATH_MSG_DEBUG(larTCMFitterTool->getEmean()<<" "<<larTCMFitterTool->getMphyMcal()<<" "<<larTCMFitterTool->getchi2()); - - //Store fit parameters in LArWFParams -#if 0 - LArWFParams* larwfparams = new LArWFParams((unsigned int)LArWave::predFitPhys, - larTCMFitterTool->getTauexp(), - larTCMFitterTool->getf(), - larTCMFitterTool->getTaud(), - larTCMFitterTool->getW0(), - larTCMFitterTool->getTaur(), - -1.00, - larTCMFitterTool->getMphyMcal()); -#endif - - //Put predicted LArPhysWave in LArPhysWaveContainer - larpredPhysWaveContainer->setPdata(itVec.channelId(), predlarPhysWave, gain); // new container - -// //Put waveform params in LArWFParamsContainer -// larWFParamsContainer->push_back(larwfparams); - - - //for debug - dump LArWave to root file - if (m_rootrawdump) { - - unsigned int N; - TH1D* hist; - char temp[80]; - - N = larPhysWave.getSize(); - - hist=new TH1D("hist","hist",N,0,N); - for(unsigned int i=0; i<N; ++i) { - hist->SetBinContent(i+1,larPhysWave.getSample(i)); - } - memset(temp,0,sizeof(temp)); - sprintf(temp,"phys%.1u%.1d%.2d%.3d", gain,emId->sampling(id),emId->phi(id),emId->eta(id)); - hist->SetName(temp); - hist->SetTitle(temp); - f.cd(); - hist->Write(); - delete hist; - - N = larCaliWave.getSize(); - hist=new TH1D("hist","hist",N,0,N); - for(unsigned int i=0; i<N; ++i) { - hist->SetBinContent(i+1,larCaliWave.getSample(i)); - } - int Mcal=static_cast<int> (hist->GetMaximum()); - memset(temp,0,sizeof(temp)); - sprintf(temp,"cali%.1u%.1d%.2d%.3d", gain,emId->sampling(id),emId->phi(id),emId->eta(id)); - hist->SetName(temp); - hist->SetTitle(temp); - f.cd(); - hist->Write(); - delete hist; - - N = predlarPhysWave.getSize(); - hist=new TH1D("hist","hist",N,0,N); - for(unsigned int i=0; i<N; ++i) { - hist->SetBinContent(i+1,predlarPhysWave.getSample(i)*Mcal*larTCMFitterTool->getEmean()); //change normalization to match with delaypulse - //hist->SetBinContent(i+1,predlarPhysWave.getSample(i)); //same normalisation to what is stored in the DB - the wave is already normalized to MphyMcal value - needed for OFC - } - memset(temp,0,sizeof(temp)); - sprintf(temp,"pphys%.1u%.1d%.2d%.3d", gain,emId->sampling(id),emId->phi(id),emId->eta(id)); - hist->SetName(temp); - hist->SetTitle(temp); - f.cd(); - hist->Write(); - delete hist; - - CellMap[emId->sampling(id)][gain]->SetBinContent(emId->eta(id)+1,emId->phi(id)+1,1.); - - //Fill TTree - ts1.iGain = gain; - ts1.iLayer = emId->sampling(id); - ts1.iPhi = emId->phi(id); - ts1.iEta = emId->eta(id); - ts1.ipdg = -1; - - ts1.taud = larTCMFitterTool->getTaud(); - ts1.tauexp = larTCMFitterTool->getTauexp(); - ts1.f = larTCMFitterTool->getf(); - ts1.w0 = larTCMFitterTool->getW0(); - ts1.taur = larTCMFitterTool->getTaur(); - ts1.SplineShift = larTCMFitterTool->getcaliStart(); - ts1.top = larTCMFitterTool->getphysShift(); - ts1.Emean = larTCMFitterTool->getEmean(); - ts1.chi2 = larTCMFitterTool->getchi2(); - ts1.MphyMcal = larTCMFitterTool->getMphyMcal(); - - ts1.fit_result = -1; - - ts1.peak_residus = -1; - ts1.CAL_nb = -1; - ts1.CAL_chan = -1; - - tree->Fill(); - - } //end of rootrawdump - - } // end if check if LArPhysWave->isEmpty() - - } // end loop over DAC values - - if ( m_testmode && nchannel >= 1 ) { - ATH_MSG_DEBUG ( "Test mode selected: only 1 channel processed." ); - break ; - } - - } // end loop over cells - - } // end loop over gains - - } // end loop over LArCaliWave containers (keys in m_keylist) - - std::string keyout ; - - // Record larpredPhysWaveContainer to detector store - keyout = "LArPhysWavesTCM" ; - ATH_CHECK( detStore()->record(larpredPhysWaveContainer,keyout) ); - ATH_MSG_INFO ( "LArPhysWaveContainer has been recorded to StoreGate with key="<<keyout); - - // Record larWFParamsContainer to DetectorStore -// keyout = "LArWFParamsTCM" ; -// sc=detStore->record(larWFParamsContainer,keyout ); -// if (sc.isFailure()) { -// log << MSG::FATAL << "Cannot record LArWFParamsContainer to StoreGate! key=" << keyout << endmsg; -// return StatusCode::FAILURE; -// } else log << MSG::INFO << "LArWFParamsContainer has been recorded to StoreGate with key="<<keyout<<endmsg; - - if (m_rootrawdump) { - for(Int_t layer=0;layer<nlayer;layer++){ - for(Int_t gain=0;gain<ngain;gain++) { - CellMap[layer][gain]->Write(); - //NentriesMap[layer][gain]->Write(); - } - } - tree->Write(); - } - - return StatusCode::SUCCESS; -} diff --git a/LArCalorimeter/LArCalibUtils/src/LArTCMPredictor.cxx b/LArCalorimeter/LArCalibUtils/src/LArTCMPredictor.cxx deleted file mode 100755 index 973cdb2888ca0b59841926d86d0a323568fbfda3..0000000000000000000000000000000000000000 --- a/LArCalorimeter/LArCalibUtils/src/LArTCMPredictor.cxx +++ /dev/null @@ -1,225 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -#include "LArCalibUtils/LArTCMPredictor.h" -#include <cmath> -//#include "LArCalibUtils/TMinuit.h" - - -LArTCMPredictor::LArTCMPredictor () -{ - m_caliStart=-999.; - m_physShift=0.; //not used so far - - //Characteristic start values - // need to declareproperty - m_taud = 450.; - m_f = 0.09; - m_tauexp = 365.; - m_taur=11.; - m_w0 = 0.184; - m_alpha=0.025; -} - -// destructor -LArTCMPredictor::~LArTCMPredictor() -{ -} - - -// implement eval pour double x et interpolation - -double LArTCMPredictor::evalLArCaliWave(int x) { - return m_gCali.getSample(x); // use shifted LArCaliWave -} - - -double LArTCMPredictor::Integrand(int t, int tau) { - return evalLArCaliWave(tau)*Gt3(t-tau); //need to choose between model -} - -double LArTCMPredictor::qfast(int a, int b) { - double sum=0.; - for(int k=a; k<=b; ++k) { //basic computation of the convolution - sum+=Integrand(b,k); - } - return sum; -} - -/////////////////////////////////////////////////////////////////////////////// -LArPhysWave LArTCMPredictor::predictLArPhysWave(const LArCaliWave &larCaliWave) -{ - m_gCali = larCaliWave; - - LArWave gCali(m_gCali), gPhysPred, gPhysPredShifted; - - if (m_caliStart==-999.) { // if m_caliStart not initialized by LArTCMFitter - m_caliStart = (double)getStart(m_gCali)-2.; // - if (m_caliStart<0.) m_caliStart=0.; // for low dac the beginning of the pulse is not always well defined - } - - gCali = translate(gCali,(int)m_caliStart); // shift of gCali - m_gCali = LArCaliWave(gCali.getWave(),m_gCali.getDt(),m_gCali.getDAC(), 0x1); //create shifted LArCaliWave - - gPhysPred.setSize(larCaliWave.getSize()); - for(unsigned int i=0; i<m_gCali.getSize();++i) { - gPhysPred.setSample(i,qfast(0,0+i)); //choose between method - } - - //Translation of m_physShift - gPhysPredShifted.setSize(larCaliWave.getSize()); - for(unsigned int i=0; i<gPhysPredShifted.getSize();++i) { - if((double)i>m_physShift) { - gPhysPredShifted.setSample(i,linearInterpolate(gPhysPred,(double)i-m_physShift)); - } else gPhysPredShifted.setSample(i,0.); - } - - //temporary - - m_gPhysPred = LArPhysWave(gPhysPredShifted.getWave(),m_gCali.getDt()); - return m_gPhysPred; -} - -// implement translate avec double nshift et interpolation - - -double LArTCMPredictor::linearInterpolate(const LArWave& larWave, double t) -{ - int bin = (int)t; - double x = t - (double)bin; - return (larWave.getSample(bin+1)-larWave.getSample(bin))*x+larWave.getSample(bin); -} - -LArWave LArTCMPredictor::translate(const LArWave& larWave, int nshift) -{ - const unsigned int N = larWave.getSize(); - std::vector<double> trWave; - trWave.resize(N); - for(unsigned int i=0; i<N; ++i) { - if(i<(N-nshift)) trWave[i]=larWave.getSample(i+nshift); - else trWave[i]=0.; - } - return LArWave(trWave,larWave.getDt()); - -} - -unsigned int LArTCMPredictor::getindMax(const LArWave &larWave) -{ - const unsigned int N = larWave.getSize(); - unsigned int indmax=0; - for (unsigned int i=0; i<N; ++i) { - if(larWave.getSample(i)>larWave.getSample(indmax)) indmax=i; - //if(larWave[i]>larWave[indmax]) indmax=i; //do not work... - } - //std::cout<<"indmax"<<indmax<<std::endl; - return indmax; -} - -/* -template <class T> unsigned int LArTCMPredictor::getMax(const T &PULSE) -{ - LArWave larWave(PULSE.getWave(),PULSE.getDt()); - return PULSE.getSample(getindMax(larWave)); -} -*/ - -double LArTCMPredictor::getMax(const LArCaliWave &larCaliWave) -{ - LArWave larWave(larCaliWave.getWave(),larCaliWave.getDt()); - return larCaliWave.getSample(getindMax(larWave)); -} - -double LArTCMPredictor::getMax(const LArPhysWave &larPhysWave) -{ - LArWave larWave(larPhysWave.getWave(),larPhysWave.getDt()); - return larPhysWave.getSample(getindMax(larWave)); -} - - -unsigned int LArTCMPredictor::getStart(const LArCaliWave &larCaliWave) -{ - LArWave larWave(larCaliWave.getWave(),larCaliWave.getDt()); - return getStart(larWave); -} - - -unsigned int LArTCMPredictor::getStart(const LArWave &larWave) -{ - const unsigned int N = larWave.getSize(); - double max=larWave.getSample(getindMax(larWave)); - double treshold = max*0.01; //need to implement member/propertie for treshold - unsigned int k=0; - do k++; while((larWave.getSample(k)<treshold)&&(k<N)); - return k-1; -} - -//////////////////////////////////////// -double LArTCMPredictor::Gt2(double t) const -// Function from Fares Djema - model without dirac nor rC -// Should not be used -{ -double term1, termexp, termcos, termsin; -double w0_2 = m_w0*m_w0; - -term1 = -1./(m_f*m_taud); -termexp = (-m_f*m_tauexp*w0_2+m_tauexp*w0_2-(m_w0*m_tauexp)*(m_w0*m_tauexp)/m_taud +(m_w0*m_tauexp)*(m_w0*m_tauexp)/(m_taud*m_f))*exp(-m_f*t/m_tauexp); -termcos = (m_f*m_tauexp*w0_2-m_tauexp*w0_2+(m_w0*m_tauexp)*(m_w0*m_tauexp)/m_taud +m_f/m_taud)*cos(m_w0*t); -termsin = (w0_2*m_w0*m_tauexp*m_tauexp+m_w0-m_f*m_w0*m_tauexp/m_taud+m_tauexp*m_w0/m_taud)*sin(m_w0*t); - -return term1 + (termexp+termcos+termsin)/(m_f*m_f+w0_2*m_tauexp*m_tauexp); -//return -1; -} - -//////////////////////////////////////// -double LArTCMPredictor::Gt3(double t) const -// modele sans dirac avec rC -// -{ - -double term1, term2; - -double w0_2 = m_w0*m_w0; -double taur_2 = m_taur*m_taur; -double tauexp_2 = m_tauexp*m_tauexp; -double f_2 = m_f*m_f; - -double Sqd = sqrt(fabs(-4+taur_2*w0_2)); -double f_t_tauexp=(m_f*t)/m_tauexp; -double tauexp_taur_w0_2 = m_tauexp*m_taur*w0_2; -double t_w0_Sqd = (t*m_w0*Sqd)/2.; -double t_taur_w0_2 = (t*m_taur*w0_2)/2.; -double exp_f_t_tauexp = exp(f_t_tauexp); - -if ((-4+taur_2*w0_2)>0.) { - term1 = -(exp(-(f_t_tauexp) - t_taur_w0_2)* - (exp(t_taur_w0_2)*Sqd* - ((-1 + m_f)*m_tauexp*(m_f*m_taud + m_tauexp)*w0_2 + - exp_f_t_tauexp*(f_2 + tauexp_2*w0_2 - m_f*tauexp_taur_w0_2)) + - exp_f_t_tauexp*m_f*Sqd* - (m_tauexp*(m_taud - m_tauexp + m_taur)*w0_2 - m_f*(1 + m_taud*m_tauexp*w0_2))* - cosh(t_w0_Sqd) + - exp_f_t_tauexp*m_f*m_w0*(m_tauexp*(-2 - tauexp_taur_w0_2 + taur_2*w0_2 + - m_taud*(-2*m_tauexp + m_taur)*w0_2) + m_f*(2*m_tauexp - m_taur + m_taud*(-2 + tauexp_taur_w0_2)))* - sinh(t_w0_Sqd))); - term2 = (m_f*m_taud*(f_2 + tauexp_2*w0_2 - m_f*tauexp_taur_w0_2)*Sqd); - -} -else { - term1 = -(exp(-(f_t_tauexp) - t_taur_w0_2)* - (exp(t_taur_w0_2)*Sqd* - ((-1 + m_f)*m_tauexp*(m_f*m_taud + m_tauexp)*w0_2 + - exp_f_t_tauexp*(f_2 + tauexp_2*w0_2 - m_f*tauexp_taur_w0_2)) + - exp_f_t_tauexp*m_f*Sqd* - (m_tauexp*(m_taud - m_tauexp + m_taur)*w0_2 - m_f*(1 + m_taud*m_tauexp*w0_2))* - cos(t_w0_Sqd) + - exp_f_t_tauexp*m_f*m_w0*(m_tauexp*(-2 - tauexp_taur_w0_2 + taur_2*w0_2 + - m_taud*(-2*m_tauexp + m_taur)*w0_2) + m_f*(2*m_tauexp - m_taur + m_taud*(-2 + tauexp_taur_w0_2)))* - sin(t_w0_Sqd))); - term2 = (m_f*m_taud*(f_2 + tauexp_2*w0_2 - m_f*tauexp_taur_w0_2)*Sqd); - -} - if (term2==0) {std::cout<<"Gt3::term2 zero divide"<<std::endl;term2=1e-9;} -return term1/term2; - -} diff --git a/LArCalorimeter/LArCalibUtils/src/LArWFParamTool.cxx b/LArCalorimeter/LArCalibUtils/src/LArWFParamTool.cxx index b58de3e23f028287550e12801abb10073dec39eb..0f601cc3661e68421dab77dab90194b9addbdf73 100644 --- a/LArCalorimeter/LArCalibUtils/src/LArWFParamTool.cxx +++ b/LArCalorimeter/LArCalibUtils/src/LArWFParamTool.cxx @@ -258,7 +258,7 @@ StatusCode LArWFParamTool::getLArWaveParams(const LArCaliWave& larCaliWave, const unsigned layer=m_emId->sampling(id); if ( m_storeResOscill[ layer ] && resOscill0) { LArWave injres0 = injRespRes(gCali,wfParams.omega0(),0); - *resOscill0 = LArCaliWave(injres0.getWave(),gCali.getDt(),gCali.getDAC(), 0x1); + *resOscill0 = LArCaliWave(injres0.getWave(),gCali.getDt(),gCali.getDAC(), 0x1,LArWave::unknown); } // find m_Taur using RTM @@ -277,7 +277,7 @@ StatusCode LArWFParamTool::getLArWaveParams(const LArCaliWave& larCaliWave, if ( m_storeResOscill[ layer ] && resOscill1) { LArWave injres1 = injRespRes(gCali,wfParams.omega0(),wfParams.taur()); - *resOscill1 = LArCaliWave(injres1.getWave(),gCali.getDt(),gCali.getDAC(), 0x1); + *resOscill1 = LArCaliWave(injres1.getWave(),gCali.getDt(),gCali.getDAC(), 0x1, LArWave::unknown); } return( StatusCode::SUCCESS ); diff --git a/LArCalorimeter/LArCalibUtils/src/components/LArCalibUtils_entries.cxx b/LArCalorimeter/LArCalibUtils/src/components/LArCalibUtils_entries.cxx index 15042ac1f1e52329d4cb32a5214d4c16996926f6..24da18ab783ac6d697593dd3728c612cb5685611 100644 --- a/LArCalorimeter/LArCalibUtils/src/components/LArCalibUtils_entries.cxx +++ b/LArCalorimeter/LArCalibUtils/src/components/LArCalibUtils_entries.cxx @@ -21,8 +21,6 @@ #include "../LArRampBuilder.h" #include "LArCalibUtils/LArRTMParamExtractor.h" #include "LArCalibUtils/LArStripsCrossTalkCorrector.h" -#include "LArCalibUtils/LArTCMFitterTool.h" -#include "LArCalibUtils/LArTCMPhysWavePredictor.h" #include "LArCalibUtils/LArPhysWavePredictor.h" #include "LArCalibUtils/LArTimeTuning.h" #include "LArCalibUtils/LArTimeTuningNtuple.h" @@ -76,7 +74,6 @@ DECLARE_COMPONENT( LArPhysWaveBuilder ) DECLARE_COMPONENT( LArPhysWaveShifter ) DECLARE_COMPONENT( LArRampBuilder ) DECLARE_COMPONENT( LArRTMParamExtractor ) -DECLARE_COMPONENT( LArTCMPhysWavePredictor ) DECLARE_COMPONENT( LArStripsCrossTalkCorrector ) DECLARE_COMPONENT( LArPhysWavePredictor ) DECLARE_COMPONENT( LArTimeTuning ) @@ -104,7 +101,6 @@ DECLARE_COMPONENT( LArDigitOscillationCorrTool ) DECLARE_COMPONENT( LArDeltaRespTool ) DECLARE_COMPONENT( LArPhysWaveTool ) DECLARE_COMPONENT( LArPhysWaveHECTool ) -DECLARE_COMPONENT( LArTCMFitterTool ) DECLARE_COMPONENT( LArWFParamTool ) DECLARE_COMPONENT( LArPhaseToolConst ) DECLARE_COMPONENT( LArPhaseToolMC ) diff --git a/LArCalorimeter/LArCnv/LArCondAthenaPool/test/LArCaliWaveContainerCnv_test.cxx b/LArCalorimeter/LArCnv/LArCondAthenaPool/test/LArCaliWaveContainerCnv_test.cxx index c23b22c186d6eaea69a8e2209d99a97c4536c7f6..5861825e8f5ccdfe174d245a9511e07b913716f4 100644 --- a/LArCalorimeter/LArCnv/LArCondAthenaPool/test/LArCaliWaveContainerCnv_test.cxx +++ b/LArCalorimeter/LArCnv/LArCondAthenaPool/test/LArCaliWaveContainerCnv_test.cxx @@ -79,7 +79,8 @@ LArCaliWave makeWave (int x) std::vector<int> { 30+x, 31+x, 32+x, 33+x, 34+x }, 1.5 + x, 0x12345678 + x, - x*2); + x*2, + LArWave::unknown); } diff --git a/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p1_test.cxx b/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p1_test.cxx index 8e830d3dea73af56480711c4a6344aea71568f09..5c34a1b511a398711987f1bb1ee2fcce5439bd0e 100644 --- a/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p1_test.cxx +++ b/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p1_test.cxx @@ -127,7 +127,8 @@ LArCaliWave makeWave (int x) std::vector<int> { 30+x, 31+x, 32+x, 33+x, 34+x }, 1.5 + x, 0x12345678 + x, - x*2); + x*2, + LArWave::unknown); } diff --git a/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p2_test.cxx b/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p2_test.cxx index aa9a8a257e64375c03cbf80b267e8351b4215d0c..d07840aa5f453ce6d3501614c02a04dea5731548 100644 --- a/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p2_test.cxx +++ b/LArCalorimeter/LArCnv/LArCondTPCnv/test/LArCaliWaveSubsetCnv_p2_test.cxx @@ -118,7 +118,8 @@ LArCaliWave makeWave (int x) std::vector<int> { 30+x, 31+x, 32+x, 33+x, 34+x }, 1.5 + x, 0x12345678 + x, - x*2); + x*2, + LArWave::unknown); } diff --git a/LArCalorimeter/LArRawConditions/LArRawConditions/LArCaliWave.h b/LArCalorimeter/LArRawConditions/LArRawConditions/LArCaliWave.h index 581e74cbce2863ba23e7135f31f63f49051d652f..c26bf1149c6aa942e09ca4390b3be9e1d1d8931b 100755 --- a/LArCalorimeter/LArRawConditions/LArRawConditions/LArCaliWave.h +++ b/LArCalorimeter/LArRawConditions/LArRawConditions/LArCaliWave.h @@ -61,13 +61,13 @@ public: double dt, int DAC, int isPulsed, - unsigned flag=0); + unsigned flag); LArCaliWave(unsigned nSamples, double dt, int DAC, int isPulsed, - unsigned flag=0); + unsigned flag); LArCaliWave(const std::vector<double>& vAmpl, const std::vector<double>& vErr, @@ -75,7 +75,7 @@ public: double dt, int DAC, int isPulsed, - unsigned flag=0); + unsigned flag); virtual ~LArCaliWave() = default; //@}