Commit 46da2ed3 authored by Luigi Del Buono's avatar Luigi Del Buono
Browse files

Bug fixes to inmprove output similarity with normal clustering code

parent 920e7eef
......@@ -15,17 +15,17 @@ DECLARE_ALGORITHM_FACTORY( FTFPGAClusterCreator )
FTFPGAClusterCreator::FTFPGAClusterCreator(const std::string& name, ISvcLocator* pSvcLocator) :
MultiTransformer(name, pSvcLocator,
{ KeyValue{"InputLocation" , LHCb::MCFTDigitLocation::Default} },
{ KeyValue{"OutputLocation" , LHCb::FTClusterLocation::Default },
KeyValue{"MCToClusterLocation",
Links::location(LHCb::FTLiteClusterLocation::Default )},
KeyValue{"MCToClusterExtLocation",
Links::location(LHCb::FTLiteClusterLocation::Default+"WithSpillover" )},
KeyValue{"HitToClusterLocation",
Links::location(LHCb::FTLiteClusterLocation::Default + "2MCHits" )},
KeyValue{"HitToClusterExtLocation",
Links::location(LHCb::FTLiteClusterLocation::Default + "2MCHitsWithSpillover")}
} ) {}
{ KeyValue{"InputLocation" , LHCb::MCFTDigitLocation::Default} },
{ KeyValue{"OutputLocation" , LHCb::FTClusterLocation::Default },
KeyValue{"MCToClusterLocation",
Links::location(LHCb::FTLiteClusterLocation::Default )},
KeyValue{"MCToClusterExtLocation",
Links::location(LHCb::FTLiteClusterLocation::Default+"WithSpillover" )},
KeyValue{"HitToClusterLocation",
Links::location(LHCb::FTLiteClusterLocation::Default + "2MCHits" )},
KeyValue{"HitToClusterExtLocation",
Links::location(LHCb::FTLiteClusterLocation::Default + "2MCHitsWithSpillover")}
} ) {}
//=============================================================================
// Initialization
......@@ -66,12 +66,12 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
//(because of channel numbering and bit representation in C++ vs FPGA code!!
//an easy adaptation shall be made in the FPGA)
/*
static std::array<unsigned int, 16> SIZE_TABLE4 = { //original FPGA
static std::array<unsigned int, 16> SIZE_TABLE4 = { //original FPGA
0, 0, 0, 0,
0, 0, 0, 0,
1, 1, 1, 1,
2, 2, 3, 4
};
};
*/
static std::array<unsigned int, 16> SIZE_TABLE4 = {
0, 1, 0, 2,
......@@ -81,7 +81,7 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
};
/*
static std::array<unsigned int, 32> SIZE_TABLE5 = { //original FPGA
static std::array<unsigned int, 32> SIZE_TABLE5 = { //original FPGA
0, 0, 0, 0,
0, 0, 0, 0,
0, 0, 0, 0,
......@@ -90,7 +90,7 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
1, 1, 1, 1,
2, 2, 2, 2,
3, 3, 4, 5
};
};
*/
static std::array<unsigned int, 32> SIZE_TABLE5 = {
0, 1, 0, 2,
......@@ -108,7 +108,6 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
// Digits loop : Digit Container is sorted wrt channelID
auto digitIter = digits.begin();
unsigned int SipmAdcValues[nChannelsInSipm] = {};
std::bitset<nChannelsInSipm> SipmBinValues = {};
LHCb::MCFTDigit* SipmDigits[nChannelsInSipm] = {};
unsigned int lastUniqueSiPM((*digitIter)->channelID().uniqueSiPM());
......@@ -122,10 +121,9 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
if (digit->channelID().uniqueSiPM() == lastUniqueSiPM) {
unsigned int nchan = digit->channelID().channel();
SipmDigits[nchan] = digit;
SipmAdcValues[nchan] = digit->adcCount();
if( (!m_usePEnotADC && digit->adcCount() >= m_adcThreshold1) ||
(m_usePEnotADC && digit->photoElectrons() >= m_adcThreshold1) ) {
SipmBinValues.set(nchan);
(m_usePEnotADC && digit->photoElectrons() >= m_adcThreshold1) ) {
SipmBinValues.set(nchan);
}
digitIter++;
}
......@@ -147,7 +145,7 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
//half-sipm actions loop
unsigned int SipmAdcSizes[nChannelsInSipm] = {};
unsigned int SipmAdcSums[nChannelsInSipm] = {};
float SipmAdcSums[nChannelsInSipm] = {};
int SipmAdcSizes2[nChannelsInSipm] = {};
std::bitset<nChannelsInSipm> SipmBigFlags = {};
for ( int ichunk = 0; ichunk < nchunk; ichunk++) {
......@@ -161,20 +159,26 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
if (int_clusterMaxWidth == 4) SipmAdcSizes[ichan] = SIZE_TABLE4[locChunk];
else if (int_clusterMaxWidth == 5) SipmAdcSizes[ichan] = SIZE_TABLE5[locChunk];
else {
int curSize(0);
int curSize(0);
while( curSize < int_clusterMaxWidth && locbits.test(curSize) ) ++curSize;
SipmAdcSizes[ichan] = curSize;
SipmAdcSizes[ichan] = curSize;
};
ub.SipmBinHalfs[ichunk].SipmBinHalf >>= 1; //shift towards smaller channel number wrto SipmAdcValues array
ub.SipmBinHalfs[ichunk].SipmBinHalf >>= 1; //shift towards smaller channel number wrto Sipm arrays
}
//b) build amplitude sum array (in 2 chunks again)
for ( int ichan = ichunk*chunkSize; ichan < (ichunk+1)*chunkSize; ichan++ ) {
int ilocchan(0), locsum(0);
int ilocchan(0);
float locsum(0);
while( (ichan + ilocchan) < (ichunk+1)*chunkSize
&& ilocchan < int_clusterMaxWidth
&& SipmBinValues.test(ichan + ilocchan) ) {
locsum+=SipmAdcValues[ichan + ilocchan];
&& ilocchan < int_clusterMaxWidth
&& SipmBinValues.test(ichan + ilocchan) ) {
float channelWeight = ( m_usePEnotADC ? SipmDigits[ichan + ilocchan] -> photoElectrons() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold3 ? m_adcThreshold3Weight.value() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold2 ? m_adcThreshold2Weight.value() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold1 ? m_adcThreshold1Weight.value() : 0. );
locsum+=channelWeight;
ilocchan++;
};
SipmAdcSums[ichan] = locsum;
......@@ -193,13 +197,15 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
SipmBigFlags.set(ichan);
}
else if (m_keepEdgeCluster && (ichan == ichunk*chunkSize)) {
SipmAdcSizes2[ichan] = SipmAdcSizes[ichan]; //keep right edge clusters
SipmAdcSizes2[ichan] = SipmAdcSizes[ichan]; //keep right edge clusters (no conditions)
}
else if ( m_keepEdgeCluster && ((ichan + SipmAdcSizes[ichan] - (ichunk+1)*chunkSize == 0) && !SipmBinValues.test(ichan - 1)) ) {
SipmAdcSizes2[ichan] = SipmAdcSizes[ichan]; //keep left edge clusters
SipmAdcSizes2[ichan] = SipmAdcSizes[ichan]; //keep left edge clusters (no conditions)
}
else if ( (SipmAdcSizes[ichan] > m_clusterMinWidth) //inequalities changed : m_clusterMinWidth=2 in FPGA logic !!! (not 1)
&& (SipmAdcSums[ichan] < m_adcThreshold1Weight.value() + m_adcThreshold2Weight.value() - 0.5) ) {
else if ( ((ichan+int_clusterMaxWidth < (ichunk+1)*chunkSize)
&& SipmBigFlags.test(ichan+int_clusterMaxWidth)) //don't filter (1st) subcluster of large cluster (like the others)
|| ((SipmAdcSizes[ichan] > m_clusterMinWidth) //inequalities : m_clusterMinWidth=2 in FPGA logic ! (not 1)
&& (SipmAdcSums[ichan] < m_adcThreshold1Weight.value() + m_adcThreshold2Weight.value() - 0.5)) ) {
SipmAdcSizes2[ichan] = 0;
}
else if ( (SipmAdcSizes[ichan] <= m_clusterMinWidth) //inequalities changed (see above)
......@@ -223,7 +229,7 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
};
//now "SipmAdcSizes2" is an nchannels array containing zeros except when a cluster exists
//this array is parallel to "SipmDigits" , "SipmAdcValues" and "SipmBigFlags", we can go on with the (old)
//this array is parallel to "SipmDigits" and "SipmBigFlags", we can go on with the (old)
//remaining cluster infos computations
......@@ -244,69 +250,69 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
// (threshold1 < threshold2 < threshold3)
for( int ilocchan = 0; ilocchan < SipmAdcSizes2[ichan]; ilocchan++) {
float channelWeight = ( m_usePEnotADC ? SipmDigits[ichan + ilocchan] -> photoElectrons() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold3 ? m_adcThreshold3Weight.value() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold2 ? m_adcThreshold2Weight.value() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold1 ? m_adcThreshold1Weight.value() : 0. );
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold3 ? m_adcThreshold3Weight.value() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold2 ? m_adcThreshold2Weight.value() :
SipmDigits[ichan + ilocchan]->adcCount() >= m_adcThreshold1 ? m_adcThreshold1Weight.value() : 0. );
totalCharge += channelWeight;
totalChargePE += SipmDigits[ichan + ilocchan] -> photoElectrons();
// mean position will be [ (rel. pos. from left) * charge ] / totalCharge (see below)
ids.push_back(SipmDigits[ichan + ilocchan]->channelID());
ids.push_back(SipmDigits[ichan + ilocchan]->channelID());
wsumPosition += ilocchan * channelWeight;
} //end of loop over digits in cluster
// compute position : when cluster is an edge of a double or big cluster keep
// the middle of the cluster only (no weighting), otherwise use weights
// also add channelID (uint) offset
double clusPosition;
if(SipmBigFlags.test(ichan)) { //big cluster trail
clusPosition = SipmDigits[ichan]->channelID() + float(SipmAdcSizes2[ichan]-1)/2;
}
// the middle of the cluster only (no weighting), otherwise use weights
// also add channelID (uint) offset
double clusPosition;
if(SipmBigFlags.test(ichan)) { //big cluster trail
clusPosition = SipmDigits[ichan]->channelID() + float(SipmAdcSizes2[ichan]-1)/2;
}
else if( (ichan+int_clusterMaxWidth<nChannelsInSipm) && SipmBigFlags.test(ichan+int_clusterMaxWidth) ) { //big cluster start
clusPosition = SipmDigits[ichan]->channelID() + float(SipmAdcSizes2[ichan]-1)/2;
}
else {
clusPosition = SipmDigits[ichan]->channelID() + wsumPosition/totalCharge;
clusPosition = SipmDigits[ichan]->channelID() + float(SipmAdcSizes2[ichan]-1)/2;
}
// The fractional position is defined in (-0.250, 0.750)
unsigned int clusChanPosition = std::floor(clusPosition-m_lowestFraction);
double fractionChanPosition = (clusPosition-clusChanPosition);
int frac = int(2*(fractionChanPosition-m_lowestFraction));
else {
clusPosition = SipmDigits[ichan]->channelID() + wsumPosition/totalCharge;
}
// The fractional position is defined in (-0.250, 0.750)
unsigned int clusChanPosition = std::floor(clusPosition-m_lowestFraction);
double fractionChanPosition = (clusPosition-clusChanPosition);
int frac = int(2*(fractionChanPosition-m_lowestFraction));
//Define new lite cluster and add to container
LHCb::FTCluster *newCluster = new LHCb::FTCluster(clusChanPosition, fractionChanPosition, frac,
ids, SipmBigFlags.test(ichan),
m_storePECharge ? totalChargePE : totalCharge);
clusterCont.insert(newCluster);
//--Setup MCParticle to cluster link
// (also for spillover / not accepted clusters, as long as it's not noise)
// check contributions from MCHits, for linking
for( int ilocchan = 0; ilocchan < SipmAdcSizes2[ichan]; ilocchan++) {
for (const auto& deposit : SipmDigits[ichan + ilocchan]->deposits() ) {
const auto& mcHit = deposit->mcHit();
if( mcHit != nullptr ) {
contributingMCHits.insert(mcHit);
LHCb::FTCluster *newCluster = new LHCb::FTCluster(clusChanPosition, fractionChanPosition, frac,
SipmAdcSizes2[ichan], ids, 1 ,//SipmBigFlags.test(ichan),
m_storePECharge ? totalChargePE : totalCharge);
clusterCont.insert(newCluster);
//--Setup MCParticle to cluster link
// (also for spillover / not accepted clusters, as long as it's not noise)
// check contributions from MCHits, for linking
for( int ilocchan = 0; ilocchan < SipmAdcSizes2[ichan]; ilocchan++) {
for (const auto& deposit : SipmDigits[ichan + ilocchan]->deposits() ) {
const auto& mcHit = deposit->mcHit();
if( mcHit != nullptr ) {
contributingMCHits.insert(mcHit);
if( mcHit->mcParticle() != nullptr )
contributingMCParticles.insert(mcHit->mcParticle());
}
}
} //end of loop over digits in cluster
for(const auto& i : contributingMCParticles) {
mcToClusterLinkExtended.link(clusChanPosition, i ) ; // Needed for xdigi / xdst
if ( i->parent()->registry()->identifier() ==
"/Event/"+LHCb::MCParticleLocation::Default)
mcToClusterLink.link(clusChanPosition, i ) ;
}
// Setup MCHit to cluster link
for(const auto& i : contributingMCHits ) {
hitToClusterLinkExtended.link(clusChanPosition, i ) ; // Needed for xdigi / xdst
if ( i->parent()->registry()->identifier() == "/Event/"+LHCb::MCHitLocation::FT)
hitToClusterLink.link(clusChanPosition, i ) ;
}
for(const auto& i : contributingMCParticles) {
mcToClusterLinkExtended.link(clusChanPosition, i ) ; // Needed for xdigi / xdst
if ( i->parent()->registry()->identifier() ==
"/Event/"+LHCb::MCParticleLocation::Default)
mcToClusterLink.link(clusChanPosition, i ) ;
}
// Setup MCHit to cluster link
for(const auto& i : contributingMCHits ) {
hitToClusterLinkExtended.link(clusChanPosition, i ) ; // Needed for xdigi / xdst
if ( i->parent()->registry()->identifier() == "/Event/"+LHCb::MCHitLocation::FT)
hitToClusterLink.link(clusChanPosition, i ) ;
}
} //end of existing cluster
......@@ -316,7 +322,6 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
//----------------------> refresh for new sipm
lastUniqueSiPM = digit->channelID().uniqueSiPM();
for( int i = 0; i < nChannelsInSipm; i++) {
SipmAdcValues[i] = 0;
SipmDigits[i] = 0;
};
SipmBinValues.reset();
......@@ -328,16 +333,16 @@ FTFPGAClusterCreator::operator()(const LHCb::MCFTDigits& digits) const {
if ( msgLevel( MSG::DEBUG) )
debug() << "Number of FTClusters created: " << clusterCont.size() << endmsg;
// Extra checks
// Extra checks
static_assert(std::is_move_constructible<LHCb::FTClusters>::value,
"FTClusters must be move constructible for this to work.");
return std::make_tuple( std::move(clusterCont),
std::move(mcToClusterLink.links()),
std::move(mcToClusterLinkExtended.links()),
std::move(hitToClusterLink.links()),
std::move(hitToClusterLinkExtended.links()));
std::move(mcToClusterLink.links()),
std::move(mcToClusterLinkExtended.links()),
std::move(hitToClusterLink.links()),
std::move(hitToClusterLinkExtended.links()));
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment