//#include "tracking/PndOnlineSTTTripletFinder.h" #include "PndOnlineSttTripletFinder.h" #include "PndOnlineSttTriplet.h" //#include "PndDetectorList.h" #include "PndOnline.h" #include "PndSttTube.h" #include "PndSttHit.h" #include "PndSttPoint.h" #include "FairEventHeader.h" #include #include #include #include #include #include // ----- Public method Init (abstract in base class) -------------------- InitStatus PndOnlineSttTripletFinder::Init() { cout << "In PndOnlineSTTTripletFinder::Init()..." << endl; std::pair range1,range2,range3; // set up some ranges corresponding to certain // straw numbers are ranges excluding the border straws, i.e. first..first+1...second-1..second range1.first = 104; range1.second = 215; range2.first = 714; range2.second = 855; range3.first = 2956; range3.second = 3175; pivot_tube_ranges.push_back(range1); pivot_tube_ranges.push_back(range2); pivot_tube_ranges.push_back(range3); pair skew1, skew2, skew3, skew4; skew1.first = 1000; skew1.second = 1395; skew2.first = 1394; skew2.second = 1813; skew3.first = 1812; skew3.second = 2267; skew4.first = 2266; skew4.second = 2745; skew_tube_ranges.push_back(skew1); skew_tube_ranges.push_back(skew2); skew_tube_ranges.push_back(skew3); skew_tube_ranges.push_back(skew4); for (int i = 1; i <= tube_array->GetEntriesFast(); i++) { PndSttTube* leftTube = (PndSttTube*)(tube_array->At(i)); if (leftTube == 0) continue; for (int j = 1; j <= tube_array->GetEntriesFast(); j++) { if (i == j) continue; PndSttTube* rightTube = (PndSttTube*)(tube_array->At(j)); if (rightTube == 0) continue; if (GetTubeDist3D(i, j, tube_array) < 1.1) { fTubeNeighbors[i].push_back(j); } } } verbosity = 0; return kSUCCESS; } // ----- Public method Exec -------------------------------------------- void PndOnlineSttTripletFinder::Exec(Option_t* opt) { cout << "In PndOnlineSTTTripletFinder::Exec()..." << endl; if(online_manager == NULL) { return; } TStopwatch timer; timer.Start(); std::vector tripletstorage; std::vector skewlethalfstorage; std::vector skewletstorage; //find triplets for (vector< pair >::iterator it=pivot_tube_ranges.begin(); it!=pivot_tube_ranges.end(); ++it) { if (verbosity > 2) { cout << " checking tubes " << it->first << " to " << it->second << "..." << endl; } TObjArray* triplets = (TObjArray*)FindTriplets(online_manager->GetHitChannelMap(PndOnline::kSTT), it->first, it->second); tripletstorage.push_back(triplets); } //cout << "Tripletstorage Size: " << tripletstorage.size() << endl; //find skewlet halves for (vector< pair >::iterator it=skew_tube_ranges.begin(); it!=skew_tube_ranges.end(); ++it) { //vector< pair >::iterator nextit = it; //if (++nextit == skew_tube_ranges.end()) break; cout << " checking skewed tubes " << it->first << " to " << it->second << "..." << endl; //cout << " and tubes " << nextit->first << " to " << nextit->second << "..." << endl; TObjArray* skewlets = FindSkewletHalf(online_manager->GetHitChannelMap(PndOnline::kSTT), it->first, it->second); skewlethalfstorage.push_back(skewlets); } //skewlets will be added later for (int i = 0; i < skewlethalfstorage.size(); i++) { if (i + 2 > skewlethalfstorage.size()) break; cout << "Combining Skew Halves " << i << " and " << i+1 << endl; TObjArray* skewlets = CombineSkewlets(skewlethalfstorage.at(i), skewlethalfstorage.at(i+1)); skewletstorage.push_back(skewlets); } //MakeTripletTracks(tripletstorage); //MakeSkewletTracks(skewletstorage); std::vector > triplettrackcands_1_3; MakeTripletTracks(tripletstorage, 0, 2, triplettrackcands_1_3); TrackPostProcessingHitProximity(triplettrackcands_1_3, online_manager->GetHitChannelMap(PndOnline::kSTT), skewletstorage); std::vector > triplettrackcands_2_3; MakeTripletTracks(tripletstorage, 1, 2, triplettrackcands_2_3); TrackPostProcessingHitProximity(triplettrackcands_2_3, online_manager->GetHitChannelMap(PndOnline::kSTT), skewletstorage); std::vector > triplettrackcands_1_2; MakeTripletTracks(tripletstorage, 0, 1, triplettrackcands_1_2); TrackPostProcessingHitProximity(triplettrackcands_1_2, online_manager->GetHitChannelMap(PndOnline::kSTT), skewletstorage); for (int i = 0; i < tripletstorage.size(); i++) { delete tripletstorage.at(i); } for (int i = 0; i < skewlethalfstorage.size(); i++) { delete skewlethalfstorage.at(i); } for (int i = 0; i < skewletstorage.size(); i++) { delete skewletstorage.at(i); } timer.Stop(); last_run_time = timer.CpuTime(); } // search over given range of tubes, see if any have hits, then look for connected groups of hits TObjArray* PndOnlineSttTripletFinder::FindTriplets(vector< list > hit_channel_map, Int_t starttube, Int_t stoptube) { if (verbosity > 2) { cout << "PndOnlineSTTTripletFinder::FindTriplets() ... " << endl << " number of channels in hit_channel_map... " << hit_channel_map.size() << endl; } TObjArray* tripletlist = new TObjArray; tripletlist->SetOwner(kTRUE); for (Int_t tubeID = starttube+1; tubeID < stoptube; tubeID++) { //cout << " checking tube " << tubeID << "..."; if (hit_channel_map[tubeID].size() == 0) { //cout << " empty." << endl; continue; } //cout << "FindTriplets: Hits in Pivot Tube " << tubeID << endl; //process all pivot hits of the current tube in the buffer for(list::iterator iter=hit_channel_map[tubeID].begin(); iter!=hit_channel_map[tubeID].end(); iter++) { ProcessPivotTubeHit(tubeID, hit_channel_map, iter, tripletlist); } } return tripletlist; } void PndOnlineSttTripletFinder::ProcessPivotTubeHit(Int_t pivotID, vector< list > hit_channel_map, list::iterator pivotIter, TObjArray* tripletlist) { PndOnlineSttTriplet* mytriplet = new PndOnlineSttTriplet; PndSttHit* pivotHit = (PndSttHit*)((*pivotIter)->BaseHit()); PndSttTube* pivotTube = (PndSttTube*)(tube_array->At(pivotHit->GetTubeID())); mytriplet->SetPivotHit(pivotHit, pivotTube); Double_t pivotTime = pivotHit->GetTimeStamp(); // imbue MC event header information for comparison Double_t mcEventTime = 0; Int_t mcEventID = 0; Int_t mcTrackID = 0; Int_t mcSttPointID = 0; ExtractMCLinkInfo(pivotHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); if (verbosity > 3) { cout << "Pivot Event Time: " << mcEventTime << endl; cout << "Pivot Event ID: " << mcEventID << endl; cout << "Pivot Track ID: " << mcTrackID << endl; cout << "Pivot Point ID: " << mcSttPointID << endl; } // we loop over all other straws to add them to the triplet // later we will iterate only over known neighbors (MUCH faster, but needs neighbor map) //for(int rightID = 0; rightID < hit_channel_map.size(); rightID++) { for (int i = 0; i < fTubeNeighbors[pivotID].size(); i++) { int rightID = fTubeNeighbors[pivotID].at(i); if(hit_channel_map[rightID].size() == 0) { continue; } if(pivotID == rightID) { continue; } if (GetTubeDist(pivotID, rightID, tube_array) < 1.5) { // iterate over all hits in the nearby tube // later we will iterate only over "not expired" hits (MUCH faster, but needs logic to keep track of hit expiration) for(list::iterator iter=hit_channel_map[rightID].begin(); iter!=hit_channel_map[rightID].end(); iter++) { PndSttHit* rightHit = (PndSttHit*)((*iter)->BaseHit()); PndSttTube* rightTube = (PndSttTube*)(tube_array->At(rightHit->GetTubeID())); Double_t rightTime = rightHit->GetTimeStamp(); // we check only if the right hit is within the pivot time window (+-drift time) // later we will check if all hits are within a common drift time window (more precise, but results in triplet candidate splitting) if ( (pivotTime + PndOnlineSTT::DriftTime > rightTime) && (pivotTime - PndOnlineSTT::DriftTime < rightTime) ) { if (verbosity > 3) { cout << "Neighbor of tube " << pivotID << " found in tube " << rightID << " within time window: " << pivotTime - PndOnlineSTT::DriftTime << ", " << pivotTime + PndOnlineSTT::DriftTime << endl; cout << "Neighbor Timestamp: " << rightTime << endl; cout << "Acceptable Timestamps: " << mytriplet->GetMinAcceptWindow() << ", " << mytriplet->GetMaxAcceptWindow() << endl; cout << "Triplet T0 data: " << mytriplet->GetMinT0() << " " << mytriplet->GetMaxT0() << endl; } if (mytriplet->HitTooOld(rightHit)) { if (verbosity > 1) { cout << "mytriplet detected too old hit and will discard it. This should not happen." << endl; } } if (mytriplet->HitTooRecent(rightHit)) { // in the future, the triplet should be split into two valid candidates if (verbosity > 2) { cout << "mytriplet detected too recent hit and will discard it." << endl; } } mytriplet->AddNeighbor(rightHit, rightTube); //imbue MC information ExtractMCLinkInfo(rightHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); mytriplet->AddMCIndex(mcEventID); mytriplet->AddTrackIndex(mcTrackID); mytriplet->MCEventEntryNumber = mytriplet->MajorMCIndex(); mytriplet->MCTrackEntryNumber = mytriplet->MajorTrackIndex(); if (verbosity > 3) { cout << "Right Event Time: " << mcEventTime << endl; cout << "Right Event ID: " << mcEventID << endl; cout << "Right Track ID: " << mcTrackID << endl; cout << "Right Point ID: " << mcSttPointID << endl; cout << "Major MC ID: " << mytriplet->MajorMCIndex() << endl; cout << "MC ID Majority: " << mytriplet->MCIndexMajority() << endl; cout << "Major Track ID: " << mytriplet->MajorTrackIndex() << endl; cout << "Track ID Majority: " << mytriplet->TrackIndexMajority() << endl; cout << "Triplet now has " << mytriplet->GetHitCount() << " hits." << endl; } } } } } if (TripletGood(mytriplet)) { //tripletlist must be the owner of the list objects tripletlist->AddLast(mytriplet); // hack: store a copy of the triplet as onlinetrack with goodness 0 to display as onlinetrackvertex PndOnlineTrack* new_track = new PndOnlineTrack(); new_track->SetVertex(mytriplet->GetCMS(tube_array)); new_track->SetT0(mytriplet->GetMinT0()); new_track->SetGoodness(0); new_track->SetDrawMode(PndOnlineTrack::kVertexAsPoint); if (verbosity > 3) { cout << "Adding Triplet Track to OnlineManager, Timestamp: " << new_track->T0() << endl; } online_manager->AddTrack(new_track); } } Bool_t PndOnlineSttTripletFinder::TripletGood(const PndOnlineSttTriplet* const mytriplet) { bool triplet_ok = true; // to be implemented: check cluttering // check hit count validity: a good triplet should have between 3 and 5 hits if ( (mytriplet->GetHitCount() < 3) || (mytriplet->GetHitCount() > 6) ) { triplet_ok = false; } return triplet_ok; } // generate tracks from pairs of triplets void PndOnlineSttTripletFinder::MakeTripletTracks(vector& tripletstorage) { map > maxdeflection; std::vector > tripletcands; maxdeflection[0][1]=10; maxdeflection[0][2]=45; maxdeflection[1][2]=25; for (int i = 0; i < tripletstorage.size(); i++) { TObjArray*innertriplets = tripletstorage.at(i); if (innertriplets == NULL) { if (verbosity > 0) { cout << "TObjArray* to Triplets is 0" << endl; } continue; } for (int j = i+1; j < tripletstorage.size(); j++) { TObjArray*outertriplets = tripletstorage.at(j); if (outertriplets == NULL) continue; //calculate tracks for inner and outer triplet combinations CombineLayers(innertriplets, outertriplets, maxdeflection[i][j], tripletcands); //TrackPostProcessingTripletLayer(tripletcands, tripletlayer); //TrackPostProcessingHitProximity(tripletcands, online_manager->GetHitChannelMap(PndOnline::kSTT)); //TrackPostProcessingMVDProximity(tripletcands, hitmap); } } } void PndOnlineSttTripletFinder::MakeTripletTracks(vector& tripletstorage, int i, int j, std::vector >& tripletcands) { map > maxdeflection; //std::vector > tripletcands; maxdeflection[0][1]=30; maxdeflection[0][2]=90; maxdeflection[1][2]=35; TObjArray*innertriplets = tripletstorage.at(i); TObjArray*outertriplets = tripletstorage.at(j); if ((innertriplets != NULL) && (outertriplets != NULL)) { //calculate tracks for inner and outer triplet combinations CombineLayers(innertriplets, outertriplets, maxdeflection[i][j], tripletcands); //TrackPostProcessingTripletLayer(tripletcands, tripletlayer); //TrackPostProcessingHitProximity(tripletcands, online_manager->GetHitChannelMap(PndOnline::kSTT)); //TrackPostProcessingMVDProximity(tripletcands, hitmap); } } void PndOnlineSttTripletFinder::CombineLayers(TObjArray* innertriplets, TObjArray* outertriplets, double maxphi, std::vector >& tripletcands) { for (int k = 0; k < innertriplets->GetEntriesFast(); k++) { PndOnlineSttTriplet* innertriplet = (PndOnlineSttTriplet *)(innertriplets->At(k)); for (int l = 0; l < outertriplets->GetEntriesFast(); l++) { PndOnlineSttTriplet* outertriplet = (PndOnlineSttTriplet *)(outertriplets->At(l)); if ((innertriplet->GetMinT0() > outertriplet->GetMaxT0()) || (innertriplet->GetMaxT0() < outertriplet->GetMinT0()) ) { //cout << "Triplet Times too far apart. Skipping Track. Inner Triplet: " << innertriplet->GetMinT0() << " - " << innertriplet->GetMaxT0() << " Outer Triplet: " << outertriplet->GetMinT0() << " - " << outertriplet->GetMaxT0() << endl; continue; } TVector3 innercms = innertriplet->GetCMS(tube_array); TVector3 outercms = outertriplet->GetCMS(tube_array); //cout << " Inner Phi: " << i << " " << innercms->Phi()*TMath::RadToDeg() << endl; //cout << " Outer Phi: " << j << " " << outercms->Phi()*TMath::RadToDeg() << endl; //cout << "Max Deflection: " << (maxdeflection[i][j]) << endl; if (TMath::Abs(innercms.DeltaPhi(outercms)) > (maxphi)*TMath::DegToRad()) { continue; } tripletcands.push_back(std::make_pair(innertriplet, outertriplet)); } } } void PndOnlineSttTripletFinder::TrackPostProcessingHitProximity(std::vector >& tripletcands, vector< list > hit_channel_map, std::vector skewletstorage) { for (int i = 0; i < tripletcands.size(); i++) { PndOnlineSttTriplet* innertriplet = tripletcands.at(i).first; PndOnlineSttTriplet* outertriplet = tripletcands.at(i).second; Double_t windowstart = innertriplet->GetMinAcceptWindow(); if (verbosity > 2) { cout << "Checking Triplet Track Candidate " << i << endl; cout << innertriplet->GetCMS(tube_array).X() << ", " << innertriplet->GetCMS(tube_array).Y() << ", " << innertriplet->GetCMS(tube_array).Z() << ", " << endl; cout << outertriplet->GetCMS(tube_array).X() << ", " << outertriplet->GetCMS(tube_array).Y() << ", " << outertriplet->GetCMS(tube_array).Z() << ", " << endl; } if (outertriplet->GetMinAcceptWindow() > windowstart) { windowstart = outertriplet->GetMinAcceptWindow(); } Double_t windowend = innertriplet->GetMaxAcceptWindow(); if (outertriplet->GetMaxAcceptWindow() < windowend) { windowend = outertriplet->GetMaxAcceptWindow(); } if (verbosity > 2) { cout << "Time Window: " << windowstart << ", " << windowend << std::endl; } Double_t originx = 0; Double_t originy = 0; Double_t radius = 0; Int_t clockwise = -1; CalculateCircle(innertriplet, outertriplet, originx, originy, radius, clockwise); TVector2 circleorigin = TVector2(originx, originy); Double_t orthophi = circleorigin.Phi(); Double_t trackphi = orthophi + 0.5*TMath::Pi()*clockwise; if (trackphi < 0) { trackphi += 2*TMath::Pi(); } if (trackphi > 2*TMath::Pi()) { trackphi -= 2*TMath::Pi(); } Int_t innerrow = fStrawMap.GetRow(innertriplet->GetPivotStraw()); Int_t outerrow = fStrawMap.GetRow(outertriplet->GetPivotStraw()); if (verbosity > 2) { cout << "Pivot Row Range: " << innerrow << ", " << outerrow << std::endl; } Int_t innersector = fStrawMap.GetSector(innertriplet->GetPivotStraw()); Int_t outersector = fStrawMap.GetSector(outertriplet->GetPivotStraw()); if (verbosity > 2) { cout << "Pivot Sector Range: " << innersector << ", " << outersector << std::endl; } Int_t minrow, maxrow; if (innerrow < outerrow) { minrow = innerrow - 1; maxrow = outerrow + 2; } else { minrow = outerrow - 1; maxrow = innerrow + 2; } Int_t minsector, maxsector; if (innersector < outersector) { minsector = innersector; maxsector = outersector; } else { minsector = outersector; maxsector = innersector; } /* for (all layers to check) { count checked layers for (hits in hit channel map of layer) { CalculateCircleProximity(originx, originy, radius, clockwise, stthit); cout checked hits } } */ Int_t trackhits = 0; PndOnlineSttTriplet bookkeeper; //cout << "Hit Channel Map Size: " << hit_channel_map.size() << std::endl; for(int isector = minsector; isector <= maxsector; isector++) { for (int irow = minrow; irow <= maxrow; irow++) { for (int istraw = 0; istraw < fStrawMap.GetStrawRow(isector, irow).size(); istraw++) { int rightID = fStrawMap.GetStrawRow(isector, irow).at(istraw); //for(int rightID = 0; rightID < hit_channel_map.size(); rightID++) { //cout << "Channel Size: " << hit_channel_map[rightID].size() << std::endl; if ( (fStrawMap.GetRow(rightID) < minrow) || (fStrawMap.GetRow(rightID) > maxrow) ) { continue; } if(fStrawMap.IsSkewedStraw(rightID)) { continue; } if(hit_channel_map[rightID].size() == 0) { continue; } if (HitInRange(originx, originy, radius, clockwise, rightID, 0.5)) { // iterate over all hits in the nearby tube // later we will iterate only over "not expired" hits (MUCH faster, but needs logic to keep track of hit expiration) Int_t tubehits = 0; for(list::iterator iter=hit_channel_map[rightID].begin(); iter!=hit_channel_map[rightID].end(); iter++) { PndSttHit* rightHit = (PndSttHit*)((*iter)->BaseHit()); PndSttTube* rightTube = (PndSttTube*)(tube_array->At(rightHit->GetTubeID())); Double_t rightTime = rightHit->GetTimeStamp(); // we check only if the right hit is within the pivot time window (+-drift time) // later we will check if all hits are within a common drift time window (more precise, but results in triplet candidate splitting) if ( (windowend > rightTime) && (windowstart < rightTime) ) { //std::cout << "Hit near track found in tube " << rightID << " within time window: " << windowstart << ", " << windowend << std::endl; ++tubehits; //cout << "Found " << tubehits << " tube hits." << endl; if (!bookkeeper.IsActive()) { bookkeeper.SetPivotHit(rightHit, rightTube); } else { bookkeeper.AddNeighbor(rightHit, rightTube, true); } Double_t mcEventTime; Int_t mcEventID; Int_t mcTrackID; Int_t mcSttPointID; ExtractMCLinkInfo(rightHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); bookkeeper.AddMCIndex(mcEventID); bookkeeper.AddTrackIndex(mcTrackID); } } if (tubehits > 0) { //cout << "Nearby tube had " << tubehits << " hits. Adding hit to track." << endl; ++trackhits; } } }}} //check skewlet proximity for (int j = 0; j < skewletstorage.size(); j++) { TObjArray* skewlets = skewletstorage.at(j); if (skewlets == NULL) { if (verbosity > 0) { cout << "TObjArray* to Skewlets is 0" << endl; } continue; } Int_t skewhits = 0; for (int k = 0; k < skewlets->GetEntriesFast(); k++) { PndOnlineSttTriplet* skewlet = (PndOnlineSttTriplet*)(skewlets->At(k)); if (HitInRange(originx, originy, radius, clockwise, TVector2(skewlet->GetPoca().X(), skewlet->GetPoca().Y()), 0.5)) { skewhits++; } } if (skewhits > 0) { //todo: add mc indexing of skewed straws to bookkeeper ++trackhits; } } if (verbosity > 2) { cout << "Major MC ID: " << bookkeeper.MajorMCIndex() << endl; cout << "MC ID Majority: " << bookkeeper.MCIndexMajority() << endl; cout << "Major Track ID: " << bookkeeper.MajorTrackIndex() << endl; cout << "Track ID Majority: " << bookkeeper.TrackIndexMajority() << endl; cout << "Triplet now has " << bookkeeper.GetHitCount() << " hits." << endl; cout << "This track has passed " << maxrow - minrow << " layers" << endl; cout << "This track has " << trackhits << " hits." << endl; } map minimumhits; //std::vector > tripletcands; minimumhits[8]=7; minimumhits[14]=12; minimumhits[19]=12; if (trackhits < minimumhits[maxrow - minrow]) continue; TVector3 innercms = innertriplet->GetCMS(tube_array); TVector3 outercms = outertriplet->GetCMS(tube_array); // hack: add track with goodness 10 for track display // vertex and momentum are (ab)used as the two points on the circle PndOnlineTrack* new_track = new PndOnlineTrack(); new_track->SetVertex(innercms); new_track->SetMomentum(outercms); new_track->SetGoodness(10); new_track->SetDrawMode(PndOnlineTrack::kVertexMomentumAsCircleHack); if (innertriplet->GetMinT0() < outertriplet->GetMinT0()) { new_track->SetT0(outertriplet->GetMinT0()); } else { new_track->SetT0(innertriplet->GetMinT0()); } //imbue MC information //new_track->SetT0Min(innertriplet->MCEventTime); //new_track->SetT0Max(outertriplet->MCEventTime); new_track->SetMCEntryNumber(bookkeeper.MajorMCIndex()); new_track->SetMCEntryMajority(bookkeeper.MCIndexMajority()); new_track->SetMCTrackNumber(bookkeeper.MajorTrackIndex()); new_track->SetMCTrackMajority(bookkeeper.TrackIndexMajority()); new_track->SetHitCount(trackhits); //new_track->SetMCEventTime(mcEventTime); new_track->SetT0Min(bookkeeper.GetMinT0()); new_track->SetT0Max(bookkeeper.GetMaxT0()); new_track->SetCircleOrigin(circleorigin); new_track->SetCircleRadius(radius); new_track->SetCircleClockwise(clockwise); new_track->SetCirclePhi(trackphi); //new_track->SetMCEventTime(innertriplet->MCEventTime); if (verbosity > 2) { cout << "Adding Found Track to OnlineManager, Timestamp: " << new_track->T0() << endl; } online_manager->AddTrack(new_track); } } //{ // for (int i = 0; i < tripletcands.size(); i++) { // PndOnlineSttTriplet* innertriplet = tripletcands.at(i).first; // PndOnlineSttTriplet* outertriplet = tripletcands.at(i).second; // Double_t windowstart = innertriplet->GetMinAcceptWindow(); // cout << "Checking Triplet Track Candidate " << i << endl; // cout << innertriplet->GetCMS(tube_array).X() << ", " << innertriplet->GetCMS(tube_array).Y() << ", " << innertriplet->GetCMS(tube_array).Z() << ", " << endl; // cout << outertriplet->GetCMS(tube_array).X() << ", " << outertriplet->GetCMS(tube_array).Y() << ", " << outertriplet->GetCMS(tube_array).Z() << ", " << endl; // if (outertriplet->GetMinAcceptWindow() > windowstart) { // windowstart = outertriplet->GetMinAcceptWindow(); // } // Double_t windowend = innertriplet->GetMaxAcceptWindow(); // if (outertriplet->GetMaxAcceptWindow() < windowend) { // windowend = outertriplet->GetMaxAcceptWindow(); // } // cout << "Time Window: " << windowstart << ", " << windowend << std::endl; // Double_t originx = 0; // Double_t originy = 0; // Double_t radius = 0; // Int_t clockwise = -1; // CalculateCircle(innertriplet, outertriplet, originx, originy, radius, clockwise); // Int_t innerrow = fStrawMap.GetRow(innertriplet->GetPivotStraw()); // Int_t outerrow = fStrawMap.GetRow(outertriplet->GetPivotStraw()); // cout << "Pivot Row Range: " << innerrow << ", " << outerrow << std::endl; // Int_t innersector = fStrawMap.GetSector(innertriplet->GetPivotStraw()); // Int_t outersector = fStrawMap.GetSector(outertriplet->GetPivotStraw()); // cout << "Pivot Sector Range: " << innersector << ", " << outersector << std::endl; // Int_t minrow, maxrow; // if (innerrow < outerrow) { // minrow = innerrow - 1; // maxrow = outerrow + 2; // } else { // minrow = outerrow - 1; // maxrow = innerrow + 2; // } // Int_t minsector, maxsector; // if (innersector < outersector) { // minsector = innersector; // maxsector = outersector; // } else { // minsector = outersector; // maxsector = innersector; // } // // for (Int_t mySector = minsector; mySector < maxsector; mySector++) { // for (Int_t myRow = minrow; myRow < maxrow; myRow++) { // // } // } // // /* // for (all layers to check) { // count checked layers // for (hits in hit channel map of layer) { // CalculateCircleProximity(originx, originy, radius, clockwise, stthit); // cout checked hits // } // } // */ // Int_t trackhits = 0; // //cout << "Hit Channel Map Size: " << hit_channel_map.size() << std::endl; // for(int rightID = 0; rightID < hit_channel_map.size(); rightID++) { // //cout << "Channel Size: " << hit_channel_map[rightID].size() << std::endl; // if(hit_channel_map[rightID].size() == 0) { // continue; // } // if (HitInRange(originx, originy, radius, clockwise, rightID, 0.5)) { // // iterate over all hits in the nearby tube // // later we will iterate only over "not expired" hits (MUCH faster, but needs logic to keep track of hit expiration) // Int_t tubehits = 0; // for(list::iterator iter=hit_channel_map[rightID].begin(); iter!=hit_channel_map[rightID].end(); iter++) { // PndSttHit* rightHit = (PndSttHit*)((*iter)->BaseHit()); // PndSttTube* rightTube = (PndSttTube*)(tube_array->At(rightHit->GetTubeID())); // Double_t rightTime = rightHit->GetTimeStamp(); // // we check only if the right hit is within the pivot time window (+-drift time) // // later we will check if all hits are within a common drift time window (more precise, but results in triplet candidate splitting) // if ( (windowend > rightTime) && (windowstart < rightTime) ) { // //std::cout << "Hit near track found in tube " << rightID << " within time window: " << windowstart << ", " << windowend << std::endl; // ++tubehits; // //cout << "Found " << tubehits << " tube hits." << endl; // } // } // if (tubehits > 0) { // //cout << "Nearby tube had " << tubehits << " hits. Adding hit to track." << endl; // ++trackhits; // } // } // } // cout << "This track has " << trackhits << " hits." << endl; // } //} ////////////////////////////////////////////////////////////////////// //////////////////// SKEWED STRAWS ////////////////////////////////////////////////////////////////////// TObjArray* PndOnlineSttTripletFinder::FindSkewletHalf(vector< list > hit_channel_map, Int_t starttube, Int_t stoptube) { if (verbosity > 2) { cout << "PndOnlineSTTTripletFinder::FindSkewletHalf() ... " << endl << " number of channels in hit_channel_map... " << hit_channel_map.size() << endl; } TObjArray* triplethalflist = new TObjArray; triplethalflist->SetOwner(kTRUE); for (Int_t tubeID = starttube+1; tubeID < stoptube; tubeID++) { //cout << " checking tube " << tubeID << "..."; if (hit_channel_map[tubeID].size() == 0) { //cout << " empty." << endl; continue; } //cout << "FindTriplets: Hits in Pivot Tube " << tubeID << endl; //process all pivot hits of the current tube in the buffer for(list::iterator iter=hit_channel_map[tubeID].begin(); iter!=hit_channel_map[tubeID].end(); iter++) { ProcessSkewPivotTubeHit(tubeID, hit_channel_map, iter, triplethalflist); } } return triplethalflist; } void PndOnlineSttTripletFinder::ProcessSkewPivotTubeHit(Int_t pivotID, vector< list > hit_channel_map, list::iterator pivotIter, TObjArray* tripletlist) { PndOnlineSttTriplet* mytriplet = new PndOnlineSttTriplet; PndSttHit* pivotHit = (PndSttHit*)((*pivotIter)->BaseHit()); PndSttTube* pivotTube = (PndSttTube*)(tube_array->At(pivotHit->GetTubeID())); mytriplet->SetPivotHit(pivotHit, pivotTube); Double_t pivotTime = pivotHit->GetTimeStamp(); // imbue MC event header information for comparison Double_t mcEventTime = 0; Int_t mcEventID = 0; Int_t mcTrackID = 0; Int_t mcSttPointID = 0; ExtractMCLinkInfo(pivotHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); if (verbosity > 3) { cout << "Pivot Event Time: " << mcEventTime << endl; cout << "Pivot Event ID: " << mcEventID << endl; cout << "Pivot Track ID: " << mcTrackID << endl; cout << "Pivot Point ID: " << mcSttPointID << endl; } // we loop over all other straws to add them to the triplet // later we will iterate only over known neighbors (MUCH faster, but needs neighbor map) //for(int rightID = 0; rightID < hit_channel_map.size(); rightID++) { for (int i = 0; i < fTubeNeighbors[pivotID].size(); i++) { int rightID = fTubeNeighbors[pivotID].at(i); if(hit_channel_map[rightID].size() == 0) { continue; } if(pivotID == rightID) { continue; } if (GetTubeDist3D(pivotID, rightID, tube_array) < 1.1) { // iterate over all hits in the nearby tube // later we will iterate only over "not expired" hits (MUCH faster, but needs logic to keep track of hit expiration) for(list::iterator iter=hit_channel_map[rightID].begin(); iter!=hit_channel_map[rightID].end(); iter++) { PndSttHit* rightHit = (PndSttHit*)((*iter)->BaseHit()); PndSttTube* rightTube = (PndSttTube*)(tube_array->At(rightHit->GetTubeID())); Double_t rightTime = rightHit->GetTimeStamp(); Double_t wireanglex = TMath::Abs(pivotTube->GetWireDirection().X() - rightTube->GetWireDirection().X()); Double_t wireangley = TMath::Abs(pivotTube->GetWireDirection().Y() - rightTube->GetWireDirection().Y()); //match only parallel tubes if ((wireanglex > 0.005) || (wireangley > 0.005)) { continue; } // we check only if the right hit is within the pivot time window (+-drift time) // later we will check if all hits are within a common drift time window (more precise, but results in triplet candidate splitting) if ( (pivotTime + PndOnlineSTT::DriftTime > rightTime) && (pivotTime - PndOnlineSTT::DriftTime < rightTime) ) { if (verbosity > 3) { cout << "Neighbor of tube " << pivotID << " found in tube " << rightID << " within time window: " << pivotTime - PndOnlineSTT::DriftTime << ", " << pivotTime + PndOnlineSTT::DriftTime << endl; cout << "Neighbor Timestamp: " << rightTime << endl; cout << "Acceptable Timestamps: " << mytriplet->GetMinAcceptWindow() << ", " << mytriplet->GetMaxAcceptWindow() << endl; cout << "Triplet T0 data: " << mytriplet->GetMinT0() << " " << mytriplet->GetMaxT0() << endl; } if (mytriplet->HitTooOld(rightHit)) { if (verbosity > 1) { cout << "mytriplet detected too old hit and will discard it. This should not happen." << endl; } } if (mytriplet->HitTooRecent(rightHit)) { // in the future, the triplet should be split into two valid candidates if (verbosity > 2) { cout << "mytriplet detected too recent hit and will discard it." << endl; } } mytriplet->AddNeighbor(rightHit, rightTube); //imbue MC information ExtractMCLinkInfo(rightHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); mytriplet->AddMCIndex(mcEventID); mytriplet->AddTrackIndex(mcTrackID); mytriplet->MCEventEntryNumber = mytriplet->MajorMCIndex(); mytriplet->MCTrackEntryNumber = mytriplet->MajorTrackIndex(); if (verbosity > 3) { cout << "Right Event Time: " << mcEventTime << endl; cout << "Right Event ID: " << mcEventID << endl; cout << "Right Track ID: " << mcTrackID << endl; cout << "Right Point ID: " << mcSttPointID << endl; cout << "Major MC ID: " << mytriplet->MajorMCIndex() << endl; cout << "MC ID Majority: " << mytriplet->MCIndexMajority() << endl; cout << "Major Track ID: " << mytriplet->MajorTrackIndex() << endl; cout << "Track ID Majority: " << mytriplet->TrackIndexMajority() << endl; cout << "Triplet now has " << mytriplet->GetHitCount() << " hits." << endl; } } } } }/* if (TripletGood(mytriplet)) { //tripletlist must be the owner of the list objects tripletlist->AddLast(mytriplet); // hack: store a copy of the triplet as onlinetrack with goodness 0 to display as onlinetrackvertex PndOnlineTrack* new_track = new PndOnlineTrack(); new_track->SetVertex(mytriplet->GetCMS3D(tube_array)); new_track->SetT0(mytriplet->GetMinT0()); new_track->SetGoodness(0); new_track->SetDrawMode(PndOnlineTrack::kVertexAsPoint); if (verbosity > 3) { cout << "Adding Triplet Track to OnlineManager, Timestamp: " << new_track->T0() << endl; } online_manager->AddTrack(new_track); }*/ tripletlist->AddLast(mytriplet); } TObjArray* PndOnlineSttTripletFinder::CombineSkewlets(TObjArray* innerlist, TObjArray* outerlist) { TObjArray* fulltripletlist = new TObjArray; fulltripletlist->SetOwner(kTRUE); //cout << "Combine Skewlets. Inner list: " << innerlist->GetEntriesFast() << ", Outer list: " << outerlist->GetEntriesFast() << endl; for (int i = 0; i < innerlist->GetEntriesFast(); i++) { PndOnlineSttTriplet* innertriplet = (PndOnlineSttTriplet *)(innerlist->At(i)); for (int j = 0; j < outerlist->GetEntriesFast(); j++) { //cout << "Inner: " << i << " Outer: " << j << endl; PndOnlineSttTriplet* outertriplet = (PndOnlineSttTriplet *)(outerlist->At(j)); if ((innertriplet->GetMinT0() > outertriplet->GetMaxT0()) || (innertriplet->GetMaxT0() < outertriplet->GetMinT0()) ) { //cout << "Triplet Times too far apart. Skipping Track. Inner Triplet: " << innertriplet->GetMinT0() << " - " << innertriplet->GetMaxT0() << " Outer Triplet: " << outertriplet->GetMinT0() << " - " << outertriplet->GetMaxT0() << endl; continue; } //determine z range of inner and outer skewlet //find poca of inner and outer skewlet //cout << "Inner Pivot: " << innertriplet->GetPivotStraw() << ", Outer Pivot: " << outertriplet->GetPivotStraw() << endl; //cout << tube_array->At(innertriplet->GetPivotStraw()) << endl << tube_array->At(outertriplet->GetPivotStraw()) << endl; if (innertriplet->GetHitCount() + outertriplet->GetHitCount() < 4) { continue; } TVector2 innerzminmax; TVector3 innerpoint = innertriplet->GetCMS3D(tube_array, &innerzminmax); TVector3 innerdirection = ((PndSttTube*)tube_array->At(innertriplet->GetPivotStraw()))->GetWireDirection(); TVector2 outerzminmax; TVector3 outerpoint = outertriplet->GetCMS3D(tube_array, &outerzminmax); TVector3 outerdirection = ((PndSttTube*)tube_array->At(outertriplet->GetPivotStraw()))->GetWireDirection(); TVector3 poca = TVector3(0.0, 0.0, 0.0); if ((innerpoint - outerpoint).Mag() > 10) { continue; } FindCrossing(innerpoint, innerdirection, outerpoint, outerdirection, poca); /* cout << "Skewlet Combiner: " << endl; innerzminmax.Print(); outerzminmax.Print(); innerpoint.Print(); outerpoint.Print(); innerdirection.Print(); outerdirection.Print(); poca.Print(); cout << (innerpoint-outerpoint).Mag() << endl; */ //check if poca is valid Double_t zmin = innerzminmax.X(); Double_t zmax = innerzminmax.Y(); if (outerzminmax.X() > zmin) { zmin = outerzminmax.X(); } if (outerzminmax.Y() < zmax) { zmax = outerzminmax.Y(); } if ((poca.Z() > zmin) && (poca.Z() < zmax) ) { //cout << "Skewlet valid. Adding to list." << endl; //poca.Print(); //generate full triplet //add full triplet to list PndOnlineSttTriplet* fulltriplet = new PndOnlineSttTriplet; *fulltriplet = *innertriplet; fulltriplet->AddTriplet(*outertriplet); fulltriplet->SetPoca(poca); fulltripletlist->AddLast(fulltriplet); //add full triplet to onlinemanager PndOnlineTrack* new_track = new PndOnlineTrack(); new_track->SetVertex(poca); //new_track->SetT0(mytriplet->GetMinT0()); if (innertriplet->GetMinT0() < outertriplet->GetMinT0()) { new_track->SetT0(innertriplet->GetMinT0()); } else { new_track->SetT0(outertriplet->GetMinT0()); } new_track->SetGoodness(0); new_track->SetDrawMode(PndOnlineTrack::kVertexAsPoint); if (verbosity > 3) { cout << "Adding Skewlet Track to OnlineManager, Timestamp: " << new_track->T0() << endl; } //cout << "Adding Skewlet Track to OnlineManager, Timestamp: " << new_track->T0() << endl; online_manager->AddTrack(new_track); } //tripletvalid } //innerloop } //outerloop return fulltripletlist; } // find skewlets TObjArray* PndOnlineSttTripletFinder::FindSkewlets(vector< list > hit_channel_map, Int_t starttube, Int_t stoptube, Int_t starttubenext, Int_t stoptubenext) { cout << "Finding Skewlets." << endl; TObjArray* skewletlist = NULL; TObjArray sf_innerlist; sf_innerlist.SetName("sf_innerlist"); sf_innerlist.SetOwner(kTRUE); TObjArray sf_innerdirection; sf_innerdirection.SetName("sf_innerdirection"); sf_innerdirection.SetOwner(kTRUE); TObjArray sf_outerlist; sf_outerlist.SetName("sf_outerlist"); sf_outerlist.SetOwner(kTRUE); TObjArray sf_outerdirection; sf_outerdirection.SetName("sf_outerdirection"); sf_outerdirection.SetOwner(kTRUE); std::vector sf_innerhits; std::vector sf_outerhits; std::vector sf_innert0; std::vector sf_outert0; TObjArray sf_innerzminmax; sf_innerzminmax.SetName("sf_innerzminmax"); sf_innerzminmax.SetOwner(kTRUE); TObjArray sf_outerzminmax; sf_outerzminmax.SetName("sf_outerzminmax"); sf_outerzminmax.SetOwner(kTRUE); // find layer skewlets for (int t = starttube; t < stoptube; t++) { //cout << " checking tube " << the_tube << "..."; if (hit_channel_map[t].size() == 0) { //cout << " empty." << endl; continue; } PndOnlineTrackHit* onlinehit = hit_channel_map[t].front(); PndSttHit* drawhit = (PndSttHit*)(onlinehit->BaseHit()); if (!drawhit) { cout << "Error: No Hit Pointer in Stack." << endl; continue; } //cout << "Tube ID: " << drawhit->GetTubeID() << endl; PndSttTube* tube = (PndSttTube*)(tube_array->At(drawhit->GetTubeID())); if (tube == 0) { cout << "Error: Hit pointed to invalid tube." << endl; continue; } // cout << tube << endl; TVector3 tubeposition = tube->GetPosition(); if ((tubeposition.Z() > 36) || (tubeposition.Z() < 34)) { // cout << "Preliminary Short Tube Processing Engaged." << endl; //continue; } if ((drawhit->GetTubeID()>starttube) && (drawhit->GetTubeID()GetTubeID(), hit_channel_map, tube_array, cmspoint, zminmax); sf_innerlist.AddLast(cmspoint.Clone()); sf_innerdirection.AddLast(tube->GetWireDirection().Clone()); sf_innerhits.push_back(neighbors); sf_innert0.push_back(drawhit->GetTimeStamp()); sf_innerzminmax.AddLast(zminmax.Clone()); cout << "Inner Skewlet found. Adding to list." << endl; //skewletlist->AddLast(cmspoint.Clone()); } } for (int t = starttubenext; t < stoptubenext; t++) { //cout << " checking tube " << the_tube << "..."; if (hit_channel_map[t].size() == 0) { //cout << " empty." << endl; continue; } PndOnlineTrackHit* onlinehit = hit_channel_map[t].front(); PndSttHit* drawhit = (PndSttHit*)(onlinehit->BaseHit()); if (!drawhit) { cout << "Error: No Hit Pointer in Stack." << endl; continue; } //cout << "Tube ID: " << drawhit->GetTubeID() << endl; PndSttTube* tube = (PndSttTube*)(tube_array->At(drawhit->GetTubeID())); if (tube == 0) { cout << "Error: Hit pointed to invalid tube." << endl; continue; } // cout << tube << endl; TVector3 tubeposition = tube->GetPosition(); if ((tubeposition.Z() > 36) || (tubeposition.Z() < 34)) { // cout << "Preliminary Short Tube Processing Engaged." << endl; //continue; } if ((drawhit->GetTubeID()>starttubenext) && (drawhit->GetTubeID()GetTubeID(), hit_channel_map, tube_array, cmspointnext, zminmax); sf_outerlist.AddLast(cmspointnext.Clone()); sf_outerdirection.AddLast(tube->GetWireDirection().Clone()); sf_outerhits.push_back(neighbors); sf_outert0.push_back(drawhit->GetTimeStamp()); sf_outerzminmax.AddLast(zminmax.Clone()); cout << "Outer Skewlet found. Adding to list." << endl; //skewletlist->AddLast(cmspointnext.Clone()); } } cout << "Found " << sf_innerlist.GetEntriesFast() << " inner skewlets and " << sf_outerlist.GetEntriesFast() << " outer skewlets." << endl; // combine layer skewlets to 3D skewlets for (int i = 0; i < sf_innerlist.GetEntriesFast(); i++) { for (int j = 0; j < sf_outerlist.GetEntriesFast(); j++) { TVector3 poca; TVector3* innerpoint = (TVector3*)(sf_innerlist.At(i)); TVector3* innerdirection = (TVector3*)(sf_innerdirection.At(i)); TVector3* outerpoint = (TVector3*)(sf_outerlist.At(j)); TVector3* outerdirection = (TVector3*)(sf_outerdirection.At(j)); if (sf_innerhits.at(i) + sf_outerhits.at(j) < 4) { continue; } if (((*innerpoint) - (*outerpoint)).Mag() > 10) { continue; } FindCrossing(*innerpoint, *innerdirection, *outerpoint, *outerdirection, poca); Double_t zmin = ((TVector2*)sf_innerzminmax.At(i))->X(); Double_t zmax = ((TVector2*)sf_innerzminmax.At(i))->Y(); Double_t zminouter = ((TVector2*)sf_outerzminmax.At(j))->X(); Double_t zmaxouter = ((TVector2*)sf_outerzminmax.At(j))->Y(); if (zminouter > zmin) { zmin = zminouter; } if (zmaxouter < zmax) { zmax = zmaxouter; } if ((poca.Z() > zmin) && (poca.Z() < zmax) ) { //cout << "Skewlet valid. Adding to list." << endl; //poca.Print(); if (skewletlist == NULL) { skewletlist = new TObjArray(); skewletlist->SetOwner(kTRUE); } skewletlist->AddLast(poca.Clone()); //temporary direct add of skewlets to track list PndOnlineTrack* new_track = new PndOnlineTrack(); new_track->SetVertex(poca); new_track->SetGoodness(0); if (sf_innert0.at(i) < sf_outert0.at(j)) { new_track->SetT0(sf_innert0.at(i)); } else { new_track->SetT0(sf_outert0.at(j)); } cout << "Adding Skewlet Track to OnlineManager, Timestamp: " << new_track->T0() << endl; online_manager->AddTrack(new_track); } } } // clean up sf_innerlist.Clear(); sf_innerdirection.Clear(); sf_outerlist.Clear(); sf_outerdirection.Clear(); sf_innerzminmax.Clear(); sf_outerzminmax.Clear(); return skewletlist; } // search over given range of tubes, see if any have hits, then look for connected groups of hits TObjArray* PndOnlineSttTripletFinder::FindSkewlets2(vector< list > hit_channel_map, Int_t starttube, Int_t stoptube, Int_t starttubenext, Int_t stoptubenext) { if (verbosity > 2) { cout << "PndOnlineSTTTripletFinder::FindSkewlets2() ... " << endl << " number of channels in hit_channel_map... " << hit_channel_map.size() << endl; } cout << "Finding Skewlets." << endl; TObjArray* skewletlist = NULL; TObjArray sf_innerlist; sf_innerlist.SetName("sf_innerlist"); sf_innerlist.SetOwner(kTRUE); TObjArray sf_innerdirection; sf_innerdirection.SetName("sf_innerdirection"); sf_innerdirection.SetOwner(kTRUE); TObjArray sf_outerlist; sf_outerlist.SetName("sf_outerlist"); sf_outerlist.SetOwner(kTRUE); TObjArray sf_outerdirection; sf_outerdirection.SetName("sf_outerdirection"); sf_outerdirection.SetOwner(kTRUE); std::vector sf_innerhits; std::vector sf_outerhits; std::vector sf_innert0; std::vector sf_outert0; TObjArray sf_innerzminmax; sf_innerzminmax.SetName("sf_innerzminmax"); sf_innerzminmax.SetOwner(kTRUE); TObjArray sf_outerzminmax; sf_outerzminmax.SetName("sf_outerzminmax"); sf_outerzminmax.SetOwner(kTRUE); // find layer skewlets for (int t = starttube; t < stoptube; t++) { //cout << " checking tube " << the_tube << "..."; if (hit_channel_map[t].size() == 0) { //cout << " empty." << endl; continue; } PndOnlineTrackHit* onlinehit = hit_channel_map[t].front(); PndSttHit* drawhit = (PndSttHit*)(onlinehit->BaseHit()); if (!drawhit) { cout << "Error: No Hit Pointer in Stack." << endl; continue; } //cout << "Tube ID: " << drawhit->GetTubeID() << endl; PndSttTube* tube = (PndSttTube*)(tube_array->At(drawhit->GetTubeID())); if (tube == 0) { cout << "Error: Hit pointed to invalid tube." << endl; continue; } list::iterator iter=hit_channel_map[t].begin(); // cout << tube << endl; TVector3 tubeposition = tube->GetPosition(); if ((tubeposition.Z() > 36) || (tubeposition.Z() < 34)) { // cout << "Preliminary Short Tube Processing Engaged." << endl; //continue; } if ((drawhit->GetTubeID()>starttube) && (drawhit->GetTubeID()GetTubeID(), hit_channel_map, iter, tube_array, cmspoint, zminmax); sf_innerlist.AddLast(cmspoint.Clone()); sf_innerdirection.AddLast(tube->GetWireDirection().Clone()); sf_innerhits.push_back(neighbors); sf_innert0.push_back(drawhit->GetTimeStamp()); sf_innerzminmax.AddLast(zminmax.Clone()); cout << "Inner Skewlet found. Adding to list." << endl; //skewletlist->AddLast(cmspoint.Clone()); } } for (int t = starttubenext; t < stoptubenext; t++) { //cout << " checking tube " << the_tube << "..."; if (hit_channel_map[t].size() == 0) { //cout << " empty." << endl; continue; } PndOnlineTrackHit* onlinehit = hit_channel_map[t].front(); PndSttHit* drawhit = (PndSttHit*)(onlinehit->BaseHit()); if (!drawhit) { cout << "Error: No Hit Pointer in Stack." << endl; continue; } //cout << "Tube ID: " << drawhit->GetTubeID() << endl; PndSttTube* tube = (PndSttTube*)(tube_array->At(drawhit->GetTubeID())); if (tube == 0) { cout << "Error: Hit pointed to invalid tube." << endl; continue; } // cout << tube << endl; TVector3 tubeposition = tube->GetPosition(); if ((tubeposition.Z() > 36) || (tubeposition.Z() < 34)) { // cout << "Preliminary Short Tube Processing Engaged." << endl; //continue; } if ((drawhit->GetTubeID()>starttubenext) && (drawhit->GetTubeID()GetTubeID(), hit_channel_map, tube_array, cmspointnext, zminmax); sf_outerlist.AddLast(cmspointnext.Clone()); sf_outerdirection.AddLast(tube->GetWireDirection().Clone()); sf_outerhits.push_back(neighbors); sf_outert0.push_back(drawhit->GetTimeStamp()); sf_outerzminmax.AddLast(zminmax.Clone()); cout << "Outer Skewlet found. Adding to list." << endl; //skewletlist->AddLast(cmspointnext.Clone()); } } cout << "Found " << sf_innerlist.GetEntriesFast() << " inner skewlets and " << sf_outerlist.GetEntriesFast() << " outer skewlets." << endl; // combine layer skewlets to 3D skewlets for (int i = 0; i < sf_innerlist.GetEntriesFast(); i++) { for (int j = 0; j < sf_outerlist.GetEntriesFast(); j++) { TVector3 poca; TVector3* innerpoint = (TVector3*)(sf_innerlist.At(i)); TVector3* innerdirection = (TVector3*)(sf_innerdirection.At(i)); TVector3* outerpoint = (TVector3*)(sf_outerlist.At(j)); TVector3* outerdirection = (TVector3*)(sf_outerdirection.At(j)); if (sf_innerhits.at(i) + sf_outerhits.at(j) < 4) { continue; } if (((*innerpoint) - (*outerpoint)).Mag() > 10) { continue; } FindCrossing(*innerpoint, *innerdirection, *outerpoint, *outerdirection, poca); Double_t zmin = ((TVector2*)sf_innerzminmax.At(i))->X(); Double_t zmax = ((TVector2*)sf_innerzminmax.At(i))->Y(); Double_t zminouter = ((TVector2*)sf_outerzminmax.At(j))->X(); Double_t zmaxouter = ((TVector2*)sf_outerzminmax.At(j))->Y(); if (zminouter > zmin) { zmin = zminouter; } if (zmaxouter < zmax) { zmax = zmaxouter; } if ((poca.Z() > zmin) && (poca.Z() < zmax) ) { //cout << "Skewlet valid. Adding to list." << endl; //poca.Print(); if (skewletlist == NULL) { skewletlist = new TObjArray(); skewletlist->SetOwner(kTRUE); } skewletlist->AddLast(poca.Clone()); //temporary direct add of skewlets to track list PndOnlineTrack* new_track = new PndOnlineTrack(); new_track->SetVertex(poca); new_track->SetGoodness(0); if (sf_innert0.at(i) < sf_outert0.at(j)) { new_track->SetT0(sf_innert0.at(i)); } else { new_track->SetT0(sf_outert0.at(j)); } cout << "Adding Skewlet Track to OnlineManager, Timestamp: " << new_track->T0() << endl; online_manager->AddTrack(new_track); } } } // clean up sf_innerlist.Clear(); sf_innerdirection.Clear(); sf_outerlist.Clear(); sf_outerdirection.Clear(); sf_innerzminmax.Clear(); sf_outerzminmax.Clear(); return skewletlist; } //void PndOnlineSttTripletFinder::ProcessSkewPivotTubeHit(Int_t pivotID, vector > hit_channel_map, list::iterator pivotIter, TObjArray* tripletlist) //{ // UInt_t neighbors = 0; // UInt_t rightID = 0; // Double_t cmsx = 0; // Double_t cmsy = 0; // Double_t cmsz = 0; // cmspoint.SetXYZ(cmsx, cmsy, cmsz); // Double_t zmin = 35-75; // Double_t zmax = 35+75; // zminmax.Set(zmin, zmax); // // PndOnlineSttTriplet* myskewlet = new PndOnlineSttTriplet; // PndSttHit* pivotHit = (PndSttHit*)((*pivotIter)->BaseHit()); // PndSttTube* pivotTube = (PndSttTube*)(tube_array->At(pivotHit->GetTubeID())); // myskewlet->SetPivotHit(pivotHit, pivotTube); // Double_t pivotTime = pivotHit->GetTimeStamp(); // // // imbue MC event header information for comparison // Double_t mcEventTime = 0; // Int_t mcEventID = 0; // Int_t mcTrackID = 0; // Int_t mcSttPointID = 0; // ExtractMCLinkInfo(pivotHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); // // for(int t = 0; t < hit_channel_map.size(); t++) { // if(hit_channel_map[t].size() == 0) { // continue; // } // //enable only if start parameters are set to leftID data // //if(leftID == the_tube) { // // continue; // //} // //PndOnlineTrackHit* onlinehit = hit_channel_map[t].front(); // //PndSttHit* stthit = (PndSttHit*)(onlinehit->BaseHit()); // rightID = t; // PndSttTube* leftTube = (PndSttTube*)stttubes->At(leftID); // PndSttTube* rightTube = (PndSttTube*)stttubes->At(rightID); //// cout << "Checking distance between tubes: " << leftID << ", " << rightID << endl; // Double_t wireanglex = TMath::Abs(leftTube->GetWireDirection().X() - rightTube->GetWireDirection().X()); // Double_t wireangley = TMath::Abs(leftTube->GetWireDirection().Y() - rightTube->GetWireDirection().Y()); //// cout << "Angle between: " << wireangle << endl; //// cout << "Rad In, Rad Out for leftID: " << leftTube->GetRadIn() << " " << leftTube->GetRadOut() << endl; //// cout << "Rad In, Rad Out for rightID: " << rightTube->GetRadIn() << " " << rightTube->GetRadOut() << endl; // if ((wireanglex < 0.005) && (wireangley < 0.005)) { // if (GetTubeDist3D(leftID, rightID, stttubes) < 1.1) { // //cout << neighbors << " Same Skew Neighbor found between tubes: " << leftID << ", " << rightID << ", " << GetTubeDist3D(leftID, rightID, stttubes) << endl; // neighbors++; //// rightTube->GetPosition().Print(); //// rightTube->GetWireDirection().Print(); // Double_t tubex = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().X(); // Double_t tubey = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().Y(); // Double_t tubez = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().Z(); // if ((tubez > 36) || (tubez < 34)) { // Double_t dz = ((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Z(); // Double_t k = (35.0 - tubez)/dz; // Double_t halflen = ((PndSttTube*)(stttubes->At(rightID)))->GetHalfLength(); // if (tubez - dz*halflen > zmin) { // zmin = tubez - dz*halflen; // } // if (tubez + dz*halflen < zmax) { // zmax = tubez + dz*halflen; // } // tubex += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().X(); // tubey += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Y(); // tubez += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Z(); // } // cmsx += tubex; // cmsy += tubey; // cmsz += tubez; // } // } // } // zminmax.Set(zmin, zmax); //// cout << "Returning CMS Point." << endl; // cmspoint.SetXYZ(cmsx/neighbors, cmsy/neighbors, cmsz/neighbors); //// cmspoint->Print(); // return neighbors; // // //TEllipse mylipse; // //mylipse.SetLineColor(kBlue); // //mylipse.SetFillColor(kBlue); // //mylipse.SetFillStyle(1001); // //Double_t myx = ((PndSttTube*)(stttubes->At(leftID)))->GetPosition().X(); // //Double_t myy = ((PndSttTube*)(stttubes->At(leftID)))->GetPosition().Y(); // //mylipse.DrawEllipse(myx, myy, 0.5, 0.5, 0, 360, 0); // //mylipse.SetX1(cmsx/neighbors); // //mylipse.SetY1(cmsy/neighbors); // //mylipse.SetR1(0.4); // //mylipse.SetR2(0.4); // //mylipse.DrawClone(); //} ////////////////////////////////////////////////////////////////////// //////////////////// SKEW HELPERS ////////////////////////////////////////////////////////////////////// UInt_t PndOnlineSttTripletFinder::FindSameSkewNeighborCMS3D(UInt_t leftID, vector< list > hit_channel_map, TObjArray* stttubes, TVector3& cmspoint, TVector2& zminmax) { UInt_t neighbors = 0; UInt_t rightID = 0; Double_t cmsx = 0; Double_t cmsy = 0; Double_t cmsz = 0; cmspoint.SetXYZ(cmsx, cmsy, cmsz); Double_t zmin = 35-75; Double_t zmax = 35+75; zminmax.Set(zmin, zmax); PndOnlineSttTriplet* myskewlet = new PndOnlineSttTriplet; //PndSttHit* pivotHit = (PndSttHit*)((*pivotIter)->BaseHit()); PndOnlineTrackHit* onlinehit = hit_channel_map[leftID].front(); PndSttHit* pivotHit = (PndSttHit*)(onlinehit->BaseHit()); PndSttTube* pivotTube = (PndSttTube*)(tube_array->At(pivotHit->GetTubeID())); myskewlet->SetPivotHit(pivotHit, pivotTube); Double_t pivotTime = pivotHit->GetTimeStamp(); // imbue MC event header information for comparison Double_t mcEventTime = 0; Int_t mcEventID = 0; Int_t mcTrackID = 0; Int_t mcSttPointID = 0; ExtractMCLinkInfo(pivotHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); for(int t = 0; t < hit_channel_map.size(); t++) { if(hit_channel_map[t].size() == 0) { continue; } //enable only if start parameters are set to leftID data //if(leftID == the_tube) { // continue; //} //PndOnlineTrackHit* onlinehit = hit_channel_map[t].front(); //PndSttHit* stthit = (PndSttHit*)(onlinehit->BaseHit()); rightID = t; PndSttTube* leftTube = (PndSttTube*)stttubes->At(leftID); PndSttTube* rightTube = (PndSttTube*)stttubes->At(rightID); // cout << "Checking distance between tubes: " << leftID << ", " << rightID << endl; Double_t wireanglex = TMath::Abs(leftTube->GetWireDirection().X() - rightTube->GetWireDirection().X()); Double_t wireangley = TMath::Abs(leftTube->GetWireDirection().Y() - rightTube->GetWireDirection().Y()); // cout << "Angle between: " << wireangle << endl; // cout << "Rad In, Rad Out for leftID: " << leftTube->GetRadIn() << " " << leftTube->GetRadOut() << endl; // cout << "Rad In, Rad Out for rightID: " << rightTube->GetRadIn() << " " << rightTube->GetRadOut() << endl; if ((wireanglex < 0.005) && (wireangley < 0.005)) { if (GetTubeDist3D(leftID, rightID, stttubes) < 1.1) { //cout << neighbors << " Same Skew Neighbor found between tubes: " << leftID << ", " << rightID << ", " << GetTubeDist3D(leftID, rightID, stttubes) << endl; neighbors++; // rightTube->GetPosition().Print(); // rightTube->GetWireDirection().Print(); Double_t tubex = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().X(); Double_t tubey = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().Y(); Double_t tubez = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().Z(); if ((tubez > 36) || (tubez < 34)) { Double_t dz = ((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Z(); Double_t k = (35.0 - tubez)/dz; Double_t halflen = ((PndSttTube*)(stttubes->At(rightID)))->GetHalfLength(); if (tubez - dz*halflen > zmin) { zmin = tubez - dz*halflen; } if (tubez + dz*halflen < zmax) { zmax = tubez + dz*halflen; } tubex += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().X(); tubey += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Y(); tubez += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Z(); } cmsx += tubex; cmsy += tubey; cmsz += tubez; } } } zminmax.Set(zmin, zmax); // cout << "Returning CMS Point." << endl; cmspoint.SetXYZ(cmsx/neighbors, cmsy/neighbors, cmsz/neighbors); // cmspoint->Print(); return neighbors; //TEllipse mylipse; //mylipse.SetLineColor(kBlue); //mylipse.SetFillColor(kBlue); //mylipse.SetFillStyle(1001); //Double_t myx = ((PndSttTube*)(stttubes->At(leftID)))->GetPosition().X(); //Double_t myy = ((PndSttTube*)(stttubes->At(leftID)))->GetPosition().Y(); //mylipse.DrawEllipse(myx, myy, 0.5, 0.5, 0, 360, 0); //mylipse.SetX1(cmsx/neighbors); //mylipse.SetY1(cmsy/neighbors); //mylipse.SetR1(0.4); //mylipse.SetR2(0.4); //mylipse.DrawClone(); } UInt_t PndOnlineSttTripletFinder::FindSameSkewNeighborCMS3D2(UInt_t leftID, vector< list > hit_channel_map, list::iterator pivotIter, TObjArray* stttubes, TVector3& cmspoint, TVector2& zminmax) { UInt_t neighbors = 0; UInt_t rightID = 0; Double_t cmsx = 0; Double_t cmsy = 0; Double_t cmsz = 0; cmspoint.SetXYZ(cmsx, cmsy, cmsz); Double_t zmin = 35-75; Double_t zmax = 35+75; zminmax.Set(zmin, zmax); PndOnlineSttTriplet* myskewlet = new PndOnlineSttTriplet; PndSttHit* pivotHit = (PndSttHit*)((*pivotIter)->BaseHit()); PndSttTube* pivotTube = (PndSttTube*)(tube_array->At(pivotHit->GetTubeID())); myskewlet->SetPivotHit(pivotHit, pivotTube); Double_t pivotTime = pivotHit->GetTimeStamp(); // imbue MC event header information for comparison Double_t mcEventTime = 0; Int_t mcEventID = 0; Int_t mcTrackID = 0; Int_t mcSttPointID = 0; ExtractMCLinkInfo(pivotHit, mcEventTime, mcEventID, mcTrackID, mcSttPointID); for(int t = 0; t < hit_channel_map.size(); t++) { if(hit_channel_map[t].size() == 0) { continue; } //enable only if start parameters are set to leftID data //if(leftID == the_tube) { // continue; //} //PndOnlineTrackHit* onlinehit = hit_channel_map[t].front(); //PndSttHit* stthit = (PndSttHit*)(onlinehit->BaseHit()); rightID = t; PndSttTube* leftTube = (PndSttTube*)stttubes->At(leftID); PndSttTube* rightTube = (PndSttTube*)stttubes->At(rightID); // cout << "Checking distance between tubes: " << leftID << ", " << rightID << endl; Double_t wireanglex = TMath::Abs(leftTube->GetWireDirection().X() - rightTube->GetWireDirection().X()); Double_t wireangley = TMath::Abs(leftTube->GetWireDirection().Y() - rightTube->GetWireDirection().Y()); // cout << "Angle between: " << wireangle << endl; // cout << "Rad In, Rad Out for leftID: " << leftTube->GetRadIn() << " " << leftTube->GetRadOut() << endl; // cout << "Rad In, Rad Out for rightID: " << rightTube->GetRadIn() << " " << rightTube->GetRadOut() << endl; if ((wireanglex < 0.005) && (wireangley < 0.005)) { if (GetTubeDist3D(leftID, rightID, stttubes) < 1.1) { //cout << neighbors << " Same Skew Neighbor found between tubes: " << leftID << ", " << rightID << ", " << GetTubeDist3D(leftID, rightID, stttubes) << endl; neighbors++; // rightTube->GetPosition().Print(); // rightTube->GetWireDirection().Print(); Double_t tubex = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().X(); Double_t tubey = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().Y(); Double_t tubez = ((PndSttTube*)(stttubes->At(rightID)))->GetPosition().Z(); if ((tubez > 36) || (tubez < 34)) { Double_t dz = ((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Z(); Double_t k = (35.0 - tubez)/dz; Double_t halflen = ((PndSttTube*)(stttubes->At(rightID)))->GetHalfLength(); if (tubez - dz*halflen > zmin) { zmin = tubez - dz*halflen; } if (tubez + dz*halflen < zmax) { zmax = tubez + dz*halflen; } tubex += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().X(); tubey += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Y(); tubez += k*((PndSttTube*)(stttubes->At(rightID)))->GetWireDirection().Z(); } cmsx += tubex; cmsy += tubey; cmsz += tubez; } } } zminmax.Set(zmin, zmax); // cout << "Returning CMS Point." << endl; cmspoint.SetXYZ(cmsx/neighbors, cmsy/neighbors, cmsz/neighbors); // cmspoint->Print(); return neighbors; //TEllipse mylipse; //mylipse.SetLineColor(kBlue); //mylipse.SetFillColor(kBlue); //mylipse.SetFillStyle(1001); //Double_t myx = ((PndSttTube*)(stttubes->At(leftID)))->GetPosition().X(); //Double_t myy = ((PndSttTube*)(stttubes->At(leftID)))->GetPosition().Y(); //mylipse.DrawEllipse(myx, myy, 0.5, 0.5, 0, 360, 0); //mylipse.SetX1(cmsx/neighbors); //mylipse.SetY1(cmsy/neighbors); //mylipse.SetR1(0.4); //mylipse.SetR2(0.4); //mylipse.DrawClone(); } Double_t PndOnlineSttTripletFinder::FindCrossing(const TVector3& p1, const TVector3& u1, const TVector3& p2, const TVector3& u2, TVector3& poca) { // cout << "Crosspoint Input Vectors: " << endl; // p1.Print(); // u1.Print(); // p2.Print(); // u2.Print(); TVector3 p21 = p2 - p1; if (p21.Mag() < 0.000001) { poca = p1; return 0; } // p21.Print(); TVector3 m = u2.Cross(u1); if (m.Mag() < 0.000001) { poca = p1 + 0.5*p21; return p21.Mag(); } // m.Print(); TVector3 mnorm = (1.0/m.Mag2())*m; TVector3 r = p21.Cross(mnorm); Double_t t1 = r.Dot(u2); Double_t t2 = r.Dot(u1); TVector3 q1 = p1+t1*u1; TVector3 q2 = p2+t2*u2; TVector3 q21 = q2 - q1; poca = q1 + 0.5*q21; // cout << "Crosspoints: " << endl; // q1.Print(); // q2.Print(); // cout << "Poca: " << endl; // poca.Print(); return (TMath::Abs(p21.Dot(m))/m.Mag()); } // get tube distance in the z = 35 cm xy projection Double_t PndOnlineSttTripletFinder::GetTubeDist3D(Int_t leftID, Int_t rightID, TObjArray* stttubes) { PndSttTube* lefttube = (PndSttTube*)(stttubes->At(leftID)); PndSttTube* righttube = (PndSttTube*)(stttubes->At(rightID)); Double_t leftx = lefttube->GetPosition().X(); Double_t lefty = lefttube->GetPosition().Y(); Double_t leftz = lefttube->GetPosition().Z(); Double_t rightx = righttube->GetPosition().X(); Double_t righty = righttube->GetPosition().Y(); Double_t rightz = righttube->GetPosition().Z(); if ((leftz > 36) || (leftz < 34)) { Double_t dz = lefttube->GetWireDirection().Z(); Double_t k = (35.0 - leftz)/dz; leftx += k*lefttube->GetWireDirection().X(); lefty += k*lefttube->GetWireDirection().Y(); leftz += k*lefttube->GetWireDirection().Z(); } if ((rightz > 36) || (rightz < 34)) { Double_t dz = righttube->GetWireDirection().Z(); Double_t k = (35.0 - rightz)/dz; rightx += k*righttube->GetWireDirection().X(); righty += k*righttube->GetWireDirection().Y(); rightz += k*righttube->GetWireDirection().Z(); } Double_t dx = rightx - leftx; Double_t dy = righty - lefty; return (TMath::Sqrt(dx*dx+dy*dy)); } ////////////////////////////////////////////////////////////////////// //////////////////// HELPERS ////////////////////////////////////////////////////////////////////// void PndOnlineSttTripletFinder::CalculateCircle(const PndOnlineSttTriplet* const innertriplet, const PndOnlineSttTriplet* const outertriplet, Double_t& originx, Double_t& originy, Double_t& radius, Int_t& clockwise) { Double_t x1= innertriplet->GetCMS(tube_array).X(); Double_t y1= innertriplet->GetCMS(tube_array).Y(); Double_t x2= outertriplet->GetCMS(tube_array).X(); Double_t y2= outertriplet->GetCMS(tube_array).Y(); Double_t a= x1; Double_t b= y1; Double_t c= x2; Double_t d= y2; Double_t px = 0.5*(b*c*c-a*a*d-b*b*d+b*d*d)/(b*c-a*d); Double_t py = 0.5*(a*a*c+b*b*c-a*c*c-a*d*d)/(b*c-a*d); Double_t r = TMath::Sqrt(px*px+py*py); //circleorigin->Set(px, py); //*circleradius = r; originx = px; originy = py; radius = r; TVector2 invec(x1, y1); TVector2 outvec(x2, y2); TVector2 rotvec = outvec.Rotate(-1*invec.Phi()); Bool_t clockwisecheck = rotvec.Phi() > TMath::Pi(); Bool_t isclockwise = (y1 > x1*y2/x2); if (x2 < 0) { isclockwise = !(isclockwise); } if (isclockwise) { //*clockwise = 1; clockwise = 1; } else { //*clockwise = -1; clockwise = -1; } if (clockwisecheck != isclockwise) { if (verbosity > 0) { cout << "ERROR: Clockwise Determination Mismatch!" << endl; } } } Bool_t PndOnlineSttTripletFinder::HitInRange(Double_t originx, Double_t originy, Double_t radius, Int_t clockwise, PndSttHit* stthit, Double_t maxdistance) { PndSttTube* tube = (PndSttTube*)(tube_array->At(stthit->GetTubeID())); TVector2 hitpos(tube->GetPosition().X(), tube->GetPosition().Y()); TVector2 circleorigin(originx, originy); TVector2 distvec = hitpos - circleorigin; TVector2 originvec = (-1)* circleorigin; if (TMath::Abs(distvec.Mod() - radius) < maxdistance) { //cout << "Clockwise: " << circleclockwise << endl; if (originvec.DeltaPhi(distvec)*clockwise > 0) { //hitcount++; return kTRUE; } } return kFALSE; } Bool_t PndOnlineSttTripletFinder::HitInRange(Double_t originx, Double_t originy, Double_t radius, Int_t clockwise, Int_t tubeID, Double_t maxdistance) { PndSttTube* tube = (PndSttTube*)(tube_array->At(tubeID)); TVector2 hitpos(tube->GetPosition().X(), tube->GetPosition().Y()); TVector2 circleorigin(originx, originy); TVector2 distvec = hitpos - circleorigin; TVector2 originvec = (-1)* circleorigin; /* cout << "Checking Hit Distance" << std::endl; cout << "Hit Position: "; hitpos.Print(); cout << "Circle Origin: "; circleorigin.Print(); cout << "Distance Vector: "; distvec.Print(); */ if (TMath::Abs(distvec.Mod() - radius) < maxdistance) { //cout << "Clockwise: " << circleclockwise << endl; if (originvec.DeltaPhi(distvec)*clockwise > 0) { //hitcount++; return kTRUE; } } return kFALSE; } Bool_t PndOnlineSttTripletFinder::HitInRange(Double_t originx, Double_t originy, Double_t radius, Int_t clockwise, const TVector2& hitpos, Double_t maxdistance) { //PndSttTube* tube = (PndSttTube*)(tube_array->At(tubeID)); //TVector2 hitpos(tube->GetPosition().X(), tube->GetPosition().Y()); TVector2 circleorigin(originx, originy); TVector2 distvec = hitpos - circleorigin; TVector2 originvec = (-1)* circleorigin; /* cout << "Checking Hit Distance" << std::endl; cout << "Hit Position: "; hitpos.Print(); cout << "Circle Origin: "; circleorigin.Print(); cout << "Distance Vector: "; distvec.Print(); */ if (TMath::Abs(distvec.Mod() - radius) < maxdistance) { //cout << "Clockwise: " << circleclockwise << endl; if (originvec.DeltaPhi(distvec)*clockwise > 0) { //hitcount++; return kTRUE; } } return kFALSE; } Double_t PndOnlineSttTripletFinder::CalculateCircleProximity(Double_t originx, Double_t originy, Double_t radius, Int_t clockwise, PndSttHit* stthit) { PndSttTube* tube = (PndSttTube*)(tube_array->At(stthit->GetTubeID())); TVector2 hitpos(tube->GetPosition().X(), tube->GetPosition().Y()); TVector2 circleorigin(originx, originy); TVector2 distvec = hitpos - circleorigin; TVector2 originvec = (-1)* circleorigin; if (TMath::Abs(distvec.Mod() - radius) < 0.6) { //cout << "Clockwise: " << circleclockwise << endl; if (originvec.DeltaPhi(distvec)*clockwise > 0) { //hitcount++; } } } // Calculates distance in cm between two Straw Tubes /* Int_t leftID: Index of the first Tube Int_t rightID: Index of the second Tube TObjArray* stttubes: Array storing the setup of the PndSttTube objects comprising the STT returns: Distance in cm between the two given tubes in xy Projection */ Double_t PndOnlineSttTripletFinder::GetTubeDist(Int_t leftID, Int_t rightID, TObjArray* stttubes) { PndSttTube* lefttube = (PndSttTube*)(stttubes->At(leftID)); PndSttTube* righttube = (PndSttTube*)(stttubes->At(rightID)); Double_t dx = righttube->GetPosition().X()-lefttube->GetPosition().X(); Double_t dy = righttube->GetPosition().Y()-lefttube->GetPosition().Y(); return (TMath::Sqrt(dx*dx+dy*dy)); } void PndOnlineSttTripletFinder::ExtractMCLinkInfo(FairMultiLinkedData* myHit, Double_t& MCEventTime, Int_t& MCEventID, Int_t& MCTrackID, Int_t& MCSttPointID) { FairMultiLinkedData myData = myHit->GetLinksWithType(FairRootManager::Instance()->GetBranchId("EventHeader.")); if (myData.GetNLinks() > 0) { FairEventHeader* myHeader = (FairEventHeader*)(FairRootManager::Instance()->GetCloneOfLinkData(myData.GetLink(0))); MCEventTime = myHeader->GetEventTime(); MCEventID = myHeader->GetMCEntryNumber(); } else { if (verbosity > 0) { cout << "Ooops, link to Event Header not found." << endl; } MCEventTime = -2; MCEventID = -2; } PndSttPoint* mySttPoint = 0; FairMultiLinkedData mySttMCLink = myHit->GetLinksWithType(FairRootManager::Instance()->GetBranchId("STTPoint")); if (mySttMCLink.GetNLinks() > 0) { mySttPoint = (PndSttPoint*)(FairRootManager::Instance()->GetCloneOfLinkData(mySttMCLink.GetLink(0))); MCSttPointID = mySttMCLink.GetLink(0).GetIndex(); if (mySttPoint > 0) { MCTrackID = mySttPoint->GetTrackID(); } else { if (verbosity > 0) { cout << "Ooops, SttPoint Pointer is 0." << endl; } MCTrackID = -2; } } else { if (verbosity > 0) { cout << "Ooops, link to STTPoint not found." << endl; } MCTrackID = -2; MCSttPointID = -2; } } // ----- Public method Finish ------------------------------------------ void PndOnlineSttTripletFinder::Finish() { } ClassImp(PndOnlineSttTripletFinder);