#include "hshowerhit.h" #include "hshowerhitsim.h" #include "hshowerpid.h" #include "hshowerhitheader.h" #include "hshowerhitfinder.h" #include "hshowercriterium.h" #include "TArrayI.h" #include "hruntimedb.h" #include "hevent.h" #include "heventheader.h" #include "hspectrometer.h" #include "hdetector.h" #include "hshowerdetector.h" #include "hratreeext.h" #include "hcategory.h" #include "hmatrixcategory.h" #include "hlinearcatiter.h" #include "hlocation.h" #include "hshowercal.h" #include "hshowerhitfpar.h" #include "hshowergeometry.h" #include "hshowerpad.h" #include "hiterator.h" #include "hdebug.h" #include "hades.h" #include "hgeomvector.h" #include "hgeomvector2.h" #include "hspecgeompar.h" #include "hgeomvolume.h" #include "showerdef.h" #include "hgeantdef.h" #include using namespace std; ClassImp(HShowerHitFinder) //_HADES_CLASS_DESCRIPTION /////////////////////////////////////////////////////////////////////// //HShowerHitFinder // // HShowerHitFinder reconstructor searches local maxima in CAL level // The Shower task is split into several tasks as shown // in the flow diagram below. // // ---------------------- // | HShowerUnpacker | // // | (embedding mode) | \ // // | | \ ------------------ // // ---------------------- | | HGeantShower | // // | ------------------\ // // | \ // // | ------------------ \-------------> ---------------------- // // | | HGeantWire | <------------ | HShowerHitDigitizer | // // | ------------------\ ---------------------- // // | \ // // ------------- ------------------ \-------------> ----------------------- // // -- | HShowerRaw | | HShowerRawMatr | <------------ | HShowerPadDigitizer | // // | ------------- ------------------\ |( creates track objects| // // | \ | for real tracks in | // // ---------------------- ------------------ \ | embedding mode too) | // // | HShowerCalibrater | | HShowerTrack | \<------------ ----------------------- // // | (embedding mode) | ------------------\ \ // // ---------------------- \ \ ----------------------- // // | ------------------ \ ----------> | HShowerCopy | // // -------------------> | HShowerCal | \<------------ |(add charge of real hit| // // ------------------\ \ | in embedding too ) | // // \ \ ----------------------- // // ------------------ ----\---------> ----------------------- // // | HShowerHitHdr | <--\--------- | HShowerHitFinder | // // ------------------ \ ----------------------- // // ------------------ \ | // // | HShowerPID | <-----\-------| // // ------------------ \ | // // ------------------ \ | // // | HShowerHit | <--------\----| // // ------------------ < \ // // \ \ // // \-------------> ------------------------ // // | HShowerHitTrackMatcher | // // ------------------------ // // // // In the case of TRACK EMBEDDING of simulated tracks into // experimental data the real data are written by the HShowerUnpacker into // HShowerRaw category. The real hits are taken into // account by the digitizer (adding of charges). The embedding mode is recognized // automatically by analyzing the // gHades->getEmbeddingMode() flag. // Mode ==0 means no embedding // ==1 realistic embedding (first real or sim hit makes the game) // ==2 keep GEANT tracks (as 1, but GEANT track numbers will always // win against real data. besides the tracknumber the output will // be the same as in 1) // ////////////////////////////////////////////////////////////////////// HShowerHitFinder::HShowerHitFinder(const Text_t *name,const Text_t *title) : HReconstructor(name,title) { fIter = NULL; m_Loc.set(4, 0, 0, 0, 0); setFillPID(); setSortFlag(kTRUE); m_pCriterium = new HShowerCriterium; } HShowerHitFinder::HShowerHitFinder() { fIter = NULL; m_Loc.set(4, 0, 0, 0, 0); setFillPID(); setSortFlag(kTRUE); m_pCriterium = new HShowerCriterium; } HShowerHitFinder::~HShowerHitFinder(void) { if (m_pCellArr) { delete m_pCellArr; } if (fIter) delete fIter; if (m_pCriterium) delete m_pCriterium; } Bool_t HShowerHitFinder::init() { if(gHades->getCurrentEvent()->getCategory(catGeantKine)) isSim = kTRUE ; else isSim = kFALSE; TArrayI arr(4); TArrayI arr1(1); HShowerDetector *pShowerDet = (HShowerDetector*)gHades->getSetup()->getDetector("Shower"); arr[0] = pShowerDet->getMaxSectors(); arr[1] = pShowerDet->getMaxModules(); arr[2] = pShowerDet->getRows(); arr[3] = pShowerDet->getColumns(); arr1[0] = arr[0] * arr[1] * arr[2] * arr[3]; m_pCalCat=gHades->getCurrentEvent()->getCategory(catShowerCal); if (m_pCalCat) { fIter = (HIterator*)getCalCat()->MakeIterator("native"); } m_pHitCat=gHades->getCurrentEvent()->getCategory(catShowerHit); if (!m_pHitCat) { if(isSim) m_pHitCat = pShowerDet->buildLinearCat("HShowerHitSim"); else m_pHitCat = pShowerDet->buildLinearCat("HShowerHit"); if (!m_pHitCat) return kFALSE; else gHades->getCurrentEvent()->addCategory(catShowerHit, m_pHitCat, "Shower"); } m_pHitHdrCat=gHades->getCurrentEvent()->getCategory(catShowerHitHdr); if (!m_pHitHdrCat) { m_pHitHdrCat=pShowerDet->buildCategory(catShowerHitHdr); if (!m_pHitHdrCat) return kFALSE; else gHades->getCurrentEvent()->addCategory(catShowerHitHdr, m_pHitHdrCat, "Shower"); } if (IsFillPID()) { m_pPIDCat=gHades->getCurrentEvent()->getCategory(catShowerPID); if (!m_pPIDCat) { m_pPIDCat=pShowerDet->buildCategory(catShowerPID); if (!m_pPIDCat) return kFALSE; else gHades->getCurrentEvent()->addCategory(catShowerPID, m_pPIDCat, "Shower"); } } else m_pPIDCat = NULL; m_pCellArr = new HRaTreeExt(m_pCalCat, &arr); return initParameters(); } Bool_t HShowerHitFinder::initParameters() { HRuntimeDb* rtdb=gHades->getRuntimeDb(); m_pGeometry = (HShowerGeometry*)rtdb->getContainer("ShowerGeometry"); if(!m_pGeometry) return kFALSE; m_pHitFPar = (HShowerHitFPar*)rtdb->getContainer("ShowerHitFPar"); if (!m_pHitFPar) return kFALSE; m_pCriterium->setParams(m_pHitFPar); m_pHSpecGeomPar = (HSpecGeomPar *)rtdb->getContainer("SpecGeomPar"); if(!m_pHSpecGeomPar) return kFALSE; return kTRUE; } void HShowerHitFinder::setCriterium(HShowerCriterium* pCrit) { if (m_pCriterium) delete m_pCriterium; m_pCriterium = pCrit; } Bool_t HShowerHitFinder::finalize() { return kTRUE; } HShowerHitFinder& HShowerHitFinder::operator=(HShowerHitFinder &c) { return c; } Int_t HShowerHitFinder::execute() { if(!fIter) return 0; // no input HShowerCal *pCal; m_pCellArr->update(); Int_t n = 0; Int_t hits = 0; fIter->Reset(); while((pCal = (HShowerCal *)fIter->Next())) { m_Loc[0] = pCal->getSector(); m_Loc[1] = pCal->getModule(); m_Loc[2] = pCal->getRow(); m_Loc[3] = pCal->getCol(); if (lookForHit(pCal, m_Loc)) hits++; n++; } if (IsSortFlagSet()) { m_pHitCat->sort(); if(m_pPIDCat) m_pPIDCat->sort(); } return 0; } Bool_t HShowerHitFinder::lookForHit(HShowerCal* cal, HLocation& fLoc) { #if DEBUG_LEVEL>2 gDebuger->enterFunc("HShowerHitFinder::execute"); gDebuger->message("Cal cat points to %p",cal); #endif HShowerHit *hit = NULL; HShowerPID *pid = NULL; Bool_t isHit = kFALSE; if (cal && isLocalMax(fLoc)) { Float_t fCharge, fShower = 0.0; Int_t ret = 0; hit=(HShowerHit *)m_pHitCat->getNewSlot(fLoc); if (hit) { if(isSim){ hit = new (hit) HShowerHitSim; } else { hit = new (hit) HShowerHit; } fillSums(hit, fLoc); //calcCoord(hit, fLoc); calcCoordWithSigma(hit, fLoc,1); fCharge = cal->getCharge(); hit->setCharge(fCharge); hit->setSector(cal->getSector()); hit->setModule(cal->getModule()); hit->setRow(cal->getRow()); hit->setCol(cal->getCol()); if (IsSortFlagSet()) hit->calcAddress(); hit->copyToTrueAddress(); // later used in low efficiency mode! fShower = m_pCriterium->showerCriterium(hit, ret,m_pHitFPar); hit->setShower(fShower); if (m_pPIDCat && fShower > 0) { pid=(HShowerPID *)m_pPIDCat->getNewSlot(fLoc); if (pid) { pid=new(pid) HShowerPID; fillPID(hit, pid); } } isHit = kTRUE; updateLocalMax(fLoc); } } #if DEBUG_LEVEL>2 gDebuger->leaveFunc("HShowerHitFinder::execute"); #endif updateFiredCells(fLoc); return isHit; } void HShowerHitFinder::fillPID(HShowerHit* hit, HShowerPID* pid) { Float_t fX, fY, fZ; Float_t fR, fPhi, fTheta; pid->setCharge(hit->getCharge()); pid->setSector(hit->getSector()); pid->setModule(hit->getModule()); pid->setRow(hit->getRow()); pid->setCol(hit->getCol()); pid->setAddress(hit->getAddress()); pid->setShower(hit->getShower()); hit->getLabXYZ(&fX, &fY, &fZ); pid->setXYZ(fX, fY, fZ); hit->getSphereCoord(&fR, &fPhi, &fTheta); pid->setSphereCoord(fR, fPhi, fTheta); // if (IsSortFlagSet()) // getPIDCat()->sort(); } void HShowerHitFinder::fillSums(HShowerHit* hit, HLocation &fLoc) { HLocation fTmpLoc; Int_t nModules = m_pHitFPar->getModules(); Float_t sum, var; Int_t cs; fTmpLoc = fLoc; for(fTmpLoc[1] = 0; fTmpLoc[1] < nModules; fTmpLoc[1]++) { sum = calculateSum(fTmpLoc, 1, &cs); var = calculateVar(fTmpLoc, 1, sum/9.0); hit->setSum(fTmpLoc[1], sum); hit->setVar(fTmpLoc[1], var); hit->setClusterSize(fTmpLoc[1], cs); } calculateSum(fLoc, 1, &cs); if (cs > 1) updateClusters(fLoc); if (nModules == 3) { fTmpLoc[1] = 2; hit->setSum25(calculateSum(fTmpLoc, 2)); } hit->updateCalc(); } void HShowerHitFinder::calcCoord(HShowerHit* hit, HLocation &fLoc) { HGeomVector vLocalCoord; HGeomVector2 vLabCoord; HGeomVector targetPos = m_pHSpecGeomPar->getTarget(0)->getTransform().getTransVector(); m_pGeometry->getLocalCoord(fLoc, vLocalCoord); m_pGeometry->getLabCoord(fLoc, vLabCoord); m_pGeometry->getSphereCoord(fLoc, vLabCoord, &targetPos); hit->setXY(vLocalCoord.getX(), vLocalCoord.getY()); hit->setLabXYZ(vLabCoord.getX(), vLabCoord.getY(), vLabCoord.getZ()); hit->setSphereCoord(vLabCoord.getRad(), vLabCoord.getPhi(), vLabCoord.getTheta()); } //////////////////////////////////////////////////////////////////////////// // // HShowerHitFinder::calcCoordWithSigma // Calculate the Position of the hit and SigmaX,SigmaY. // SigmaX,SigmaY are determined based on the geometry of the pad. // Weighted mean value for X,Y position determination is used. // Procedure: // 1.) calculate weighted X_Loc, Y_Loc in local coord. system // 2.) determine Sigma_X, Sigma_Y (based on geometry of the pad) // 3.) transform X_Loc, Y_Loc -> X_Lab, Y_Lab, Z_Lab // 4.) transform X_Lab, Y_Lab, Z_Lab -> Sphere Coord. Sys. // /////////////////////////////////////////////////////////////////////////// void HShowerHitFinder::calcCoordWithSigma(HShowerHit* hit, HLocation &fLoc, Int_t nRange) { HGeomVector vLocalCoord; HShowerPadTab *pPads = m_pGeometry->getPadParam(fLoc[1]); // Access to the pad plane (module== fLoc[1]) Int_t rowL, rowU, colL, colR,iRow1, iRow2, iCol1, iCol2; Int_t row = fLoc[2]; Int_t col = fLoc[3]; // Check pads around Loc_Max m_pHitFPar->getRowBord(fLoc[0], fLoc[1], &rowL, &rowU); m_pHitFPar->getColBord(fLoc[0], fLoc[1], &colL, &colR); if((row < rowL) || (row > rowU) || (col < colL) || (col > colR)) cout<<"--We are outside of the detector. "< rowU) iRow2 = rowU; if((iCol1 = col - nRange) < colL) iCol1 = colL; if((iCol2 = col + nRange) > colR) iCol2 = colR; Float_t fSum = 0; Float_t fWSumX_Loc = 0; Float_t fWSumY_Loc = 0; Float_t sigma_X_tmp = 0; Float_t sigma_X = 0; Float_t sigma_Y_tmp = 0; Float_t sigma_Y = 0; HLocation fTmpLoc; fTmpLoc = fLoc; HShowerCal* pCell; ///////////////////////////////////////////////////////////// // actual calculation of weighted_X_Loc and weighted_Y_Loc // ///////////////////////////////////////////////////////////// for(fTmpLoc[2] = iRow1; fTmpLoc[2] <= iRow2; fTmpLoc[2]++) for(fTmpLoc[3] = iCol1; fTmpLoc[3] <= iCol2; fTmpLoc[3]++) { pCell = (HShowerCal*)m_pCellArr->getObject(fTmpLoc); if((float)pCell->getCharge()<=0.0) continue; // position of the pad in Local coord. sys. m_pGeometry->getLocalCoord(fTmpLoc, vLocalCoord); //position in Local Sys. fSum = fSum + pCell->getCharge(); fWSumX_Loc = fWSumX_Loc + ((pCell->getCharge()) * (vLocalCoord.getX()) ); //weighted X_Loc fWSumY_Loc = fWSumY_Loc + ((pCell->getCharge()) * (vLocalCoord.getY()) ); //weighted Y_Loc // Determine Sigma_X (geometry), pad_size/sqrt(12), sqrt(12)-standard deviation of a uniform density Float_t fXLD = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getXld_mm(); Float_t fXRD = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getXrd_mm(); Float_t fXLU = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getXlu_mm(); Float_t fXRU = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getXru_mm(); Float_t x1=fabs(fXLU-fXRD); Float_t x2=fabs(fXLD-fXRU); sigma_X_tmp = (x1>x2)? x1:x2; sigma_X_tmp = sigma_X_tmp/sqrt(12.); sigma_X_tmp = ((float)pCell->getCharge()) * sigma_X_tmp; sigma_X_tmp = sigma_X_tmp * sigma_X_tmp; // so I have Q(i)^2 * Sigma(i)^2 sigma_X = sigma_X + sigma_X_tmp; // and here: Q(i)^2 * Sigma(i)^2 + Q(j)^2 * Sigma(j)^2 // Determine Sigma_Y (geometry) Float_t fYLD = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getYld_mm(); Float_t fYRD = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getYrd_mm(); Float_t fYLU = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getYlu_mm(); Float_t fYRU = pPads->getPad(fTmpLoc[2],fTmpLoc[3])->getYru_mm(); x1=fabs(fYLU-fYRD); x2=fabs(fYLD-fYRU); sigma_Y_tmp = (x1>x2)? x1:x2; sigma_Y_tmp = sigma_Y_tmp/sqrt(12.); sigma_Y_tmp = ((float)pCell->getCharge()) * sigma_Y_tmp; sigma_Y_tmp = sigma_Y_tmp * sigma_Y_tmp; // so I have Q(i)^2 * Sigma(i)^2 sigma_Y = sigma_Y + sigma_Y_tmp; //and here: Q(i)^2 * Sigma(i)^2 + Q(j)^2 * Sigma(j)^2 } if(fSum<=0.0) cout<<"--Shw_Hit_F-- Chagrge = 0, should not happen"<transVectToLab(fLoc, vv, vvv); hit->setXY((float)fWSumX_Loc/fSum,(float)fWSumY_Loc/fSum ); // Hit in Local System. hit->setSigmaXY(sqrt(sigma_X/(fSum*fSum)),sqrt(sigma_Y/(fSum*fSum))); hit->setLabXYZ(vvv.getX(), vvv.getY(), vvv.getZ()); // Hit in Lab. Sysstem HGeomVector2 sphereOut; // Sphere coordin. HGeomVector targetPos = m_pHSpecGeomPar->getTarget(0)->getTransform().getTransVector(); m_pGeometry->transLabToSphereCoord(vvv, sphereOut, &targetPos); hit->setSphereCoord(sphereOut.getRad(), sphereOut.getPhi(), sphereOut.getTheta()); /* if(0){ // create control file Int_t nEvt=gHades->getCurrentEvent()->getHeader()->getEventSeqNumber(); // get Event Seq. Nb. ofstream out1("results.txt", ios::app); HGeomVector vLocalCoord_Tmp; m_pGeometry->getLocalCoord(fLoc, vLocalCoord_Tmp); HGeomVector2 vLabCoord; m_pGeometry->getSphereCoord(fLoc, vLabCoord, &targetPos); out1<<" evnt_Nb "<getClusterSize(fLoc[1])<< " old_X "<<(float)(vLocalCoord_Tmp.getX())<<" new_X "<< (float)fWSumX_Loc/fSum<< " old_Y "<<(float)(vLocalCoord_Tmp.getY())<<" new_Y "<< ((float)fWSumY_Loc/fSum)<< " sig_X "<getSigmaX()<<" sigY "<getSigmaY()<< //" sig_X "<getObject(fLoc); Float_t charge, charge1; Int_t rowL, rowU, colL, colR, iRow1, iRow2, iCol1, iCol2;; Int_t row = fLoc[2]; Int_t col = fLoc[3]; HLocation fTmpLoc; m_pHitFPar->getRowBord(fLoc[0], fLoc[1], &rowL, &rowU); m_pHitFPar->getColBord(fLoc[0], fLoc[1], &colL, &colR); if((row < rowL) || (row > rowU) || (col < colL) || (col > colR)) return 0; if((iRow1 = row - 1) < rowL) iRow1 = rowL; if((iRow2 = row + 1) > rowU) iRow2 = rowU; if((iCol1 = col - 1) < colL) iCol1 = colL; if((iCol2 = col + 1) > colR) iCol2 = colR; fTmpLoc = fLoc; charge = pRefCell->getCharge(); for(fTmpLoc[2] = iRow1; fTmpLoc[2] <= iRow2; fTmpLoc[2]++) for(fTmpLoc[3] = iCol1; fTmpLoc[3] <= iCol2; fTmpLoc[3]++) { if((fTmpLoc[2] == row) && (fTmpLoc[3] == col)) continue; pCell = (HShowerCal*)m_pCellArr->getObject(fTmpLoc); if (pCell->isLocalMax() != 0) return 0; charge1 = pCell->getCharge(); if(charge < charge1) return 0; } pRefCell->setLocalMax(); return 1; } Float_t HShowerHitFinder::calculateVar(HLocation &fLoc, Int_t nRange, Float_t avg) { Int_t rowL, rowU, colL, colR, iRow1, iRow2, iCol1, iCol2; Int_t row = fLoc[2]; Int_t col = fLoc[3]; m_pHitFPar->getRowBord(fLoc[0], fLoc[1], &rowL, &rowU); m_pHitFPar->getColBord(fLoc[0], fLoc[1], &colL, &colR); if((row < rowL) || (row > rowU) || (col < colL) || (col > colR)) return 0.0f; if((iRow1 = row - nRange) < rowL) iRow1 = rowL; if((iRow2 = row + nRange) > rowU) iRow2 = rowU; if((iCol1 = col - nRange) < colL) iCol1 = colL; if((iCol2 = col + nRange) > colR) iCol2 = colR; Float_t fDiff = 0.0, fVar = 0.0; HLocation fTmpLoc; fTmpLoc = fLoc; HShowerCal* pCell; // for(fTmpLoc[2] = row - nRange; fTmpLoc[2] <= row + nRange; fTmpLoc[2]++) // for(fTmpLoc[3] = col - nRange; fTmpLoc[3] <= col + nRange; fTmpLoc[3]++) for(fTmpLoc[2] = iRow1; fTmpLoc[2] <= iRow2; fTmpLoc[2]++) for(fTmpLoc[3] = iCol1; fTmpLoc[3] <= iCol2; fTmpLoc[3]++) { pCell = (HShowerCal*)m_pCellArr->getObject(fTmpLoc); fDiff += (pCell->getCharge() - avg); fVar += fDiff * fDiff; } return fVar; } Float_t HShowerHitFinder::calculateSum(HLocation &fLoc, Int_t nRange, Int_t* pncs) { Int_t nThreshold = m_pHitFPar->getThreshold(); Int_t rowL, rowU, colL, colR,iRow1, iRow2, iCol1, iCol2; Int_t row = fLoc[2]; Int_t col = fLoc[3]; m_pHitFPar->getRowBord(fLoc[0], fLoc[1], &rowL, &rowU); m_pHitFPar->getColBord(fLoc[0], fLoc[1], &colL, &colR); if((row < rowL) || (row > rowU) || (col < colL) || (col > colR)) return 0.0f; if((iRow1 = row - nRange) < rowL) iRow1 = rowL; if((iRow2 = row + nRange) > rowU) iRow2 = rowU; if((iCol1 = col - nRange) < colL) iCol1 = colL; if((iCol2 = col + nRange) > colR) iCol2 = colR; Float_t fSum = 0; HLocation fTmpLoc; fTmpLoc = fLoc; Int_t cs = 0; HShowerCal* pCell; // for(fTmpLoc[2] = row - nRange; fTmpLoc[2] <= row + nRange; fTmpLoc[2]++) // for(fTmpLoc[3] = col - nRange; fTmpLoc[3] <= col + nRange; fTmpLoc[3]++) for(fTmpLoc[2] = iRow1; fTmpLoc[2] <= iRow2; fTmpLoc[2]++) for(fTmpLoc[3] = iCol1; fTmpLoc[3] <= iCol2; fTmpLoc[3]++) { pCell = (HShowerCal*)m_pCellArr->getObject(fTmpLoc); fSum += pCell->getCharge(); if (pCell->getCharge() > nThreshold) cs++; } if (pncs) *pncs = cs; return fSum; } HShowerHitHeader* HShowerHitFinder::getHitHeader(HLocation &fLoc) { HLocation fTmpLoc; fTmpLoc.set(2, fLoc.getIndex(0), fLoc.getIndex(1)); fTmpLoc[0] = fLoc[0]; fTmpLoc[1] = fLoc[1]; HShowerHitHeader *pHitHdr = (HShowerHitHeader*) ((HMatrixCategory*)m_pHitHdrCat)->getObject(fTmpLoc); if (!pHitHdr) { pHitHdr = (HShowerHitHeader*)((HMatrixCategory*)m_pHitHdrCat) ->getSlot(fTmpLoc); if (pHitHdr != NULL) { pHitHdr = new(pHitHdr) HShowerHitHeader; pHitHdr->setSector(fTmpLoc[0]); pHitHdr->setModule(fTmpLoc[1]); } } return pHitHdr; } void HShowerHitFinder::updateClusters(HLocation &fLoc) { HShowerHitHeader *pHitHdr = getHitHeader(fLoc); if (pHitHdr) pHitHdr->incClusters(); } void HShowerHitFinder::updateLocalMax(HLocation &fLoc) { HShowerHitHeader *pHitHdr = getHitHeader(fLoc); if (pHitHdr) pHitHdr->incLocalMax(); } void HShowerHitFinder::updateFiredCells(HLocation &fLoc) { HShowerHitHeader *pHitHdr = getHitHeader(fLoc); if (pHitHdr) pHitHdr->incFiredCells(); }