//_HADES_CLASS_DESCRIPTION ////////////////////////////////////////////////////////////////////////////// // // HParticleBtRingF // // This class predicts fired pads based on theta and phi angle of a track. // Predicted pads and rich cal hits are stored and used to search for // fired pads matching to prediction area. // // // // ////////////////////////////////////////////////////////////////////////////// #include "hparticlebtringf.h" #include "TF2.h" #include "hcategory.h" #include "hparticlecand.h" #include "heventheader.h" #include "hrichcal.h" #include "hades.h" #include "hruntimedb.h" #include "hcategorymanager.h" #include "hparticlebtangletrafo.h" // ---------------------------------------------------------------- ClassImp(HParticleBtRingF) // ---------------------------------------------------------------- HParticleBtRingF::HParticleBtRingF( ){ } HParticleBtRingF::~HParticleBtRingF( ){ } Bool_t HParticleBtRingF::init() { fNSector = 6; fBtPar =(HParticleBtPar*) gHades->getRuntimeDb()->getContainer("ParticleBtPar"); if(!fBtPar) { Error ("init()", "Retrieve ZERO pointer for ParticleBtPar!"); return kFALSE; } //Init constants fNVertex = fBtPar->getNVertex(); fNParMean = fBtPar->getNParMean(); fNParSigma = fBtPar->getNParSigma(); fNRingSeg = fBtPar->getNRingSegments(); fNRingSegStep = fBtPar->getRingSegStep(); fNRingSegOffset = fBtPar->getRingSegOffset(); fNRichSeg = fBtPar->getNRichSeg(); fRichSegBorderX = fBtPar->getRichSegBorderX(); fRichSegBorderY = fBtPar->getRichSegBorderY(); fThetaAngleMin = fBtPar->getThetaAngleMin(); fThetaAngleMax = fBtPar->getThetaAngleMax(); fSizeMatrix = fBtPar->getSizeMatrix(); memcpy(&fChargeLimit[0] ,fBtPar->getChargeThres() , sizeof(fChargeLimit) ); memcpy(&fChargeLimitMaximum[0] ,fBtPar->getChargeThresMax(), sizeof(fChargeLimitMaximum) ); fNSigma = fBtPar->getSigmaRange(); fMaxSigmaRange = fBtPar->getSigmaRangeMax(); fVertexPosMin = fBtPar->getVertexPosMin(); fVertexPosMax = fBtPar->getVertexPosMax(); fVertexStep = fBtPar->getVertexStep(); memcpy(&fPhiOff[0] , fBtPar->getPhiOffset() , sizeof(fPhiOff ) ); memcpy(&fPhiOff2[0], fBtPar->getPhiOffset2(), sizeof(fPhiOff2 ) ); fPhiOffsetPar = fBtPar->getOffsetPar(); //TODO TF2 fPol2DMean.resize(fNRingSeg); fPol2DSigma.resize(fNRingSeg); for(Int_t bin=0; binSetParameter(par, fBtPar->getTF2ParMean( bin, vertex, par)); for(Int_t par=0; parSetParameter(par, fBtPar->getTF2ParSigma(bin, vertex, par)); } fRad2Deg.push_back(TMath::DegToRad()*(bin*fNRingSegStep+fNRingSegOffset)); fRad2DegX.push_back(TMath::Sin(fRad2Deg[bin])); fRad2DegY.push_back(-1*TMath::Cos(fRad2Deg[bin])); } fFiredPads.resize( fNSector , vector( fNRichSeg , 0 ) ); return kTRUE; } //-------------------------------------------------------------------- // -- Helper functions void HParticleBtRingF::addressToColRow(const Int_t address, Int_t &sec, Int_t &row, Int_t &col) { //Converts richCal address to sector, column and row numbers sec = (Int_t)address/10000; row = (Int_t)(address-sec*10000)/100; col = address-sec*10000-row*100; sec = sec==6 ? 0 : sec; } Int_t HParticleBtRingF::correctPhi(const Int_t sec, const Float_t phi) { //Correct phi angle if phi angle disagrees with sector number //Agreement is necessary for transformations from angle to RICH pad plane if(sec==0){ if(phi>120) return 120; else if(phi<60) return 60; else return phi; }else if(sec==1){ if(phi<120) return 120; else if(phi>180) return 180; else return phi; }else if(sec==2){ if(phi<180) return 180; else if(phi>240) return 240; else return phi; }else if(sec==3){ if(phi<240) return 240; else if(phi>300) return 300; else return phi; }else if(sec==4){ if(phi>360 || (phi>=0 && phi < 60)) return 360; else if(phi<300) return 300; else return phi; }else{ //sec==5 if(phi>300 || phi<0) return 0; else if(phi >60 && phi<=300) return 60; else return phi; } } void HParticleBtRingF::sortElements(Double_t &entry1 ,Double_t &entry2) { //Sorts 2 Double_t values by size (ascending) if(entry1 >entry2){ Float_t tmp = entry1; entry1 = entry2; entry2 = tmp; } } Int_t HParticleBtRingF::getVertexNum(const Float_t vertex) { //Returns vertex number to given z-vertex position for(Int_t i=0; i < fNVertex;i++){ if((vertex>-fVertexPosMin+(i*fVertexStep))&& (vertex<-fVertexPosMax+(i*fVertexStep))) return i; } if(vertex<=fVertexPosMin) return 0; else //v>-2.5 return 14; } //-------------------------------------------------------------------- // -- Called once per event void HParticleBtRingF::fillRichCal(HCategory* catRichCal) { //Fills rich cal hit addresses and charges into vectors //Bad hits and hits from bad events are exluded Int_t sizeRichCal = catRichCal->getEntries(); HRichCal* richCal; Int_t sec; Int_t col; Int_t row; for(Int_t l=0; l < sizeRichCal;l++){ richCal = HCategoryManager::getObject(richCal,catRichCal,l); if(!richCal) continue; if(!richCal->getIsCleanedSector() && !richCal->getIsCleanedHigh()){ addressToColRow(richCal->getAddress(),sec, row, col); if(richCal->getCharge() > fChargeLimit[sec]){ fRichHitAdd.push_back(richCal->getAddress()); fRichHitCharge.push_back(richCal->getCharge()); } if(row ringMatrixVec; xPad -= fSizeMatrix*0.5; yPad -= fSizeMatrix*0.5; address = 10000 * (sec ? sec : 6) + 100 * yPad + xPad; for(Int_t posX=0; posX getTheta(); Double_t phi = cand->getPhi(); Int_t sec = cand->getSector(); Float_t vertexZ = vertex.getZ(); Int_t vertexNum = getVertexNum(vertexZ); Double_t radius; Double_t sigma; Double_t sigmaX0; Double_t sigmaY0; Double_t radiusSeg0; Double_t sigmaSeg0; Double_t sigmaX0Seg0; Double_t sigmaY0Seg0; Double_t radiusNext; Double_t sigmaNext; Double_t sigmaX0Next; Double_t sigmaY0Next; Double_t sigNum; vector predictionVec; phi = correctPhi(sec,phi); if(theta>fThetaAngleMin && theta < fThetaAngleMax && vertexZ > -80. && vertexZ < 20.){ //Correct theta and phi angle Double_t thetaCor; Double_t phiCor; if(fAngleCor.alignRichRing(theta,phi, thetaCor, phiCor)){ theta = thetaCor; phi = phiCor; } thetaCor = richAngles.zTheta2dTheta(vertex.getZ(),theta, -(phi-fPhiOff[sec])); theta-=thetaCor; //Transform angles to x,y Padplane coordinates Double_t xPad = richAngles.angles2xPad(theta,-(phi-fPhiOff[sec])); Double_t yPad = richAngles.angles2yPad(theta,-(phi-fPhiOff[sec])); Double_t posX = richAngles.angles2x(theta,-(phi-fPhiOff[sec])); Double_t posY = richAngles.angles2y(theta,-(phi-fPhiOff[sec])); //z vertex correction Double_t yPadCorr = richAngles.zTheta2dYPad(vertex.getZ(),theta,-(phi-fPhiOff2[sec])); Double_t xPadCorr = richAngles.zTheta2dXPad(); Double_t posXCorr = richAngles.zTheta2dY(vertex.getZ(),theta,-(phi-fPhiOff2[sec])); Double_t posYCorr = richAngles.zTheta2dX(); posX += posXCorr; posY += posYCorr; xPad += xPadCorr; yPad += yPadCorr; //Coordinates for n sigma areas Double_t sigmaX[2] = {0.,0.}; Double_t sigmaY[2] = {0.,0.}; Double_t padX[2] = {0.,0.}; Double_t padY[2] = {0.,0.}; Int_t padXRound[2] = {0,0}; Int_t padYRound[2] = {0,0}; sigNum=fNSigma; //Loop over different phi angles of the ring for( Int_t bin = 0; bin < fNRingSeg; bin++ ) { //Find neighbouring segments Int_t binPrev; Int_t binNext; if(bin==0){ binPrev = fNRingSeg-1; binNext = 2; } else if(bin==fNRingSeg-1){ binPrev = fNRingSeg-2; binNext = 0; } else{ binPrev = bin-1; binNext = bin+1; } radius = fPol2DMean[bin][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigma = fPol2DSigma[bin][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaX0 = fRad2DegX[bin]*sigNum*sigma; sigmaY0 = fRad2DegY[bin]*sigNum*sigma; //Smooth ring prediction area width neighbour sigma and radius radiusSeg0 = fPol2DMean[binPrev][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaSeg0 = fPol2DSigma[binPrev][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaX0Seg0 = fRad2DegX[binPrev]*sigNum*sigmaSeg0; sigmaY0Seg0 = fRad2DegY[binPrev]*sigNum*sigmaSeg0; radiusNext = fPol2DMean[binNext][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaNext = fPol2DSigma[binNext][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaX0Next = fRad2DegX[binNext]*sigNum*sigmaNext; sigmaY0Next = fRad2DegY[binNext]*sigNum*sigmaNext; radius = (radius + radiusSeg0 + radiusNext )/3.; if(sigmaEval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaSeg0 = fPol2DSigma[0][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaX0Seg0 = fRad2DegX[0]*sigNum*sigmaSeg0; sigmaY0Seg0 = fRad2DegY[0]*sigNum*sigmaSeg0; radius = (radius + radiusSeg0 )/2.; if(sigmaEval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaSeg0 = fPol2DSigma[fNRingSeg-1][vertexNum]->Eval(theta,-(phi-fPhiOff[sec])+fPhiOffsetPar); sigmaX0Seg0 = fRad2DegX[fNRingSeg-1]*sigNum*sigmaSeg0; sigmaY0Seg0 = fRad2DegY[fNRingSeg-1]*sigNum*sigmaSeg0; radius = (radius + radiusSeg0 )/2.; if(sigmagetIndex()); fillMatrix(TMath::Floor(xPad),TMath::Floor(yPad),sec); } else { predictionVec.clear(); predictionVec.push_back(-1); fPrediction.push_back(predictionVec); fTrackTheta.push_back(-1); fTrackPhi.push_back(-1); fTrackVertex.push_back(-1); fTrackSec.push_back(-1); fTrackPCandIdx.push_back(-1); fRingMatrix.push_back(predictionVec); } } Float_t HParticleBtRingF::getRingMatrix(const Int_t trackNo) { //Returns fraction of pads fired on ring prediction compared to pads fired in ring matrix around ring center Int_t padsFired = 0; Int_t padsFiredPred = 0; for(UInt_t i=0; isetPrediction ( prediction ); btRingInfo->setRingMatrix ( ringMatrix ); btRingInfo->setRichHitAdd ( richHitAdd ); btRingInfo->setRichHitCharge ( richHitCharge); btRingInfo->setTrackTheta ( trackTheta ); btRingInfo->setTrackPhi ( trackPhi ); btRingInfo->setTrackVertex ( trackVertex ); btRingInfo->setTrackSec ( trackSec ); btRingInfo->setTrackPCandIdx ( trackPCandIdx); return kTRUE; }else{ return kFALSE; } } Bool_t HParticleBtRingF::hasNoisyRichSeg(Bool_t *trackInSec) { Bool_t secUsed[6]; for(Int_t i=0;i<6;i++) secUsed[0+i] = trackInSec[i]; Bool_t isNoisy = kFALSE; for(Int_t i = 0; i < fNSector; i++){ for(Int_t j = 0; j < fNRichSeg; j++){ if( fFiredPads[i][j] > 100 && secUsed[i]) isNoisy = kTRUE; } } if(isNoisy) cout << "Noisy RICH segment detected and removed" << endl; return isNoisy; }