21 #include "CoralBase/Blob.h"
27 #include "CaloDetDescr/CaloDetDescrElement.h"
35 #include "tbb/parallel_for.h"
52 return StatusCode::FAILURE;
58 return StatusCode::FAILURE;
68 return StatusCode::FAILURE;
78 msg(MSG::ERROR) <<
"Could not get LArOnlineID helper !" <<
endmsg;
79 return StatusCode::FAILURE;
89 msg(MSG::ERROR) <<
"Could not get LArOnlineID helper !" <<
endmsg;
90 return StatusCode::FAILURE;
109 return StatusCode::SUCCESS;
132 ATH_CHECK(caloSuperCellMgrHandle.isValid());
150 ATH_MSG_INFO(
" Will use helper class for start time." );
175 m_DSPConfig = std::make_unique<LArDSPConfig>(attrList);
181 return StatusCode::FAILURE;
206 std::unique_ptr<tbb::global_control> tbbgc;
209 tbbgc=std::make_unique<tbb::global_control>( tbb::global_control::max_allowed_parallelism,
m_nThreads);
216 tbb::parallel_for(tbb::blocked_range<size_t>(0,
m_allChannelData.size()),looper);
231 std::unique_ptr<LArOFCComplete> larOFCComplete=std::make_unique<LArOFCComplete>();
233 if (
sc.isFailure()) {
234 ATH_MSG_ERROR(
"Failed to set groupingType for LArOFCComplete object" );
238 if (
sc.isFailure()) {
244 std::unique_ptr<LArOFCComplete> larOFCCompleteV2=std::make_unique<LArOFCComplete>();
246 if (
sc.isFailure()) {
247 ATH_MSG_ERROR(
"Failed to set groupingType for LArOFCComplete object" );
251 if (
sc.isFailure()) {
257 std::unique_ptr<LArOFCBinComplete> larOFCBinComplete;
259 larOFCBinComplete=std::make_unique<LArOFCBinComplete>();
261 if (
sc.isFailure()) {
262 ATH_MSG_ERROR(
"Failed to set groupingType for LArOFCBinComplete object" );
266 if (
sc.isFailure() ) {
267 ATH_MSG_ERROR(
"Could not initialize LArOFCComplete data object - exit!" ) ;
273 std::unique_ptr<LArShapeComplete> larShapeComplete;
275 larShapeComplete = std::make_unique<LArShapeComplete>();
277 if (
sc.isFailure()) {
278 ATH_MSG_ERROR(
"Failed to set groupingType for LArShapeComplete object" );
282 if (
sc.isFailure()) {
283 ATH_MSG_ERROR(
"Failed initialize LArShapeComplete object" );
295 if (chanData.faultyOFC) ++nFailed;
298 const int gain=chanData.gain;
300 const std::vector<std::vector<float> > & allPhaseOFCa=chanData.ofc_a;
301 const std::vector<std::vector<float> > & allPhaseOFCb=chanData.ofc_b;
302 const std::vector<std::vector<float> > & allPhaseShape=chanData.shape;
303 const std::vector<std::vector<float> > & allPhaseShapeDer=chanData.shapeDer;
304 larOFCComplete->
set(ch_id,
gain,allPhaseOFCa,allPhaseOFCb,chanData.tstart,chanData.timeBinWidthOFC);
305 if (larOFCBinComplete) larOFCBinComplete->
set(ch_id,
gain,chanData.phasewMaxAt3);
306 if (
m_fillShape ) larShapeComplete->
set(ch_id,
gain,allPhaseShape,allPhaseShapeDer,chanData.tstart,chanData.timeBinWidthOFC);
310 const std::vector<std::vector<float> > & allPhaseOFCV2a=chanData.ofcV2_a;
311 const std::vector<std::vector<float> > & allPhaseOFCV2b=chanData.ofcV2_b;
312 larOFCCompleteV2->
set(ch_id,
gain,allPhaseOFCV2a,allPhaseOFCV2b,chanData.tstart,chanData.timeBinWidthOFC);
319 ATH_MSG_ERROR(
"Number of channels * gains with failed OFC verification: " << nFailed );
327 if (
sc.isFailure()) {
329 return StatusCode::FAILURE;
334 if (
sc.isFailure()) {
335 ATH_MSG_ERROR(
"Could not symlink ILArOFC with LArOFCComplete." );
336 return StatusCode::FAILURE;
344 if (
sc.isFailure()) {
346 return StatusCode::FAILURE;
351 if (
sc.isFailure()) {
352 ATH_MSG_ERROR(
"Could not symlink ILArOFC with LArOFCComplete." );
353 return StatusCode::FAILURE;
358 if (larOFCBinComplete) {
360 if (
sc.isFailure()) {
362 return StatusCode::FAILURE;
370 if (
sc.isFailure()) {
372 return StatusCode::FAILURE;
375 ATH_MSG_DEBUG(
"Trying to symlink ILArShape with LArShapeComplete");
377 if (
sc.isFailure()) {
378 ATH_MSG_ERROR(
"Could not symlink ILArShape with LArShapeComplete." );
379 return StatusCode::FAILURE;
381 ATH_MSG_INFO(
"ILArShape symlink with LArShapeComplete successfully" ) ;
388 ATH_MSG_INFO(
"Reverted corrections of non-cost wave container");
391 return StatusCode::SUCCESS;
414 <<
") is too small to fit your OFC request (" <<
m_nPoints <<
" points)" ) ;
442 aWave = aWave * (1./peak);
462 float maxSampleValAt3=-1;
477 Eigen::MatrixXd acInverseV2;
491 for (
unsigned iPhase=0;iPhase<
m_nPhases;iPhase++) {
494 <<
", Gain = " <<
gain <<
", Phase = " << iPhase <<
":");
498 std::vector<float>& theSamples=chanData.
shape[iPhase];
499 std::vector<float>& theSamplesDer=chanData.
shapeDer[iPhase];
505 for (
unsigned iSample=0;iSample<
m_nSamples;++iSample){
507 theSamples.push_back( aWave.
getSample(tbin) );
508 theSamplesDer.push_back( aDerivedWave.
getSample(tbin) );
513 maxSampleValAt3=theSamples[2];
519 Eigen::VectorXd delta;
521 if (thisChanUseDelta || thisChanUseDeltaV2) {
522 std::vector<float> theSamples32;
523 theSamples32.reserve(32);
524 for (
unsigned iSample=0;iSample<32 ;++iSample){
526 if (tbin>=aWave.
getSize())
continue;
527 theSamples32.push_back( aWave.
getSample(tbin) );
534 std::vector<float>& vOFC_a= chanData.
ofc_a[iPhase];
535 std::vector<float>& vOFC_b= chanData.
ofc_b[iPhase];
538 optFiltPed(theSamples,theSamplesDer,acInverse,vOFC_a,vOFC_b);
540 if (thisChanUseDelta) {
541 optFiltDelta(theSamples,theSamplesDer,acInverse,delta,vOFC_a,vOFC_b);
544 optFilt(theSamples,theSamplesDer,acInverse,vOFC_a,vOFC_b);
557 std::vector<float>& vOFCV2_a= chanData.
ofcV2_a[iPhase];
558 std::vector<float>& vOFCV2_b= chanData.
ofcV2_b[iPhase];
561 optFiltPed(theSamples,theSamplesDer,acInverseV2,vOFCV2_a,vOFCV2_b);
563 if (thisChanUseDeltaV2) {
564 optFiltDelta(theSamples,theSamplesDer,acInverseV2,delta,vOFCV2_a,vOFCV2_b);
567 optFilt(theSamples,theSamplesDer,acInverseV2,vOFCV2_a,vOFCV2_b);
598 if (
sc.isFailure()) {
607 WAVEIT it_e=waveCnt->
end(
gain);
608 for (;
it!=it_e;++
it) {
610 if (
cabling->isOnlineConnected (chid)){
619 return StatusCode::SUCCESS;
636 ATH_MSG_INFO(
"LArCaliWaveContainer: Corrections not yet applied, applying them now..." );
638 ATH_MSG_ERROR(
"Failed to apply corrections to LArCaliWaveContainer!" );
639 return StatusCode::FAILURE;
642 ATH_MSG_INFO(
"Applied corrections to non-const Wave container");
650 return StatusCode::FAILURE;
659 WAVEIT it_e=waveCnt->
end(
gain);
660 for (;
it!=it_e;++
it) {
662 for (
const auto& cw : wVec) {
671 return StatusCode::SUCCESS;
675 void LArOFCAlg::optFilt(
const std::vector<float> &gWave,
const std::vector<float> &gDerivWave,
const Eigen::MatrixXd& acInverse,
676 std::vector<float>& vecOFCa, std::vector<float>& vecOFCb) {
677 assert(gWave.size()==gDerivWave.size());
679 const int optNpt = gWave.size();
681 Eigen::VectorXd gResp(optNpt), gDerivResp(optNpt);
682 for (
int i=0;
i<optNpt;
i++) {
684 gDerivResp[
i] = gDerivWave[
i];
687 Eigen::Matrix2d isol;
689 (gResp.transpose()*acInverse*gResp)[0],
690 (gResp.transpose()*acInverse*gDerivResp)[0],
691 (gDerivResp.transpose()*acInverse*gResp)[0],
692 (gDerivResp.transpose()*acInverse*gDerivResp)[0];
695 Eigen::Vector2d Atau;
696 Eigen::Vector2d Ktemp;
697 Eigen::Matrix2d isolInv = isol.inverse();
706 Atau = isolInv*Ktemp;
709 Eigen::VectorXd OFCa = Amp[0]*acInverse*gResp + Amp[1]*acInverse*gDerivResp;
710 Eigen::VectorXd OFCb = Atau[0]*acInverse*gResp + Atau[1]*acInverse*gDerivResp;
713 vecOFCa.resize(optNpt);
714 vecOFCb.resize(optNpt);
715 for (
int i=0;
i<optNpt;
i++) {
721 void LArOFCAlg::optFiltPed(
const std::vector<float> &gWave,
const std::vector<float> &gDerivWave,
const Eigen::MatrixXd& acInverse,
722 std::vector<float>& vecOFCa, std::vector<float>& vecOFCb) {
723 assert(gWave.size()==gDerivWave.size());
725 const int optNpt = gWave.size();
727 Eigen::VectorXd gResp(optNpt), gDerivResp(optNpt);
728 for (
int i=0;
i<optNpt;
i++) {
730 gDerivResp[
i] = gDerivWave[
i];
733 Eigen::Matrix3d isol;
734 Eigen::Vector3d Kunit(1.,1.,1.);
735 auto s3=(gDerivResp.transpose()*acInverse*gResp)[0];
736 auto s4=(gResp.transpose()*acInverse*Kunit)[0];
737 auto s5=(gDerivResp.transpose()*acInverse*Kunit)[0];
739 (gResp.transpose()*acInverse*gResp)[0],
s3,
s4,
740 s3, (gDerivResp.transpose()*acInverse*gResp)[0], s5,
741 s4, s5, (Kunit.transpose()*acInverse*Kunit)[0];
744 Eigen::Vector3d Atau;
745 Eigen::Vector3d Ktemp;
746 Eigen::Matrix3d isolInv = isol.inverse();
757 Atau = isolInv*Ktemp;
760 Eigen::VectorXd OFCa = Amp[0]*acInverse*gResp + Amp[1]*acInverse*gDerivResp + Amp[2]*acInverse*Kunit;
761 Eigen::VectorXd OFCb = Atau[0]*acInverse*gResp + Atau[1]*acInverse*gDerivResp + Atau[2]*acInverse*Kunit;
764 vecOFCa.resize(optNpt);
765 vecOFCb.resize(optNpt);
766 for (
int i=0;
i<optNpt;
i++) {
774 const Eigen::MatrixXd& acInverse,
const Eigen::VectorXd& delta,
775 std::vector<float>& vecOFCa, std::vector<float>& vecOFCb) {
779 assert(gWave.size()==gDerivWave.size());
781 const int optNpt = gWave.size();
783 Eigen::VectorXd gResp(optNpt), gDerivResp(optNpt);
784 for (
int i=0;
i<optNpt;
i++) {
786 gDerivResp[
i] = gDerivWave[
i];
801 Eigen::Matrix3d isol;
803 (gResp.transpose()*acInverse*gResp)[0],
804 (gResp.transpose()*acInverse*gDerivResp),
805 (gResp.transpose()*acInverse*delta)[0],
807 (gDerivResp.transpose()*acInverse*gResp)[0],
808 (gDerivResp.transpose()*acInverse*gDerivResp)[0],
809 (gDerivResp.transpose()*acInverse*delta)[0],
811 (delta.transpose()*acInverse*gResp)[0],
812 (delta.transpose()*acInverse*gDerivResp)[0],
813 (delta.transpose()*acInverse*delta)[0];
817 Eigen::Vector3d Atau;
818 Eigen::Vector3d Ktemp;
819 Eigen::Matrix3d isolInv = isol.inverse();
831 Atau = isolInv*Ktemp;
834 Eigen::VectorXd OFCa = Amp[0]*acInverse*gResp + Amp[1]*acInverse*gDerivResp + Amp[2]*acInverse * delta;
835 Eigen::VectorXd OFCb = Atau[0]*acInverse*gResp + Atau[1]*acInverse*gDerivResp + Atau[2]*acInverse * delta ;
839 vecOFCa.resize(optNpt);
840 vecOFCb.resize(optNpt);
841 for (
int i=0;
i<optNpt;
i++) {
849 const float LArOFCAlg::m_fcal3Delta[5] ={0.0790199937765, 0.0952000226825, 0.0790199937765, 0.0952000226825, 0.0790199937765};
850 const float LArOFCAlg::m_fcal2Delta[5]={-0.01589001104, -0.0740399733186, -0.01589001104, -0.0740399733186, -0.01589001104};
851 const float LArOFCAlg::m_fcal1Delta[5] ={0.0679600232979, -0.139479996869, 0.0679600232979, -0.139479996869, 0.0679600232979};
881 for (
unsigned int i = 0;
i<samples.size();++
i) {
920 if (
cabling->isOnlineConnected (chid)){
924 ATH_MSG_ERROR(
" dde = 0 , onl_id, ofl_id= "<< chid<<
" "<<ofl_id );
933 if (fabs(dde->
eta())>4.0){
943 if (fabs(dde->
eta())>4.0){
959 const std::vector<float>& Shape,
const char* ofcversion,
const unsigned phase)
const {
962 float recAmpl=0, recTime=0;
963 for (
unsigned iSample=0;iSample<
m_nSamples;++iSample){
965 #ifdef LAROFCALG_DEBUGOUTPUT
966 ATH_MSG_VERBOSE(
"a["<<iSample<<
"]="<<vOFC_a[iSample] <<
" b["<<iSample<<
"]="<<vOFC_b[iSample]
967 <<
" Sample=" << aWave.getSample(tbin));
969 recAmpl += OFCa[iSample] * Shape[iSample];
970 recTime += OFCb[iSample] * Shape[iSample];
978 ATH_MSG_WARNING(
"Applying phase " <<
phase <<
" of " << ofcversion <<
" to original wave yields an Amplitude of "<< recAmpl
984 ATH_MSG_WARNING(
"Applying phase " <<
phase <<
" of " << ofcversion <<
" to original wave yields a time offset of " << recTime
993 mLog << MSG::WARNING <<
"OFCs";
1000 #ifdef LAROFCALG_DEBUGOUTPUT
1001 #undef LAROFCALG_DEBUGOUTPUT