/************************************************************************** * This file is property of and copyright by the ALICE HLT Project * * All rights reserved. * * * * Primary Authors: * * Indranil Das * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id$ */ /////////////////////////////////////////////// //Author : Indranil Das, SINP, INDIA // Sukalyan Chattopadhyay, SINP, INDIA // //Email : indra.das@saha.ac.in // sukalyan.chattopadhyay@saha.ac.in // // This class implements a hit reconstruction algorithm for the dimuon // high level trigger. // The algorithm finds 3 pad clusters by looking for unique pads with a charge // above a certain threshold. A centre of gravity type calculation is applied // to the three pads forming the cluster to find the hit's X or Y coordinate // along the non-bending and bending planes individually. // The sepperate X and Y coordinates are then merged to give the full coordinate // of the hit. ///////////////////////////////////////////////// #include "AliHLTMUONHitReconstructor.h" #include "AliHLTMUONRecHitsBlockStruct.h" #include #include const AliHLTInt32_t AliHLTMUONHitReconstructor::fgkDetectorId = 0xA00; const AliHLTInt32_t AliHLTMUONHitReconstructor::fgkDDLOffSet = 12; const AliHLTInt32_t AliHLTMUONHitReconstructor::fgkNofDDL = 8; const AliHLTInt32_t AliHLTMUONHitReconstructor::fgkDDLHeaderSize = 8; const AliHLTInt32_t AliHLTMUONHitReconstructor::fgkLutLine = 59648 + 1; AliHLTMUONHitReconstructor::AliHLTMUONHitReconstructor() : AliHLTLogging(), fHLTMUONDecoder(), fkBlockHeaderSize(8), fkDspHeaderSize(8), fkBuspatchHeaderSize(4), fDCCut(0), fPadData(NULL), fLookUpTableData(NULL), fRecPoints(NULL), fRecPointsCount(NULL), fMaxRecPointsCount(0), fCentralCountB(0), fCentralCountNB(0), fDigitPerDDL(0), fCentralChargeB(NULL), fCentralChargeNB(NULL), fRecX(NULL), fRecY(NULL), fAvgChargeX(NULL), fAvgChargeY(NULL), fNofBChannel(NULL), fNofNBChannel(NULL), fNofFiredDetElem(0), fIdToEntry() { /// Default constructor fkBlockHeaderSize = 8; fkDspHeaderSize = 8; fkBuspatchHeaderSize = 4; try { fPadData = new AliHLTMUONPad[fgkLutLine]; } catch (const std::bad_alloc&) { HLTError("Dynamic memory allocation failed for fPadData in constructor."); throw; } fPadData[0].fDetElemId = 0; fPadData[0].fIX = 0 ; fPadData[0].fIY = 0 ; fPadData[0].fRealX = 0.0 ; fPadData[0].fRealY = 0.0 ; fPadData[0].fRealZ = 0.0 ; fPadData[0].fHalfPadSize = 0.0 ; fPadData[0].fPlane = -1 ; fPadData[0].fCharge = 0 ; bzero(fGetIdTotalData, 336*237*2*sizeof(int)); } AliHLTMUONHitReconstructor::~AliHLTMUONHitReconstructor() { /// Default destructor if (fPadData) { delete [] fPadData; fPadData = NULL; } } void AliHLTMUONHitReconstructor::SetLookUpTable( const AliHLTMUONHitRecoLutRow* lookupTable, const IdManuChannelToEntry* idToEntry ) { /// Sets the Lookup table (LUT) containing the position of each pad with /// electronic channel associated with it. Also the appropriate manu /// channel ID mapping to LUT row is also set. assert( lookupTable != NULL ); assert( idToEntry != NULL ); fIdToEntry = idToEntry; fLookUpTableData = lookupTable; } bool AliHLTMUONHitReconstructor::Run( const AliHLTUInt32_t* rawData, AliHLTUInt32_t rawDataSize, AliHLTMUONRecHitStruct* recHit, AliHLTUInt32_t& nofHit ) { // main function called by HLTReconstructor to perform DHLT Hitreconstruction fRecPoints = recHit; fMaxRecPointsCount = nofHit; fRecPointsCount = &nofHit; *fRecPointsCount = 0; fDigitPerDDL = 0; if (not DecodeDDL(rawData, rawDataSize)) { // Dont need to log any message again. Already done so in DecodeDDL. return false; } if (fDigitPerDDL == 1) { // There are no digits to process so stop here. return true; } if (not FindRecHits()) { HLTError("Failed to generate RecHits"); return false; } return true; } bool AliHLTMUONHitReconstructor::DecodeDDL(const AliHLTUInt32_t* rawData,AliHLTUInt32_t rawDataSize) { //function to decode Raw Data AliHLTMUONRawDecoder& handler = reinterpret_cast(fHLTMUONDecoder.GetHandler()); UInt_t bufferSize = UInt_t(rawDataSize*sizeof(AliHLTUInt32_t)); handler.SetDCCut(fDCCut); handler.SetPadData(fPadData); handler.SetLookUpTable(fLookUpTableData); handler.SetIdManuChannelToEntry(fIdToEntry); handler.SetNofFiredDetElemId(fNofFiredDetElem); handler.SetMaxFiredPerDetElem(fMaxFiredPerDetElem); if(!fHLTMUONDecoder.Decode(rawData,bufferSize)) return false; fDigitPerDDL = handler.GetDataCount(); fMaxFiredPerDetElem[fNofFiredDetElem-1] = handler.GetDataCount(); if(fDigitPerDDL == 1){ HLTInfo("An Empty DDL File found"); } return true; } bool AliHLTMUONHitReconstructor::FindRecHits() { // fuction that calls hit reconstruction detector element-wise. assert( fCentralChargeB == NULL ); assert( fCentralChargeNB == NULL ); assert( fRecX == NULL ); assert( fRecY == NULL ); assert( fAvgChargeX == NULL ); assert( fAvgChargeY == NULL ); assert( fNofBChannel == NULL ); assert( fNofNBChannel == NULL ); bool resultOk = false; for(int iDet=0; iDet0) FindCentralHits(fMaxFiredPerDetElem[iDet-1],fMaxFiredPerDetElem[iDet]); else // minimum value is 1 because dataCount in ReadDDL starts from 1 instead of 0; FindCentralHits(1,fMaxFiredPerDetElem[iDet]); try { fRecY = new float[fCentralCountB]; HLTDebug("Allocated fRecY with %d elements.", fCentralCountB); fRecX = new float[fCentralCountNB]; HLTDebug("Allocated fRecX with %d elements.", fCentralCountNB); fAvgChargeY = new float[fCentralCountB]; HLTDebug("Allocated fAvgChargeY with %d elements.", fCentralCountB); fAvgChargeX = new float[fCentralCountNB]; HLTDebug("Allocated fAvgChargeX with %d elements.", fCentralCountNB); fNofBChannel = new int[fCentralCountB]; HLTDebug("Allocated fNofBChannel with %d elements.", fCentralCountB); fNofNBChannel = new int[fCentralCountNB]; HLTDebug("Allocated fNofNBChannel with %d elements.", fCentralCountNB); resultOk = true; } catch(const std::bad_alloc&){ HLTError("Dynamic memory allocation failed for internal arrays."); resultOk = false; } } if (resultOk) RecXRecY(); if (resultOk) { resultOk = MergeRecHits(); } if (resultOk) { if(iDet==0) // minimum value in loop is 1 because dataCount in ReadDDL starts from 1 instead of 0; for(int i=1;i (0+1) && fPadData[iPad].fIY < (79 - 1)){ //if(fPadData[iPad].fIY > 0){ if(fCentralCountB>0){ hasFind = false; for(b = 0;b fPadData[idManuChannelCentral].fCharge){ fCentralChargeB[b] = iPad; }// if condn on pad charge }// if condon on pad position }// for loop over b if(!hasFind){ fCentralChargeB[fCentralCountB] = iPad; fCentralCountB++; } } else{ fCentralChargeB[fCentralCountB] = iPad; fCentralCountB++; }// check the size of centralHitB for(b = 0;b 2 (to avoid edge value pb) }// B/Nb checking else{ if(fCentralCountNB>0){ hasFind = false; for(nb = 0;nb fPadData[idManuChannelCentral].fCharge){ fCentralChargeNB[nb] = iPad; }// if condn over to find higher charge }// if condn over to find position }// for loop over presently all nb values if(!hasFind){ fCentralChargeNB[fCentralCountNB] = iPad; fCentralCountNB++; } }// centralHitNB size test else{ fCentralChargeNB[fCentralCountNB] = iPad; fCentralCountNB++; }// centralHitNB size test }// fill for bending and nonbending hit }// detElemId loop } void AliHLTMUONHitReconstructor::RecXRecY() { // find reconstructed X and Y for each plane separately assert( fRecX != NULL ); assert( fRecY != NULL ); assert( fAvgChargeX != NULL ); assert( fAvgChargeY != NULL ); assert( fNofBChannel != NULL ); assert( fNofNBChannel != NULL ); int b,nb; int idCentral; int idLower = 0; int idUpper = 0; int idRight = 0; int idLeft = 0; for(b=0;b0) fNofBChannel[b]++ ; if(fPadData[idCentral].fCharge>0) fNofBChannel[b]++ ; if(fPadData[idUpper].fCharge>0) fNofBChannel[b]++ ; HLTDebug("lower : %d, middle : %d, upper : %d, nofChannel : %d",fPadData[idLower].fCharge, fPadData[idCentral].fCharge,fPadData[idUpper].fCharge,fNofBChannel[b]); HLTDebug("RecY[%d] : %f",b,fRecY[b]); } for(nb=0;nb0) fNofNBChannel[nb]++ ; if(fPadData[idCentral].fCharge>0) fNofNBChannel[nb]++ ; if(fPadData[idRight].fCharge>0) fNofNBChannel[nb]++ ; HLTDebug("left : %d, middle : %d, right : %d, nofChannel : %d",fPadData[idLeft].fCharge, fPadData[idCentral].fCharge,fPadData[idRight].fCharge,fNofNBChannel[nb]); HLTDebug("RecX[%d] : %f",nb,fRecX[nb]); } } bool AliHLTMUONHitReconstructor::MergeRecHits() { // Merge reconstructed hits first over same plane then bending plane with non-bending plane assert( fRecX != NULL ); assert( fRecY != NULL ); assert( fAvgChargeX != NULL ); assert( fAvgChargeY != NULL ); assert( fNofBChannel != NULL ); assert( fNofNBChannel != NULL ); int idCentralB,idCentralNB ; float padCenterXB; float padCenterYNB; float diffX,diffY; float halfPadLengthX,halfPadLengthY; // MERGE Bending Plane hits, which are placed side by side for(int i=0;i fAvgChargeY[j]){ fRecY[i] = (fRecY[i]*fAvgChargeY[i] + fRecY[j]*fAvgChargeY[j] )/(fAvgChargeY[i] + fAvgChargeY[j]); fRecY[j] = 0.0; } else{ fRecY[j] = (fRecY[i]*fAvgChargeY[i] + fRecY[j]*fAvgChargeY[j] )/(fAvgChargeY[i] + fAvgChargeY[j]); fRecY[i] = 0.0; }// search for higher charge }//pad position }//j for loop }//if fRecY[i] != 0.0 }// i for loop // MERGE Non Bending Plane hits, which are place side by side for(int i=0;i fAvgChargeX[j]){ fRecX[i] = (fRecX[i]*fAvgChargeX[i] + fRecX[j]*fAvgChargeX[j] )/(fAvgChargeX[i] + fAvgChargeX[j]); fRecX[j] = 0.0; } else{ fRecX[j] = (fRecX[i]*fAvgChargeX[i] + fRecX[j]*fAvgChargeX[j] )/(fAvgChargeX[i] + fAvgChargeX[j]); fRecX[i] = 0.0; }// search for higher charge }//pad position }//j for loop }//if fRecX[i] != 0.0 }// i for loop // Merge bending Plane hits with Non Bending for(int b=0;b fabsf(padCenterXB)) diffX = fabsf(fRecX[nb]) - fabsf(padCenterXB); else diffX = fabsf(padCenterXB) - fabsf(fRecX[nb]); if(fabsf(padCenterYNB)>fabsf(fRecY[b])) diffY = fabsf(padCenterYNB) - fabsf(fRecY[b]); else diffY = fabsf(fRecY[b]) - fabsf(padCenterYNB); if(diffX < halfPadLengthX && diffY < halfPadLengthY ){//&& fPadData[idCentralB].fIY != 0){ if(fNofBChannel[b]==3) fRecY[b] += 0.025*sin(12.0*(fRecY[b] - fPadData[idCentralB].fRealY)) ; fRecX[nb] += 0.075*sin(9.5*(fRecX[nb] - fPadData[idCentralNB].fRealX)) ; // First check that we have not overflowed the buffer. if((*fRecPointsCount) == fMaxRecPointsCount){ HLTError("Number of RecHit (i.e. %d) exceeds the max number of RecHit limit %d.",(*fRecPointsCount),fMaxRecPointsCount); return false; } //fRecPoints[(*fRecPointsCount)].fId = idCentralB; fRecPoints[(*fRecPointsCount)].fX = fRecX[nb]; fRecPoints[(*fRecPointsCount)].fY = fRecY[b]; fRecPoints[(*fRecPointsCount)].fZ = fPadData[idCentralB].fRealZ; // fRecPoints[(*fRecPointsCount)].fXCenter = fPadData[idCentralNB].fRealX; // fRecPoints[(*fRecPointsCount)].fYCenter = fPadData[idCentralB].fRealY; // fRecPoints[(*fRecPointsCount)].fNofBChannel = fNofBChannel[b]; // fRecPoints[(*fRecPointsCount)].fNofNBChannel = fNofNBChannel[nb]; // fRecPoints[(*fRecPointsCount)].fDetElemId = (AliHLTUInt32_t)fPadData[idCentralB].fDetElemId; (*fRecPointsCount)++; }//if lies wihtin 5.0 mm }// condn over fRecX ! = 0.0 }// loop over NB side }// condn on fRecY[b] != 0.0 }// loop over B side; return true; } void AliHLTMUONHitReconstructor::Clear() { // function to clear internal arrays. for(int iPad=1;iPad> 12) & 0x1FFFF; IdManuChannelToEntry& idToEntry = * const_cast(fIdToEntry); fLutEntry = idToEntry[fIdManuChannel]; fPadCharge = int(((unsigned short)(dataWord & 0xFFF)) - fLookUpTableData[fLutEntry].fPed); fCharge = 0; if(fPadCharge > fDCCut && fPadCharge > 5.0*fLookUpTableData[fLutEntry].fSigma){ // (charge > 4) is due cut out the noise level fPadData[fDataCount].fDetElemId = fLookUpTableData[fLutEntry].fDetElemId; fPadData[fDataCount].fIX = fLookUpTableData[fLutEntry].fIX; fPadData[fDataCount].fIY = fLookUpTableData[fLutEntry].fIY; fPadData[fDataCount].fRealX = fLookUpTableData[fLutEntry].fRealX; fPadData[fDataCount].fRealY = fLookUpTableData[fLutEntry].fRealY; fPadData[fDataCount].fRealZ = fLookUpTableData[fLutEntry].fRealZ; fPadData[fDataCount].fHalfPadSize = fLookUpTableData[fLutEntry].fHalfPadSize; fPadData[fDataCount].fPlane = fLookUpTableData[fLutEntry].fPlane; if ( fPadCharge < fLookUpTableData[fLutEntry].fThres ) { fCharge = (fLookUpTableData[fLutEntry].fA0)*fPadCharge; }else{ fCharge = (fLookUpTableData[fLutEntry].fA0)*(fLookUpTableData[fLutEntry].fThres) + (fLookUpTableData[fLutEntry].fA0)*(fPadCharge-fLookUpTableData[fLutEntry].fThres) + (fLookUpTableData[fLutEntry].fA1)*(fPadCharge-fLookUpTableData[fLutEntry].fThres)*(fPadCharge-fLookUpTableData[fLutEntry].fThres); } fPadData[fDataCount].fCharge = fCharge; if(fLookUpTableData[fLutEntry].fDetElemId != fPrevDetElemId){ if((*fNofFiredDetElem)>0){ fMaxFiredPerDetElem[(*fNofFiredDetElem)-1] = fDataCount; } (*fNofFiredDetElem)++; fPrevDetElemId = fLookUpTableData[fLutEntry].fDetElemId ; } // HLTDebug("buspatch : %d, detele : %d, id : %d, manu : %d, channel : %d, iX : %d, iY: %d, (X,Y) : (%f, %f), charge : %d, padsize : %f, plane : %d", // fBusPatchId,fPadData[fDataCount].fDetElemId, // fIdManuChannel,((dataWord >> 18) & 0x7FF),((dataWord >> 12) & 0x3F), // fPadData[fDataCount].fIX,fPadData[fDataCount].fIY, // fPadData[fDataCount].fRealX,fPadData[fDataCount].fRealY, // fPadData[fDataCount].fCharge,fPadData[fDataCount].fHalfPadSize,fPadData[fDataCount].fPlane); fDataCount ++; }// if charge is more than DC Cut limit condition }