// C++ includes #include #include #include using namespace std; // ROOT includes #include "TClonesArray.h" #include "TCanvas.h" #include "TH2F.h" #include "TArc.h" #include "TDatabasePDG.h" #include "TRandom.h" // CBM includes #include "CbmMCPoint.h" #include "CbmRootManager.h" #include "CbmSttHit.h" #include "CbmSttTrack.h" #include "CbmSttTrackFinderIdeal.h" #include "CbmSttHoughAccumulatorNew.h" #include "CbmSttHoughDefines.h" // ----- Default constructor ------------------------------------------- CbmSttTrackFinderIdeal::CbmSttTrackFinderIdeal() { rootoutput = kFALSE; fMCTrackArray = NULL; fVerbose = 1; } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ CbmSttTrackFinderIdeal::CbmSttTrackFinderIdeal(Int_t verbose) { fMCTrackArray = NULL; fVerbose = verbose; if (verbose > 2) rootoutput = kTRUE; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- CbmSttTrackFinderIdeal::~CbmSttTrackFinderIdeal() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- void CbmSttTrackFinderIdeal::Init() { // Get and check CbmRootManager CbmRootManager* ioman = CbmRootManager::Instance(); if (!ioman) { cout << "-E- CbmSttTrackFinderIdeal::Init: " << "RootManager not instantised!" << endl; return; } // Get MCTrack array fMCTrackArray = (TClonesArray*) ioman->ActivateBranch("MCTrack"); if ( ! fMCTrackArray) { cout << "-E- CbmSttTrackFinderIdeal::Init: No MCTrack array!" << endl; return; } } // ------------------------------------------------------------------------- // ----- Public method DoFind ------------------------------------------ Int_t CbmSttTrackFinderIdeal::DoFind(TClonesArray* trackArray) { // Check pointers if ( !fMCTrackArray ) { cout << "-E- CbmSttTrackFinderIdeal::DoFind: " << "MCTrack array missing! " << endl; return -1; } if (fHitCollectionList.GetEntries() == 0) { cout << "-E- CbmSttTrackFinderIdeal::DoFind: " << "No hit arrays present, call AddHitCollection() first (at least once)! " << endl; return -1; } if (fPointCollectionList.GetEntries() == 0) { cout << "-E- CbmSttTrackFinderIdeal::DoFind: " << "No point arrays present, call AddHitCollection() first (at least once)! " << endl; return -1; } if ( !trackArray ) { cout << "-E- CbmSttTrackFinderIdeal::DoFind: " << "Track array missing! " << endl; return -1; } TCanvas *finderCanvas; if (rootoutput) { TH2F myRange("range", "range", 100, -50., 50., 100, -50., 50); finderCanvas = new TCanvas("findercanvas", "findercanvas", 800, 600); myRange.DrawCopy(); plotAllStraws(); } // Initialise control counters Int_t nNoMCTrack = 0; Int_t nNoTrack = 0; Int_t nNoSttPoint = 0; Int_t nNoSttHit = 0; // Create pointers to hit and SttPoint CbmSttHit* pMhit = NULL; CbmMCPoint* pMCpt = NULL; CbmMCTrack* pMCtr = NULL; CbmSttTrack* pTrck = NULL; // Number of STT hits Int_t nHits = 0; for (Int_t hitListCounter = 0; hitListCounter < fHitCollectionList.GetEntries(); hitListCounter++) { nHits += ((TClonesArray *)fHitCollectionList.At(hitListCounter))->GetEntriesFast(); } // Declare some variables outside the loops Int_t ptIndex = 0; // MCPoint index Int_t mcTrackIndex = 0; // MCTrack index Int_t trackIndex = 0; // STTTrack index // Create STL map from MCtrack index to number of valid SttHits map > hitMap; // Loop over hits for (Int_t iHit = 0; iHit < nHits; iHit++) { pMhit = GetHitFromCollections(iHit); if ( ! pMhit ) continue; ptIndex = pMhit->GetRefIndex(); if (rootoutput) { TArc myArc; myArc.SetFillStyle(0); myArc.SetLineColor(1); if (pMhit->GetIsochrone() == 0) myArc.DrawArc(pMhit->GetX(), pMhit->GetY(), 1.0); else myArc.DrawArc(pMhit->GetX(), pMhit->GetY(), pMhit->GetIsochrone()); } if (ptIndex < 0) continue; // fake or background hit pMCpt = GetPointFromCollections(iHit); if ( ! pMCpt ) continue; mcTrackIndex = pMCpt->GetTrackID(); Double_t wireX = pMhit->GetX(), wireY = pMhit->GetY(); (hitMap[mcTrackIndex])[wireX * wireX + wireY * wireY]++; } // Create STL map from MCTrack index to SttTrack index map correlationMap, trackMap; // Create STTTracks for reconstructable MCTracks Int_t nMCacc = 0; // accepted MCTracks (more than 3 points) Int_t nTracks = 0; // reconstructable MCTracks Int_t nMCTracks = fMCTrackArray->GetEntriesFast(); for (Int_t iMCTrack = 0; iMCTrack < nMCTracks; iMCTrack++) { pMCtr = (CbmMCTrack*) fMCTrackArray->At(iMCTrack); if ( ! pMCtr ) continue; if (hitMap[iMCTrack].size() < 3) continue; // temporary hack, fix this in monte carlo .... if (TDatabasePDG::Instance()->GetParticle(pMCtr->GetPdgCode()) == NULL) continue; if (!(fabs(TDatabasePDG::Instance()->GetParticle(pMCtr->GetPdgCode())->Charge() / 3.0) > 0.)) continue; nMCacc++; new((*trackArray)[nTracks]) CbmSttTrack(); if (fVerbose) cout << "-I- CbmSttTrackFinderIdeal: STTTrack " << nTracks << " created from MCTrack " << iMCTrack << " (" << pMCtr->GetSttPoints() << " STTPoints)" << endl; correlationMap[nTracks] = iMCTrack; trackMap[iMCTrack] = nTracks++; } // Loop over hits. Get corresponding MCPoint and MCTrack index for (Int_t iHit = 0; iHit < nHits; iHit++) { pMhit = GetHitFromCollections(iHit); if ( ! pMhit ) { cout << "-E- CbmSttTrackFinderIdeal::DoFind: Empty slot " << "in HitArray at position " << iHit << endl; nNoSttHit++; continue; } if (pMhit->GetIsochrone() == 0.) { cout << "mvd hit: " << pMhit->GetX() << " " << pMhit->GetY() << " " << pMhit->GetZ() << endl; } ptIndex = pMhit->GetRefIndex(); if (ptIndex < 0) continue; // fake or background hit pMCpt = GetPointFromCollections(iHit); if ( ! pMCpt ) { nNoSttPoint++; continue; } mcTrackIndex = pMCpt->GetTrackID(); if (mcTrackIndex < 0 || mcTrackIndex > nMCTracks) { cout << "-E- CbmSttTrackFinderIdeal::DoFind: " << "MCTrack index out of range. " << mcTrackIndex << " " << nMCTracks << endl; nNoMCTrack++; continue; } if (trackMap.find(mcTrackIndex) == trackMap.end()) continue; trackIndex = trackMap[mcTrackIndex]; pTrck = (CbmSttTrack*) trackArray->At(trackIndex); if ( ! pTrck ) { cout << "-E- CbmSttTrackFinderIdeal::DoFind: " << "No SttTrack pointer. " << iHit << " " << ptIndex << " " << mcTrackIndex << " " << trackIndex << endl; nNoTrack++; continue; } pTrck->AddHit(iHit, pMhit); if (rootoutput) { TArc myArc; myArc.SetFillStyle(0); myArc.SetLineColor(trackIndex + 2); if (pMhit->GetIsochrone() == 0) myArc.DrawArc(pMhit->GetX(), pMhit->GetY(), 1.0); else myArc.DrawArc(pMhit->GetX(), pMhit->GetY(), pMhit->GetIsochrone()); } if (fVerbose > 2) cout << "Hit " << iHit << " from STTPoint " << ptIndex << " (MCTrack " << mcTrackIndex << ") added to STTTrack " << trackIndex << endl; } for (Int_t trackTeller = 0; trackTeller < nTracks; trackTeller++) { // loop over pTrck = (CbmSttTrack*) trackArray->At(trackTeller); if (pTrck != NULL) { Double_t dSeed, rSeed, phiSeed; GetTrack(dSeed, phiSeed, rSeed, correlationMap[trackTeller]); Double_t xSeed = (dSeed + rSeed) * cos(phiSeed), ySeed = (dSeed + rSeed) * sin(phiSeed); Double_t xSeed_old = xSeed, ySeed_old = ySeed, rSeed_old = rSeed; // ZoomTrack(dSeed, phiSeed, rSeed, pTrck); xSeed = (dSeed + rSeed) * cos(phiSeed); ySeed = (dSeed + rSeed) * sin(phiSeed); if (rootoutput) { TArc myArc; myArc.SetFillStyle(0); //myArc.SetLineColor(trackTeller + 2); //myArc.DrawArc(xSeed_old, ySeed_old, rSeed_old); myArc.SetLineColor(trackTeller + 2); myArc.DrawArc(xSeed, ySeed, rSeed); } pTrck->GetParamLast()->SetX(dSeed); pTrck->GetParamLast()->SetY(phiSeed); pTrck->GetParamLast()->SetZ(0.); pTrck->GetParamLast()->SetTx(rSeed); pTrck->GetParamLast()->SetTy(0.); pTrck->GetParamLast()->SetQp(0.); } } if (rootoutput) { finderCanvas->Update(); finderCanvas->Show(); char waitchar; cout << "press any key to continue." << endl; cin >> waitchar; delete finderCanvas; } if (fVerbose) { cout << endl; cout << "-------------------------------------------------------" << endl; cout << "-I- Ideal STT track finding -I-" << endl; cout << "Hits: " << nHits << endl; cout << "MCTracks: total " << nMCTracks << ", accepted " << nMCacc << ", reconstructable: " << nTracks << endl; cout << "SttHits not found : " << nNoSttHit << endl; cout << "SttPoints not found : " << nNoSttPoint << endl; cout << "MCTracks not found : " << nNoMCTrack << endl; cout << "SttTracks not found : " << nNoTrack << endl; cout << "-------------------------------------------------------" << endl; } else cout << "-I- CbmSttTrackFinderIdeal: all " << nMCTracks << ", acc. " << nMCacc << ", rec. " << nTracks << endl; return nTracks; } void CbmSttTrackFinderIdeal::GetTrackletCircular(Double_t firstX, Double_t firstY, Double_t firstR, Double_t secondX, Double_t secondY, Double_t secondR, Double_t thirdX, Double_t thirdY, Double_t thirdR, Double_t *circleRadii, Double_t *circleCentersX, Double_t *circleCentersY) const { Int_t trackletCounter = 0; Double_t small_limit = 0.0001; for (Int_t sign1 = -1; sign1 < 3; sign1 += 2) { for (Int_t sign2 = -1; sign2 < 3; sign2 += 2) { for (Int_t sign3 = -1; sign3 < 3; sign3 += 2) { // if three points are collinear, shift middle point by 50 microns if ((firstX - secondX == 0) && (secondX - thirdX == 0)) { secondX += 0.005; secondY += 0.005; } else if (!((fabs(firstX - secondX) < small_limit) || (fabs(secondX - thirdX) < small_limit))) { if (fabs((firstY - secondY) / (firstX - secondX) - (secondY - thirdY) / (secondX - thirdX)) < small_limit) { secondX += 0.005; secondY += 0.005; } } Double_t a = -2. * (firstX - secondX), b = -2. * (firstY - secondY), c = 2. * (sign1 * firstR - sign2 * secondR), d = ((firstX * firstX + firstY * firstY - firstR * firstR) - (secondX * secondX + secondY * secondY - secondR * secondR)), e = -2. * (firstX - thirdX), f = -2. * (firstY - thirdY), g = 2. * (sign1 * firstR - sign3 * thirdR), h = ((firstX * firstX + firstY * firstY - firstR * firstR) - (thirdX * thirdX + thirdY * thirdY - thirdR * thirdR)); Double_t A = -1. * (f * d - b * h) / (a * f - b * e), B = -1. * (f * c - b * g) / (a * f - b * e), C = (e * d - a * h) / (a * f - e * b), D = (e * c - a * g) / (a * f - e * b), I = B * B + D * D - 1., II = 2. * A * B - 2. * B * firstX + 2. * C * D - 2. * D * firstY + sign1 * 2. * firstR, III = A * A - 2. * A * firstX + firstX * firstX + C * C - 2. * C * firstY + firstY * firstY - firstR * firstR; if ((fabs(I) > small_limit) && ((II * II - 4. * I * III) > 0.)) { Double_t r = (-1. * II - sqrt(II * II - 4. * I * III)) / (2. * I), x = A + B * r, y = C + D * r; circleRadii[trackletCounter] = sqrt(r * r); circleCentersX[trackletCounter] = x; circleCentersY[trackletCounter] = y; trackletCounter++; } else { circleRadii[trackletCounter] = -1.; circleCentersX[trackletCounter] = 0.; circleCentersY[trackletCounter] = 0.; trackletCounter++; } } } } } void CbmSttTrackFinderIdeal::ZoomTrack(Double_t &dSeed, Double_t &phiSeed, Double_t &rSeed, CbmSttTrack *track) { track->SortHits(); if (track->GetNofHits() < 3) { cout << "-W- CbmSttTrackFinderHough::GetSeed2Circular: less than three hits found, not possible to fit!" << endl; return; } else { CbmSttHoughAccumulatorNew *accumulator = new CbmSttHoughAccumulatorNew(2, kTRUE, 150, dSeed - 5., dSeed + 5., 150, phiSeed - 0.5, phiSeed + 0.5, 100, rSeed - 200., rSeed + 200.); Int_t segmentLength1 = track->GetNofHits() / 3, segmentLength2 = track->GetNofHits() / 3, segmentLength3 = track->GetNofHits() - segmentLength2 - segmentLength1; Int_t segmentStart1 = 0, segmentStart2 = segmentLength1, segmentStart3 = segmentLength1 + segmentLength2; // select n random combinations for (Int_t teller = 0; teller < 100; teller++) { Int_t firstSample = gRandom->Integer(segmentLength1), secondSample = gRandom->Integer(segmentLength2), thirdSample = gRandom->Integer(segmentLength3); Int_t i = segmentStart1 + firstSample, ii = segmentStart2 + secondSample, iii = segmentStart3 + thirdSample; CbmSttHit *iFirstHit = GetHitFromCollections(track->GetHitIndex(i)), *iSecondHit = GetHitFromCollections(track->GetHitIndex(ii)), *iThirdHit = GetHitFromCollections(track->GetHitIndex(iii)); if ((!iFirstHit) || (!iSecondHit) || (!iThirdHit)) continue; Double_t firstX = iFirstHit->GetX(), firstY = iFirstHit->GetY(), firstR = iFirstHit->GetIsochrone(), secondX = iSecondHit->GetX(), secondY = iSecondHit->GetY(), secondR = iSecondHit->GetIsochrone(), thirdX = iThirdHit->GetX(), thirdY = iThirdHit->GetY(), thirdR = iThirdHit->GetIsochrone(); Double_t circleRadii[8], circleCentersX[8], circleCentersY[8]; GetTrackletCircular(firstX, firstY, firstR, secondX, secondY, secondR, thirdX, thirdY, thirdR, circleRadii, circleCentersX, circleCentersY); // loop over all tracklets between this and next for (Int_t trackletteller = 0; trackletteller < 8; trackletteller++) { if (circleRadii[trackletteller] > 0.) { Double_t phi = atan(circleCentersY[trackletteller] / circleCentersX[trackletteller]), dist = sqrt(circleCentersX[trackletteller] * circleCentersX[trackletteller] + circleCentersY[trackletteller] * circleCentersY[trackletteller]) - circleRadii[trackletteller]; if (circleCentersX[trackletteller] < 0.) { if (circleCentersY[trackletteller] < 0.) phi -= dPi; else phi += dPi; } accumulator->AddHit(dist, i, phi, ii, circleRadii[trackletteller], iii); } } } accumulator->Cluster(); dSeed = accumulator->GetHighestPeakX(); phiSeed = accumulator->GetHighestPeakY(); rSeed = accumulator->GetHighestPeakZ(); delete accumulator; } } void CbmSttTrackFinderIdeal::GetTrack(Double_t &dSeed, Double_t &phiSeed, Double_t &rSeed, Int_t mcTrackNo) { CbmMCTrack *mcTrack = (CbmMCTrack*) fMCTrackArray->At(mcTrackNo); // TODO: read field from container rSeed = sqrt(mcTrack->GetMomentum().X() * mcTrack->GetMomentum().X() + mcTrack->GetMomentum().Y() * mcTrack->GetMomentum().Y()) / 0.006; Double_t phiStart = atan(mcTrack->GetMomentum().Y() / mcTrack->GetMomentum().X()), xStart = mcTrack->GetStartVertex().X(), yStart = mcTrack->GetStartVertex().Y(); if (mcTrack->GetMomentum().X() < 0.) { if (mcTrack->GetMomentum().Y() < 0.) phiStart -= dPi; else phiStart += dPi; } Double_t sign = 1.; if (TDatabasePDG::Instance()->GetParticle(mcTrack->GetPdgCode()) != NULL) { if (TDatabasePDG::Instance()->GetParticle(mcTrack->GetPdgCode())->Charge() < 0.) { sign = -1.; } } Double_t xCircleCenter = xStart + sign * (rSeed * sin(phiStart)), yCircleCenter = yStart - sign * (rSeed * cos(phiStart)); cerr << xCircleCenter << " " << rSeed << " " << sign << " " << phiStart << endl; dSeed = sqrt(xCircleCenter * xCircleCenter + yCircleCenter * yCircleCenter) - rSeed; if (xCircleCenter == 0) { phiSeed = dPi / 2.; } else { phiSeed = atan(yCircleCenter / xCircleCenter); } if (xCircleCenter < 0.) { if (yCircleCenter < 0.) phiSeed -= dPi; else phiSeed += dPi; } } // ------------------------------------------------------------------------- Bool_t CbmSttTrackFinderIdeal::putStraw(Double_t xpos, Double_t ypos, Double_t radius) { Double_t pipeDiam = 0.42, tubeOuterDiam = 1.032, CoverThickness = 0.1, outerRadius = 42., innerRadius = 15.; if ((sqrt(xpos * xpos + ypos * ypos) < outerRadius - CoverThickness - (tubeOuterDiam / 2.)) && (sqrt(xpos * xpos + ypos * ypos) > innerRadius + CoverThickness + (tubeOuterDiam / 2.)) && (sqrt(ypos * ypos) > (pipeDiam + tubeOuterDiam / 2.))) { TArc myArc; myArc.SetFillStyle(0); myArc.SetLineColor(17); myArc.DrawArc(xpos, ypos, radius); } else { return false; } return true; } void CbmSttTrackFinderIdeal::plotAllStraws() { Double_t tubeOuterDiam = 1.032, sttCenterX = 0., sttCenterY = 0.; Int_t ringmax = -1; Bool_t started = kFALSE, goOn = kTRUE; Int_t ringteller = 18; while (goOn) { goOn = kFALSE; Double_t sqrt3 = sqrt(3.), radius = tubeOuterDiam / 2.; Double_t zpos = radius; // place first Double_t xpos = sttCenterX + ((ringteller) * 2 * radius), ypos = sttCenterY; for(Int_t i = 0; i < ringteller; i++) { xpos -= radius; ypos += sqrt3 * radius; if (putStraw(xpos, ypos, zpos)) goOn = kTRUE; } for(Int_t i = 0; i < ringteller; i++) { xpos -= 2 * radius; if (putStraw(xpos, ypos, zpos)) goOn = kTRUE; } for(Int_t i = 0; i < ringteller; i++) { xpos -= radius; ypos -= sqrt3 * radius; if (putStraw(xpos, ypos, zpos)) goOn = kTRUE; } for(Int_t i = 0; i < ringteller; i++) { xpos += radius; ypos -= sqrt3 * radius; if (putStraw(xpos, ypos, zpos)) goOn = kTRUE; } for(Int_t i = 0; i < ringteller; i++) { xpos += 2 * radius; if (putStraw(xpos, ypos, zpos)) goOn = kTRUE; } for(Int_t i = 0; i < ringteller; i++) { xpos += radius; ypos += sqrt3 * radius; if (putStraw(xpos, ypos, zpos)) goOn = kTRUE; } if (goOn) started = kTRUE; if (!started) goOn = kTRUE; ringteller++; if ((ringmax > 0) && (ringteller == ringmax)) goOn = kFALSE; } } CbmSttHit* CbmSttTrackFinderIdeal::GetHitFromCollections(Int_t hitCounter) { CbmSttHit *retval = NULL; Int_t relativeCounter = hitCounter; for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++) { Int_t size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast(); if (relativeCounter < size) { retval = (CbmSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter); break; } else { relativeCounter -= size; } } return retval; } CbmMCPoint* CbmSttTrackFinderIdeal::GetPointFromCollections(Int_t hitCounter) { CbmMCPoint *retval = NULL; Int_t relativeCounter = hitCounter; for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++) { Int_t size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast(); if (relativeCounter < size) { Int_t tmpHit = ((CbmSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter))->GetRefIndex(); retval = (CbmMCPoint*) ((TClonesArray *)fPointCollectionList.At(collectionCounter))->At(tmpHit); break; } else { relativeCounter -= size; } } return retval; } ClassImp(CbmSttTrackFinderIdeal)