// ------------------------------------------------------------------------- // ----- CbmSttTrackFinderIdeal source file ----- // ----- Created 28/03/06 by V. Friese ----- // ------------------------------------------------------------------------- // C++ includes #include #include // ROOT includes #include "TClonesArray.h" #include "TCanvas.h" #include "TH2F.h" #include "TArc.h" #include "TDatabasePDG.h" #include "TRandom.h" // CBM includes #include "FairMCPoint.h" #include "FairRootManager.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; else rootoutput = kFALSE; // stt1 } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- CbmSttTrackFinderIdeal::~CbmSttTrackFinderIdeal() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- void CbmSttTrackFinderIdeal::Init() { // Get and check FairRootManager FairRootManager* ioman = FairRootManager::Instance(); if (!ioman) { cout << "-E- CbmSttTrackFinderIdeal::Init: " << "RootManager not instantised!" << endl; return; } // Get MCTrack array // fMCTrackArray = (TClonesArray*) ioman->ActivateBranch("MCTrack"); fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack"); // EL 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; FairMCPoint* pMCpt = NULL; FairMCTrack* 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])[(Int_t)(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 = (FairMCTrack*) 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->GetNPoints(kTOF) << " 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); // CHECK!!! 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) { FairMCTrack *mcTrack = (FairMCTrack*) fMCTrackArray->At(mcTrackNo); TVector3 vectorMomentum; mcTrack->GetMomentum(vectorMomentum); TVector3 vectorTrack; mcTrack->GetStartVertex(vectorTrack); // TODO: read field from container rSeed = sqrt(vectorMomentum.X() * vectorMomentum.X() + vectorMomentum.Y() * vectorMomentum.Y()) / 0.006; Double_t phiStart = atan(vectorMomentum.Y() / vectorMomentum.X()), xStart = vectorTrack.X(), yStart = vectorTrack.Y(); if (vectorMomentum.X() < 0.) { if (vectorMomentum.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)); dSeed = sqrt(xCircleCenter * xCircleCenter + yCircleCenter * yCircleCenter) - rSeed; 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; } FairMCPoint* CbmSttTrackFinderIdeal::GetPointFromCollections(Int_t hitCounter) { FairMCPoint *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 = (FairMCPoint*) ((TClonesArray *)fPointCollectionList.At(collectionCounter))->At(tmpHit); break; } else { relativeCounter -= size; } } return retval; } ClassImp(CbmSttTrackFinderIdeal)