/* * PndFtsCellTrackletGenerator.cxx * * Created on: May 24, 2016 * Author: kibellus * * This class is used for the generation of track-candidates. * For this a cellular automaton is used. * The cellular automaton gets all FTS-hits as input and creates objects of the class PndTrackCand. * This algorithm add all hits to one track candidate which are neighbouring. * The method setHits shold be used to add the hits to the cellular automaton. * The method findTracks starts the cellular automaton. * The method getTracklets returns the finished track candidates. * After that the method reset shold used to clear the data structures. */ #include //ClassImp(PndFtsCellTrackletGenerator); PndFtsCellTrackletGenerator::PndFtsCellTrackletGenerator() { // TODO Auto-generated constructor stub } PndFtsCellTrackletGenerator::~PndFtsCellTrackletGenerator() { // TODO Auto-generated destructor stub } /* * Clears the data structure. * This is neccessary because the map fTracklets contains the track candidates of the step before. */ void PndFtsCellTrackletGenerator::reset() { fTracklets.clear(); } /* * This mehtod needs to be used before the cellular automaton can be start. */ void PndFtsCellTrackletGenerator::setHits(std::vector hits) { fHits = hits; } /* * This is the extern interface to start the cellular automaton. * In this method the cellular automaton will be started for all double layer. */ void PndFtsCellTrackletGenerator::findTracks() { std::cout << fHits.size() << " hits Found. " << std::endl; map> splittedHits = splitLayers(fHits); //first module for(int i=0;i<4;i++) findTracks(splittedHits,i); //second module for(int i=4;i<8;i++) findTracks(splittedHits,i); //third module for(int i=8;i<12;i++) findTracks(splittedHits,i); //fourth module for(int i=12;i<16;i++) findTracks(splittedHits,i); //fifth modules for(int i=16;i<20;i++) findTracks(splittedHits,i); //sixth module for(int i=20;i<24;i++) findTracks(splittedHits,i); } /* * Starts the cellular automaton for the layer from the parameter and write the result to the map fTracklets. */ void PndFtsCellTrackletGenerator::findTracks(map> splittedHits, Int_t layer){ map tracks = FindTracklets(splittedHits[layer]); typedef map::iterator iterator; for(iterator it = tracks.begin();it!=tracks.end();it++){ fTracklets[layer].push_back(it->second); } } /* * This method order all hits by the layer. * The return value is a map: layer -> list. * The method GetLayerID returns the id of the layer. * We need to merge each two layer to one because the tubes are located in double layers. */ map> PndFtsCellTrackletGenerator::splitLayers( vector hits) { map> map; for (int i = 0; i < hits.size(); i++) { PndFtsHit* hit = fHits[i]; map[(hit->GetLayerID() - 1) / 2].push_back(hit); } return map; } /* *The cellular automaton is implemented in this method. */ map PndFtsCellTrackletGenerator::FindTracklets( vector hits) { typedef std::map >::iterator it_type; typedef std::map::iterator it_type2; map states; //states of the hits map hitMap; //map ID -> hit map > neighbors = getNeighbors(hits); //the IDs of all neighboring hits //create state Map for (int i = 0; i < hits.size(); i++) { PndFtsHit* hit = hits[i]; states[hit->GetTubeID()] = hit->GetTubeID(); //init the state with the unique tube id hitMap[hit->GetTubeID()] = hit; } //run the cellular automaton Bool_t changes = kTRUE; //saves if a state changed while (changes) { changes = kFALSE; //for all hits for (it_type2 it = states.begin(); it != states.end(); it++) { //for all neighbors vector hitNeighbors = neighbors[it->second]; //get all neighbours of the hit for (int i = 0; i < hitNeighbors.size(); i++) { Int_t state1 = states[it->first]; //get the state of the considered hit Int_t state2 = states[hitNeighbors[i]]; //get the state of the neighbour if (state1 != state2) { //check if the cellular automaton need to change a state changes = kTRUE; //save that a state will be changed Int_t min = std::min(state1, state2); //get the minimum of both states states[it->first] = min; //set the state to the minimum states[hitNeighbors[i]] = min; } } } } //create the TrackCands map result; //get all hits with the same state and add them to one track candidate for (it_type2 it = states.begin(); it != states.end(); it++) { Int_t hitNumber = result[it->second].GetNHits()+1; result[it->second].AddHit(hitMap[it->first]->GetEntryNr(),hitNumber); } return result; } /* * Creates a map: tubeID of the hit -> list of tubeIDs of all neighbouring hits */ map > PndFtsCellTrackletGenerator::getNeighbors( vector hits) { map > neighbors; for (int i = 0; i < hits.size(); i++) { PndFtsHit* hit1 = hits[i]; transform(kTRUE,hit1); for (int j = i + 1; j < hits.size(); j++) { PndFtsHit* hit2 = hits[j]; transform(kTRUE,hit2); TVector3 p1; hit1->Position(p1); TVector3 p2; hit2->Position(p2); TVector3 diff = p1 - p2; Double_t norm = TMath::Sqrt(diff.X()*diff.X()+diff.Z()*diff.Z()); Double_t limit = 1.2; //if(hit1->GetSkewed()) limit=1.5; if (norm < limit) { //are neighbors? neighbors[hit1->GetTubeID()].push_back(hit2->GetTubeID()); neighbors[hit2->GetTubeID()].push_back(hit1->GetTubeID()); } transform(kTRUE,hit2); } transform(kFALSE,hit1); } return neighbors; } /* * creates a rotation matrix to rotate the coordinate system angle degrees around the y-axis. */ TMatrix PndFtsCellTrackletGenerator::getRotationMatrix(Double_t angle){ TMatrix a(3,3); a[0][0] = TMath::Cos(angle); a[0][1] = -TMath::Sin(angle); a[0][2] = 0; a[1][0] = TMath::Sin(angle); a[1][1] = TMath::Cos(angle); a[1][2] = 0; a[2][0] = 0; a[2][1] = 0; a[2][2] = 1; return a; } /* * Transforms the position of the hit to a new coordinate system (5° to the left or to the right) */ void PndFtsCellTrackletGenerator::transform(Bool_t transToNewSystem, PndFtsHit* hit){ Double_t angle = 0; if(((hit->GetLayerID() - 1) / 2)%4==1) angle = -5; if(((hit->GetLayerID() - 1) / 2)%4==2) angle = 5; if(angle==0) return; //no rotating for layer 1 and 4 if(!transToNewSystem) angle*=-1; TMatrix rotMat = getRotationMatrix(angle*TMath::DegToRad()); TMatrix hitMat(3,1); hitMat[0][0] = hit->GetX(); hitMat[1][0] = hit->GetY(); hitMat[2][0] = hit->GetZ(); TMatrix newPos = rotMat*hitMat; hit->SetXYZ(newPos[0][0],newPos[1][0],newPos[2][0]); }