//////////////////////////////////////////////////////////// // // PndTrkLegendreTask3 // // Class for secondary track pattern recognition // // authors: Lia Lavezzi - INFN Pavia (2012) // //////////////////////////////////////////////////////////// #include "PndTrkLegendreTask3.h" // stt #include "PndSttHit.h" #include "PndSttTube.h" #include "PndSttMapCreator.h" // sds #include "PndSdsHit.h" // track(cand) #include "PndTrackCand.h" #include "PndTrackCandHit.h" #include "PndTrack.h" #include "PndTrkConformalHit.h" #include "PndTrkConformalHitList.h" #include "PndTrkConformalTransform.h" #include "PndTrkTools.h" #include "PndTrkSkewHit.h" #include "PndTrkFitter.h" // fairroot #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" // ROOT #include "TClonesArray.h" #include "TVector3.h" #include "TArc.h" #include "TLine.h" #include "TMarker.h" #include "TSpectrum2.h" #include "TSpectrum.h" #include "TStopwatch.h" // tracking #include "PndTrkClusterList.h" #include using namespace std; // ----- Default constructor ------------------------------------------- PndTrkLegendreTask3::PndTrkLegendreTask3() : FairTask("secondary track finder", 0), fDisplayOn(kFALSE), fPersistence(kTRUE), fUseMVDPix(kTRUE), fUseMVDStr(kTRUE), fUseSTT(kTRUE), fSecondary(kFALSE), fMvdPix_RealDistLimit(1000), fMvdStr_RealDistLimit(1000), fStt_RealDistLimit(1000), fMvdPix_ConfDistLimit(1000), fMvdStr_ConfDistLimit(1000), fStt_ConfDistLimit(1000) { sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); PndGeoHandling::Instance(); } PndTrkLegendreTask3::PndTrkLegendreTask3(int verbose) : FairTask("secondary track finder", verbose), fDisplayOn(kFALSE), fPersistence(kTRUE), fUseMVDPix(kTRUE), fUseMVDStr(kTRUE), fUseSTT(kTRUE), fSecondary(kFALSE), fMvdPix_RealDistLimit(1000), fMvdStr_RealDistLimit(1000), fStt_RealDistLimit(1000), fMvdPix_ConfDistLimit(1000), fMvdStr_ConfDistLimit(1000), fStt_ConfDistLimit(1000) { sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); PndGeoHandling::Instance(); } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndTrkLegendreTask3::~PndTrkLegendreTask3() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- InitStatus PndTrkLegendreTask3::Init() { fEventCounter = 0; fMvdPix_RealDistLimit = 0.5; // CHECK limits fMvdStr_RealDistLimit = 0.5; // CHECK limits fStt_RealDistLimit = 1.5 * 0.5; // CHECK limits fMvdPix_ConfDistLimit = 0.003; // CHECK limits fMvdStr_ConfDistLimit = 0.003; // CHECK limits fStt_ConfDistLimit = 0.001; // CHECK limits if(fSecondary) { fMvdPix_ConfDistLimit = 0.007; // CHECK limits fMvdStr_ConfDistLimit = 0.007; // CHECK limits } // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndTrkLegendreTask3::Init: " << "RootManager not instantiated, return!" << endl; return kFATAL; } // -- HITS ------------------------------------------------- // // STT fSttHitArray = (TClonesArray*) ioman->GetObject(fSttBranch); if ( ! fSttHitArray ) { cout << "-W- PndTrkLegendreTask3::Init: " << "No STTHit array, return!" << endl; return kERROR; } // // MVD PIXEL fMvdPixelHitArray = (TClonesArray*) ioman->GetObject(fMvdPixelBranch); if ( !fMvdPixelHitArray){ std::cout << "-W- PndTrkLegendreTask3::Init: " << "No MVD Pixel hitArray, return!" << std::endl; return kERROR; } // // MVD STRIP fMvdStripHitArray = (TClonesArray*) ioman->GetObject(fMvdStripBranch); if ( !fMvdStripHitArray){ std::cout << "-W- PndTrkLegendreTask3::Init: " << "No MVD Strip hitArray, return!" << std::endl; return kERROR; } fTrackArray = new TClonesArray("PndTrack"); fTrackCandArray = new TClonesArray("PndTrackCand"); ioman->Register("Track", "pr", fTrackArray, fPersistence); // CHECK ioman->Register("TrackCand", "pr", fTrackCandArray, fPersistence); // CHECK // ---------------------------------------- maps of STT tubes fMapper = new PndSttMapCreator(fSttParameters); fTubeArray = fMapper->FillTubeArray(); // ---------------------------------------------------- end map if(fDisplayOn) { display = new TCanvas("display", "display", 0, 0, 800, 800); // CHECK display->Divide(2, 2); } legendre = new PndTrkLegendreTransform3(); if(fSecondary == kTRUE) legendre->SearchSecondaryTracks(); conform = new PndTrkConformalTransform(); tools = new PndTrkTools(); fFitter = new PndTrkFitter(fVerbose); fClean = new PndTrkClean(fTubeArray); return kSUCCESS; } // ------------------------------------------------------------------------- void PndTrkLegendreTask3::SetParContainers() { FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb(); fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar"); } // ------------------------------------------------------------------------- void PndTrkLegendreTask3::Initialize() { stthitlist = new PndTrkSttHitList(fTubeArray); mvdpixhitlist = new PndTrkSdsHitList(MVDPIXEL); mvdstrhitlist = new PndTrkSdsHitList(MVDSTRIP); if(fUseSTT) { stthitlist->AddTCA(FairRootManager::Instance()->GetBranchId(fSttBranch), fSttHitArray); stthitlist->Instanciate(); } if(fUseMVDPix) { mvdpixhitlist->AddTCA(FairRootManager::Instance()->GetBranchId(fMvdPixelBranch), fMvdPixelHitArray); mvdpixhitlist->InstanciatePixel(); } if(fUseMVDStr) { mvdstrhitlist->AddTCA(FairRootManager::Instance()->GetBranchId(fMvdStripBranch), fMvdStripHitArray); mvdstrhitlist->InstanciateStrip(); } conformalhitlist = new PndTrkConformalHitList(); fFoundPeaks.clear(); } void PndTrkLegendreTask3::Exec(Option_t* opt) { fTrackArray->Delete(); fTrackCandArray->Delete(); if(fVerbose > 0) cout << "*********************** " << fEventCounter << " ***********************" << endl; // CHECK delete this --- // if(fEventCounter == 126 || fEventCounter == 526) { // fEventCounter++; // return; // } // --- fEventCounter++; // CHECK : it may happen that the hit of a track generate more than one track; e.g. // some of them are associated together to form a track, some are left out. If later the left outs // result in a new peak, this new peak will have also the old, already used hits, associated to it. // This happens since the used hits do not enter the new legendre transform, but can be associated // to the found cluster with the distance criterion. In this way the final cluster, which then is // fitted again with the mvd in the legendre plane, will give an existing peak. // This results in an infinite loop ---> FIX IT! (how: timeout? better hit-to-cluster association? // forbidden peak positions in legendre transform?) // CHECK THIS! if(fSttHitArray->GetEntriesFast() < 3) { Reset(); return; } Initialize(); if(fVerbose > 1) { cout << "number of stt hits " << fSttHitArray->GetEntriesFast() << endl; cout << "number of mvdpix hits " << fMvdPixelHitArray->GetEntriesFast() << endl; cout << "number of mvdstr hits " << fMvdStripHitArray->GetEntriesFast() << endl; } if(fDisplayOn) { Refresh(); char goOnChar; // cout << "Start?"; // cin >> goOnChar; display->Update(); display->Modified(); cout << " STARTING" << endl; } // translation and rotation - CHECK no rotation now Double_t delta = 0, trasl[2] = {0., 0.}; if(!fSecondary) { conform->SetOrigin(trasl[0], trasl[1], delta); Int_t nchits = FillConformalHitList(); // if(nchits == 0) return; // CHECK // cout << nchits << " " << trasl[0] << " " << trasl[1] << " " << delta << endl; } int maxpeak = 1000; int ipeak = 0; PndTrkClusterList clusterlist; cout << "STARTING CLUSTERLIST " << clusterlist.GetNofClusters() << endl; std::vector< std::pair > foundpeaks; // LOOP OVER PEAKS UNTILL ITS HEIGHT IS <= 3 fTimer = new TStopwatch(); fTimer->Start(); fTime = 0; int icounter = 0; // while(maxpeak > 3) { // LOOP WHEN DOING SECONDARY SEACH CHECK! // TIME if(fDisplayOn == kFALSE) { fTimer->Stop(); fTime += fTimer->RealTime(); if(fTime > 1.5) { if(fVerbose > 0) cerr << fTime << endl; Reset(); return; } fTimer->Start(); } // if(fSecondary) { // translation and rotation - CHECK fRefHit = FindReferenceHit(); if(fVerbose > 1) cout << "refhit " << fRefHit << endl; if(fRefHit == NULL) { Reset(); return; } ComputeTraAndRot(fRefHit, delta, trasl); conform->SetOrigin(trasl[0], trasl[1], delta); Int_t nchits = FillConformalHitList(); // if(nchits == 0) return; // CHECK if(fVerbose > 1) cout << nchits << " " << trasl[0] << " " << trasl[1] << " " << delta << endl; } PndTrkCluster cluster; if(fDisplayOn) { Refresh(); char goOnChar; cout << endl; cout << "Next peak:? "; cin >> goOnChar; display->Update(); display->Modified(); } if(fVerbose > 1) cout << "@@@@ PEAK No. " << ipeak << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << endl; ipeak++; // APPLY LEGENDRE TO STT ALONE // ------------------------------------------------------- double theta_max, r_max; maxpeak = ApplyLegendre(icounter, theta_max, r_max); cout << "GET NOF MAXIM " << legendre->GetNofMaxima() << endl; // if(maxpeak <= 3) break; for(int imax = 0; imax < legendre->GetNofMaxima() ; imax++) { double thm, rm; PndTrkLegendreCluster thislist = legendre->ExtractLegendreMaximum(imax); thm = thislist.GetTheta(); rm = thislist.GetR(); cout << "MAX " << imax << ": "; for(int ihit = 0; ihit < thislist.GetNofHits(); ihit++) { PndTrkHit *hit = thislist.GetHit(ihit); cout << " (" << hit->GetDetectorID() << ", " << hit->GetHitID() << ")"; } cout << endl; if(fDisplayOn) { char goOnChar; cout << endl; cout << "Next peak? "; cin >> goOnChar; Refresh(); display->cd(3); legendre->DrawMaximum(imax); display->cd(1); legendre->DrawMaximumHits(imax); display->Update(); display->Modified(); } } clusterlist = legendre->GetClusterList(); PndTrkClusterList tmpclusterlist = fClean->MergeClusters(&clusterlist); clusterlist = tmpclusterlist; cout << "************* MERGING DONE ***************" << endl; // treat every cluster PndTrkClusterList clusterlist2; for(int iclus = 0; iclus < clusterlist.GetNofClusters(); iclus++) { double fitm, fitq; PndTrkLegendreCluster *cluster1 = (PndTrkLegendreCluster *) clusterlist.GetCluster(iclus); legendre->ExtractLineParameters(cluster1, fitm, fitq); cout << "cluster " << iclus << " " << cluster1->GetNofHits() << " " << cluster1->GetTheta() << " " << cluster1->GetR() << endl; if(fDisplayOn) { char goOnChar; cout << endl; cout << "Next line" << endl; cin >> goOnChar; Refresh(); display->cd(2); TLine *line = new TLine(-0.07, fitq + fitm * (-0.07), 0.07, fitq + fitm * (0.07)); line->Draw("SAME"); display->Update(); display->Modified(); } // cout << "CLUS1 " << endl; // cluster1->Print(); // ADD HITS TO CLUSTER: PIXEL, STRIP AND STT PARALLEL // ------------------------------------------------------- // method of hit-to-track association: // 0 conformal: mvd and stt in conformal plane // 1 real: mvd and stt in real plane // 2 mixed: mvd in real/stt in conformal plane int method = 0; // CHECK change to AddHitToClusterByDistance PndTrkCluster cluster2 = CreateClusterByDistance(method, fitm, fitq); // CHECK // cout << "CLUS2 " << endl; // cluster2.Print(); PndTrkCluster cluster3 = *cluster1; int nhits = cluster3.MergeTo(&cluster2); // CHECK (delete thsi when AddHitToClusterByDistance) cout << "MERGED " << nhits << endl; cluster3.Print(); if(fDisplayOn) { display->cd(1); cluster3.LightUp(); display->Update(); display->Modified(); } // // APPLY LEGENDRE TO CLUSTER + MVD // // ------------------------------------------------------- // maxpeak = ApplyLegendre(iclus, &cluster, theta_max, r_max); // cout << "(with MVD) GET NOF MAXIM " << legendre->GetNofMaxima() << endl; // for(int imax = 0; imax < legendre->GetNofMaxima() ; imax++) { // double thm, rm; // PndTrkCluster thislist = legendre->ExtractLegendreMaximum(imax, thm, rm); // cout << "MAX " << imax << ": "; // for(int ihit = 0; ihit < thislist.GetNofHits(); ihit++) { // PndTrkHit *hit = thislist.GetHit(ihit); // cout << " (" << hit->GetDetectorID() << ", " << hit->GetHitID() << ")"; // } // cout << endl; // } fClean->Cleanup2(&cluster3); clusterlist2.AddCluster(&cluster3); } // } clusterlist = clusterlist2; if(fDisplayOn) { char goOnChar; display->Update(); display->Modified(); cout << "Finish? "; cin >> goOnChar; cout << "FINISH " << clusterlist.GetNofClusters() << endl; for(int iclus = 0; iclus < clusterlist.GetNofClusters(); iclus++) { Refresh(); PndTrkCluster *clusterb = clusterlist.GetCluster(iclus); clusterb->LightUp(); display->Update(); display->Modified(); cin >> goOnChar; } } Reset(); } void PndTrkLegendreTask3::Reset() { if(fDisplayOn) { char goOnChar; display->Update(); display->Modified(); cout << "Finish? "; cin >> goOnChar; cout << "FINISH" << endl; } if(stthitlist) delete stthitlist; if(mvdpixhitlist) delete mvdpixhitlist; if(mvdstrhitlist) delete mvdstrhitlist; if(fTimer) { fTimer->Stop(); fTime += fTimer->RealTime(); if(fVerbose > 0) cerr << fEventCounter << " Real time " << fTime << " s" << endl; } legendre->ResetLegendre(); } // ============================================================================================ Int_t PndTrkLegendreTask3::FillConformalHitList() { conformalhitlist->SetConformalTransform(conform); // FILL ONCE FOR ALL THE CONFORMAL HIT LIST for(int jhit = 0; jhit < mvdpixhitlist->GetNofHits(); jhit++) { PndTrkHit *hit = mvdpixhitlist->GetHit(jhit); PndTrkConformalHit * chit = conform->GetConformalHit(hit); conformalhitlist->AddHit(chit); } for(int jhit = 0; jhit < mvdstrhitlist->GetNofHits(); jhit++) { PndTrkHit *hit = mvdstrhitlist->GetHit(jhit); PndTrkConformalHit * chit = conform->GetConformalHit(hit); conformalhitlist->AddHit(chit); } for(int jhit = 0; jhit < stthitlist->GetNofHits(); jhit++) { PndTrkHit *hit = stthitlist->GetHit(jhit); if(hit->IsSttSkew()) continue; PndTrkConformalHit * chit = conform->GetConformalSttHit(hit); conformalhitlist->AddHit(chit); } return conformalhitlist->GetNofHits(); } void PndTrkLegendreTask3::FillLegendreHisto(Int_t mode) { // mode 0 STT alone // 1 STT + MVD for(int ihit = 0; ihit < conformalhitlist->GetNofHits(); ihit++) { PndTrkConformalHit *chit = conformalhitlist->GetHit(ihit); if(mode == 0 && chit->GetDetectorID() != FairRootManager::Instance()->GetBranchId(fSttBranch)) continue; PndTrkHit *hit = chit->GetHit(); if(hit->IsUsed()) continue; // CHECK temporary .... int detid = -1; if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) detid = 0; else if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) detid = 1; else if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) detid = 2; // .... legendre->FillLegendreHisto(hit, chit->GetU(), chit->GetV(), chit->GetIsochrone()); if(fDisplayOn) { DrawConfHit(chit->GetU(), chit->GetV(), chit->GetIsochrone()); } } } void PndTrkLegendreTask3::FillLegendreHisto(PndTrkCluster *cluster) { // --------------------------------------------------------------- for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { PndTrkHit *hit = cluster->GetHit(ihit); double conformal[3]; PndTrkConformalHit * chit = conform->GetConformalSttHit(hit); // CHECK temporary .... int detid = -1; if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) detid = 0; else if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) detid = 1; else if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) detid = 2; // .... legendre->FillLegendreHisto(hit, chit->GetU(), chit->GetV(), chit->GetIsochrone()); // conformalhitlist->AddHit(chit); } // add mvd hits to stt cluster for(int ihit = 0; ihit < conformalhitlist->GetNofHits(); ihit++) { PndTrkConformalHit *chit = conformalhitlist->GetHit(ihit); if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) continue; PndTrkHit *hit = chit->GetHit(); if(hit->IsUsed()) continue; // CHECK temporary .... int detid = -1; if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) detid = 0; else if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) detid = 1; else if(chit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) detid = 2; // .... // cout << "FILL " << detid << " " << chit->GetHitID() << " " << chit->GetU() << " " << chit->GetV() << " " << chit->GetIsochrone() << endl; legendre->FillLegendreHisto(hit, chit->GetU(), chit->GetV(), chit->GetIsochrone()); if(fDisplayOn) { DrawConfHit(chit->GetU(), chit->GetV(), chit->GetIsochrone()); } } } PndTrkCluster PndTrkLegendreTask3::CreateClusterByConfDistance(double fitm, double fitq) { PndTrkCluster cluster; // ADD HIT IN CONFORMAL PLANE WITH DISTANCE CRITERION // cost x + sint y - r = 0 -> |cost x0 + sint y0 - r| // cout << "CONFORMALHITLIST " << conformalhitlist->GetNofHits() << endl; for(int ihit = 0; ihit < conformalhitlist->GetNofHits(); ihit++) { PndTrkConformalHit *chit = conformalhitlist->GetHit(ihit); double dist = chit->GetDistanceFromTrack(fitm, fitq); int hitID = chit->GetHitID(); int detID = chit->GetDetectorID(); // CHECK limits double distlimit = 0; if(detID == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) distlimit = 0.003; else if(detID == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) distlimit = 0.003; else if(detID == FairRootManager::Instance()->GetBranchId(fSttBranch)) { if(chit->GetIsochrone() > 1.e-4) distlimit = 6 * chit->GetIsochrone(); else distlimit = 0.001; } // --------------------------- Bool_t accept = kFALSE; if(dist < distlimit) { if(fDisplayOn) { if(chit->GetIsochrone() > 0) { TArc *arc = new TArc(chit->GetU(), chit->GetV(), chit->GetIsochrone()); arc->SetFillStyle(0); arc->SetLineColor(5); display->cd(2); arc->Draw("SAME"); } else { TMarker *mrk; if(detID == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) mrk = new TMarker(chit->GetU(), chit->GetV(), 21); else if(detID == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) mrk = new TMarker(chit->GetU(), chit->GetV(), 25); mrk->SetMarkerColor(5); display->cd(2); mrk->Draw("SAME"); } } PndTrkHit *hit; if(detID == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) { hit = mvdpixhitlist->GetHit(hitID); } else if(detID == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) { hit = mvdstrhitlist->GetHit(hitID); } else if(detID == FairRootManager::Instance()->GetBranchId(fSttBranch)) { hit = stthitlist->GetHit(hitID); } // ********** CHECK *********** // if you put "is used" the efficiency beecomes much lower! // DO NOT: if(!hit->IsUsed()) cluster.AddHit(hit); accept = kTRUE; } // else cout << "DISCARDED " << hitID << " " << detID << " " << dist << " " << distlimit << endl; } return cluster; } PndTrkCluster PndTrkLegendreTask3::CreateClusterByRealDistance(double xc0, double yc0, double R0) { // ADD HIT IN REAL PLANE WITH DISTANCE CRITERION // sqrt((xc - x)**2 + (yc - y)**2) < distlim PndTrkCluster cluster; // ......................................................... if(fDisplayOn) { TArc *arc = new TArc(xc0, yc0, R0); arc->SetFillStyle(0); arc->SetLineColor(5); display->cd(1); arc->Draw("SAME"); } for(int ihit = 0; ihit < stthitlist->GetNofHits(); ihit++) { PndTrkHit *hit = stthitlist->GetHit(ihit); if(hit->IsSttSkew()) continue; double dist = hit->GetXYDistanceFromTrack(xc0, yc0, R0); // cout << ihit << " dist " << dist << endl; // CHECK limits double distlimit = 1.5 * 0.5; if(dist < distlimit){ // cout << "TRACK " << xc0 << " " << yc0 << " " << R0 << endl; // cout << "ADD HIT " << dist << " " << distlimit << " " << hit->GetHitID() << endl; cluster.AddHit(hit); if(fDisplayOn) { TArc *arc = new TArc(hit->GetPosition().X(), hit->GetPosition().Y(), hit->GetIsochrone()); arc->SetFillStyle(0); arc->SetLineColor(5); display->cd(1); arc->Draw("SAME"); } } } for(int ihit = 0; ihit < mvdpixhitlist->GetNofHits(); ihit++) { PndTrkHit *hit = mvdpixhitlist->GetHit(ihit); double dist = hit->GetXYDistanceFromTrack(xc0, yc0, R0); // CHECK limits double distlimit = 0.5; if(dist < distlimit) { cluster.AddHit(hit); if(fDisplayOn) { TMarker *mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 21); mrk->SetMarkerColor(5); display->cd(1); mrk->Draw("SAME"); } } } for(int ihit = 0; ihit < mvdstrhitlist->GetNofHits(); ihit++) { PndTrkHit *hit = mvdstrhitlist->GetHit(ihit); double dist = hit->GetXYDistanceFromTrack(xc0, yc0, R0); // CHECK limits double distlimit = 0.5; if(dist < distlimit){ cluster.AddHit(hit); if(fDisplayOn) { TMarker *mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 25); mrk->SetMarkerColor(5); display->cd(1); mrk->Draw("SAME"); } } } return cluster; } PndTrkCluster PndTrkLegendreTask3::CreateClusterByMixedDistance(double fitm, double fitq) { PndTrkCluster cluster; double delta = conformalhitlist->GetConformalTransform()->GetRotation(); // CHECK // ADD HIT IN MIXED MODE REAL P. FOR MVD & CONFORMAL P. FOR STT, // WITH DISTANCE CRITERION // real: fabs(R - sqrt((xc - x)**2 + (yc - y)**2)) < distlim // cofn: cost x + sint y - r = 0 -> |cost x0 + sint y0 - r| // CHECK if this needs to be kept --> change xc0 to xc etc // center and radius Double_t xc0, yc0, R0; FromConformalToRealTrack(fitm, fitq, xc0, yc0, R0); // ......................................................... if(fDisplayOn) { TArc *arc = new TArc(xc0, yc0, R0); arc->SetFillStyle(0); arc->SetLineColor(5); display->cd(1); arc->Draw("SAME"); } for(int ihit = 0; ihit < mvdpixhitlist->GetNofHits(); ihit++) { PndTrkHit *hit = mvdpixhitlist->GetHit(ihit); double dist = hit->GetXYDistanceFromTrack(xc0, yc0, R0); // CHECK limits double distlimit = 0.5; if(dist < distlimit) { cluster.AddHit(hit); if(fDisplayOn) { TMarker *mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 21); mrk->SetMarkerColor(5); display->cd(1); mrk->Draw("SAME"); } } } for(int ihit = 0; ihit < mvdstrhitlist->GetNofHits(); ihit++) { PndTrkHit *hit = mvdstrhitlist->GetHit(ihit); double dist = hit->GetXYDistanceFromTrack(xc0, yc0, R0); // CHECK limits double distlimit = 0.5; if(dist < distlimit){ cluster.AddHit(hit); if(fDisplayOn) { TMarker *mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 25); mrk->SetMarkerColor(5); display->cd(1); mrk->Draw("SAME"); } } } for(int ihit = 0; ihit < conformalhitlist->GetNofHits(); ihit++) { PndTrkConformalHit *chit = conformalhitlist->GetHit(ihit); double dist = chit->GetDistanceFromTrack(fitm, fitq); int hitID = chit->GetHitID(); int detID = chit->GetDetectorID(); // CHECK limits if(detID != FairRootManager::Instance()->GetBranchId(fSttBranch)) continue; double distlimit = 0; if(chit->GetIsochrone() > 1.e-4) distlimit = 6 * chit->GetIsochrone(); else distlimit = 0.001; // --------------------------- if(dist < distlimit) { // cout << "ADDED " << hitID << " " << detID << " " << dist << " " << distlimit << endl; if(fDisplayOn) { TArc *arc = new TArc(chit->GetU(), chit->GetV(), chit->GetIsochrone()); arc->SetFillStyle(0); arc->SetLineColor(5); display->cd(2); arc->Draw("SAME"); } PndTrkHit *hit = stthitlist->GetHit(hitID); // ********** CHECK *********** // if you put "is used" the efficiency beecomes much lower! // DO NOT: if(!hit->IsUsed()) cluster.AddHit(hit); } } return cluster; } Bool_t PndTrkLegendreTask3::DoesRealHitBelong(PndTrkHit *hit, double x0, double y0, double R) { double dist = hit->GetXYDistanceFromTrack(x0, y0, R); int detID = hit->GetDetectorID(); if(detID == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) { if(dist < fMvdPix_RealDistLimit) return kTRUE; } else if(detID == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) { if(dist < fMvdStr_RealDistLimit) return kTRUE; } else if(detID == FairRootManager::Instance()->GetBranchId(fSttBranch)) { if(dist < fStt_RealDistLimit) return kTRUE; } return kFALSE; } Bool_t PndTrkLegendreTask3::DoesConfHitBelong(PndTrkConformalHit *chit, double fitm, double fitp) { double dist = chit->GetDistanceFromTrack(fitm, fitp); int detID = chit->GetDetectorID(); if(detID == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) { if(dist < fMvdPix_ConfDistLimit) return kTRUE; } else if(detID == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) { if(dist < fMvdStr_ConfDistLimit) return kTRUE; } else if(detID == FairRootManager::Instance()->GetBranchId(fSttBranch)) { if(chit->GetIsochrone() > 1.e-4) { if(dist < (6 * chit->GetIsochrone())) return kTRUE; } else if(dist < fStt_ConfDistLimit) return kTRUE; } return kFALSE; } PndTrkCluster PndTrkLegendreTask3::CreateClusterByDistance(Int_t mode, double fitm, double fitp) { PndTrkCluster cluster; // mode: ------------------------------------------- // 0 all conformal: pix conf - str conf - stt conf // 1 all real: pix real - str real - stt real // 2 mix: pix real - str real - stt conf // -------------------------------------------------- Double_t x0, y0, R; FromConformalToRealTrack(fitm, fitp, x0, y0, R); for(int ihit = 0; ihit < conformalhitlist->GetNofHits(); ihit++) { PndTrkConformalHit *chit = conformalhitlist->GetHit(ihit); PndTrkHit *hit = chit->GetHit(); Bool_t accept = kFALSE; if(hit == fRefHit) accept = DoesRealHitBelong(hit, x0, y0, R); else { if(mode == 0) accept = DoesConfHitBelong(chit, fitm, fitp); else if(mode == 1) accept = DoesRealHitBelong(hit, x0, y0, R); else if(mode == 2) { if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) accept = DoesConfHitBelong(chit, fitm, fitp); else accept = DoesRealHitBelong(hit, x0, y0, R); } } if(accept == kTRUE && fDisplayOn) { display->cd(1); if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) { TArc *arc = new TArc(hit->GetPosition().X(), hit->GetPosition().Y(), hit->GetIsochrone()); arc->SetFillStyle(0); arc->SetLineColor(5); arc->Draw("SAME"); } else { TMarker *mrk; if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 21); else if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 25); mrk->SetMarkerColor(5); mrk->Draw("SAME"); } } // ********** CHECK *********** // if you put "is used" the efficiency beecomes much lower! // DO NOT: if(!hit->IsUsed()) if(accept == kTRUE) { // cout << "ADDING TO CLUSTER" << endl; cluster.AddHit(hit); } } return cluster; } void PndTrkLegendreTask3::FromConformalToRealTrack(double fitm, double fitp, double &x0, double &y0, double &R) { // CHECK if this needs to be kept --> change xc0 to xc etc // center and radius Double_t xcrot0, ycrot0; ycrot0 = 1 / (2 * fitp); xcrot0 = - fitm * ycrot0; R = sqrt(xcrot0 * xcrot0 + ycrot0 * ycrot0); // re-rotation and re-traslation of xc and yc // rotation x0 = TMath::Cos(conformalhitlist->GetConformalTransform()->GetRotation())*xcrot0 - TMath::Sin(conformalhitlist->GetConformalTransform()->GetRotation())*ycrot0; y0 = TMath::Sin(conformalhitlist->GetConformalTransform()->GetRotation())*xcrot0 + TMath::Cos(conformalhitlist->GetConformalTransform()->GetRotation())*ycrot0; // traslation x0 += conformalhitlist->GetConformalTransform()->GetTranslation().X(); y0 += conformalhitlist->GetConformalTransform()->GetTranslation().Y(); } void PndTrkLegendreTask3::DrawGeometry() { if(hxy == NULL) hxy = new TH2F("hxy", "xy plane", 100, -43, 43, 100, -43, 43); else hxy->Reset(); display->cd(1); hxy->Draw(); display->Update(); display->Modified(); } void PndTrkLegendreTask3::DrawZGeometry(int whichone, double phimin, double phimax, double zmin, double zmax) { if(hxz == NULL) hxz = new TH2F("hxz", "xz plane", 100, -43, 43, 100, -43, 113); else hxz->Reset(); display->cd(3); hxz->Draw(); if(whichone == 2) { if(hzphi == NULL) hzphi = new TH2F("hzphi", "z - phi plane", phimax - phimin, phimin, phimax, zmax - zmin, zmin, zmax); else { hzphi->Reset(); hzphi->GetXaxis()->SetLimits(phimin, phimax); hzphi->GetYaxis()->SetLimits(zmin, zmax); } display->cd(4); hzphi->Draw(); } display->Update(); display->Modified(); } // ============================================================================================ // DISPLAY *********************************************************************************** // ============================================================================================ void PndTrkLegendreTask3::Refresh(){ // CHECK char goOnChar; // cout << "Refresh?" << endl; // cin >> goOnChar; // cout << "REFRESHING" << endl; DrawGeometry(); if(fVerbose) cout << "Refresh stt" << endl; DrawHits(stthitlist); if(fVerbose) cout << "Refresh pixel" << endl; DrawHits(mvdpixhitlist); if(fVerbose) cout << "Refresh strip" << endl; DrawHits(mvdstrhitlist); if(fVerbose) cout << "Refresh stop" << endl; } void PndTrkLegendreTask3::DrawHits(PndTrkHitList *hitlist) { display->cd(1); hitlist->Draw(); display->Update(); display->Modified(); } void PndTrkLegendreTask3::DrawGeometryConf(double x1, double y1, double x2, double y2) { // CHECK char goOnChar; // cout << "DRAWING GEOMETRY CONF" << endl; // cin >> goOnChar; // CHECK previous calculations, now not used; if(huv == NULL) huv = new TH2F("huv", "uv plane", 100, x1, y1, 100, x2, y2); else { huv->Reset(); huv->GetXaxis()->SetLimits(x1, y1); huv->GetYaxis()->SetLimits(x2, y2); } display->cd(2); huv->Draw(); display->Update(); display->Modified(); } void PndTrkLegendreTask3::DrawLegendreHisto() { display->cd(3); legendre->Draw(); display->Update(); display->Modified(); } void PndTrkLegendreTask3::RefreshConf() { // CHECK delete // CHECK char goOnChar; // cout << "RefreshConf?" << endl; // cin >> goOnChar; // cout << "REFRESHING CONF" << endl; } void PndTrkLegendreTask3::RefreshZ() { // CHECK delete // CHECK char goOnChar; // cout << "RefreshZ?" << endl; // cin >> goOnChar; // cout << "REFRESHING Z" << endl; } void PndTrkLegendreTask3::DrawConfHit(double u, double v, double r, int marker) { display->cd(2); if(r >= 0) { TArc *arc = new TArc(u, v, r); arc->SetFillStyle(0); arc->Draw("SAME"); } else { TMarker *mrk = new TMarker(u, v, marker); mrk->Draw("SAME"); } // display->Update(); // display->Modified(); } PndTrkCluster PndTrkLegendreTask3::CreateSttCluster(PndTrkHit *firsthit) { /** THE PROCEDURE ============= 2) LOOP (int hit2) on all the hits a) consider only the ones in region iregion 2) LOOP (int hit3) on all the hits IN THE CLUSTER a) IF the hit2 you are considering in the region IS the one in the cluster (hit3), continue. ELSE try to associate it to one of the hits in the cluster (IsSttAssociate) b) if IsSttAssociate is true, add the hit2 to the cluster too. **/ PndTrkCluster cluster; cluster.AddHit(firsthit); cluster.SetIRegion(firsthit->GetIRegion()); // 2) LOOP on all the hits for(int ihit = 0; ihit < stthitlist->GetNofHits(); ihit++) { PndTrkHit *hit = stthitlist->GetHit(ihit); if(hit->IsUsed()) { // cout << "already used IV" << endl; continue; } if(!hit->IsRegion(firsthit->GetIRegion())) continue; // cout << ihit << " " << cluster.GetNofHits() << endl; // 3) LOOP (int hit3) on all the hits IN THE CLUSTER for(int jhit = 0; jhit < cluster.GetNofHits(); jhit++) { PndTrkHit *clushit = cluster.GetHit(jhit); if(*hit == *clushit) { continue; } Bool_t success = IsSttAssociate(hit, clushit); // cout << "success " << success << endl; if(!success) continue; cluster.AddHit(hit); break; } } if(cluster.GetNofHits() < 2) for(int ihit = 0; ihit < cluster.GetNofHits(); ihit++) cluster.DeleteAllHits(); return cluster; } Bool_t PndTrkLegendreTask3::IsSttAssociate(PndTrkHit *hit1, PndTrkHit *hit2) { double distance = hit1->GetXYDistance(hit2); // ComputeDistance(hit1, hit2); if(distance < STTPARALDISTANCE) { // cout << "distance " << distance << endl; return kTRUE; } return kFALSE; } void PndTrkLegendreTask3::LightCluster(PndTrkCluster *cluster) { // cout << "LIGHT UP" << endl; Refresh(); cluster->LightUp(); display->Update(); display->Modified(); } // ================================== ZFINDER PndTrkCluster PndTrkLegendreTask3::CreateSkewHitList(PndTrkTrack *track) { double xc = track->GetCenter().X(); double yc = track->GetCenter().Y(); double R = track->GetRadius(); PndTrkCluster skewhitlist; double phimin = 400, phimax = -1, zmin = 1000, zmax = -1; for(int ihit = 0; ihit < stthitlist->GetNofHits(); ihit++) { PndTrkHit *hit = stthitlist->GetHit(ihit); // if(hit->IsUsed()) { // cout << "already used IV" << endl; // continue; // } if(hit->IsSttSkew() == kFALSE) continue; int tubeID = hit->GetTubeID(); PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID); TVector3 wireDirection = tube->GetWireDirection(); Double_t halflength = tube->GetHalfLength(); TVector3 first = tube->GetPosition() + wireDirection * halflength; // CHECK TVector3 second = tube->GetPosition() - wireDirection * halflength; // CHECK // if(fDisplayOn) { // char goOnChar; // display->cd(1); // TLine *l = new TLine(first.X(), first.Y(), second.X(), second.Y()); // l->SetLineColor(kBlue); // l->Draw("SAME"); // display->Update(); // display->Modified(); // cin >> goOnChar; // } double m1 = (first - second).Y()/(first - second).X(); double q1 = first.Y() - m1 * first.X(); // 1. compute intersection between the track circle and the wire TVector2 intersection1, intersection2; Int_t nofintersections = tools->ComputeSegmentCircleIntersection(TVector2(first.X(), first.Y()), TVector2(second.X(), second.Y()), xc, yc, R, intersection1, intersection2); if(nofintersections == 0) continue; if(nofintersections >= 2) { cout << "ERROR: MORE THAN 1 INTERSECTION!!" << endl; continue; // CHECK } if(fDisplayOn) { char goOnChar; display->cd(1); TLine *l = new TLine(first.X(), first.Y(), second.X(), second.Y()); l->SetLineColor(kBlue); l->Draw("SAME"); TMarker *mrk = new TMarker(intersection1.X(), intersection1.Y(), 20); mrk->SetMarkerColor(kBlue); mrk->Draw("SAME"); display->Update(); display->Modified(); // cin >> goOnChar; } // 2. find the tangent to the track in the intersection point // tangent approximation TVector2 tangent = tools->ComputeTangentInPoint(xc, yc, intersection1); // 3. rotate clockwise the tangent/point/(wire, not explicitely) // in order to have the wire parallel to the x axis; // then translate everything to have the wire ON the x axis double beta = wireDirection.Phi(); if(beta < 0) beta += TMath::Pi(); // ... rotate the tangent double rtx = TMath::Cos(beta) * tangent.X() + TMath::Sin(beta) * tangent.Y(); double rty = TMath::Cos(beta) * tangent.Y() - TMath::Sin(beta) * tangent.X(); TVector2 rottangent(rtx, rty); rottangent = rottangent.Unit(); // ... rotate the point double rx = TMath::Cos(beta) * intersection1.X() + TMath::Sin(beta) * intersection1.Y(); double ry = TMath::Cos(beta) * intersection1.Y() - TMath::Sin(beta) * intersection1.X(); // translation Double_t deltay = ry; rty -= deltay; ry -= deltay; // rotm, rotp Double_t rotm = rottangent.Y()/rottangent.X(); Double_t rotp = ry - rotm * rx; // ellipsis double a = hit->GetIsochrone() * TMath::Cos(SKEW_ANGLE); // CHECK skew angle hard coded double b = hit->GetIsochrone(); // center of ellipsis Double_t x0a, x0b, y0; y0 = 0.; x0a = (-rotp + TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm; x0b = (-rotp - TMath::Sqrt(b * b + a * a * rotm * rotm)) / rotm; // intersection point double intxa = (x0a * b * b - rotm * rotp * a * a) / (b * b + rotm * rotm * a * a); double intya = rotm * intxa + rotp; double intxb = (x0b * b * b - rotm * rotp * a * a) / (b * b + rotm * rotm * a * a); double intyb = rotm * intxb + rotp; // 4. retraslate/rerotate all back to the original plane // retranslate y0 += deltay; intya += deltay; intyb += deltay; // rerotate double x0anew = TMath::Cos(beta) * x0a - TMath::Sin(beta) * y0; double y0anew = TMath::Cos(beta) * y0 + TMath::Sin(beta) * x0a; double x0bnew = TMath::Cos(beta) * x0b - TMath::Sin(beta) * y0; double y0bnew = TMath::Cos(beta) * y0 + TMath::Sin(beta) * x0b; double intxanew = TMath::Cos(beta) * intxa - TMath::Sin(beta) * intya; double intyanew = TMath::Cos(beta) * intya + TMath::Sin(beta) * intxa; double intxbnew = TMath::Cos(beta) * intxb - TMath::Sin(beta) * intyb; double intybnew = TMath::Cos(beta) * intyb + TMath::Sin(beta) * intxb; intxa = intxanew; intya = intyanew; intxb = intxbnew; intyb = intybnew; // now we have x0a, y0a, center of the 1st ellipse // and x0b, y0b, center of the 2nd ellipse x0a = x0anew; double y0a = y0anew; x0b = x0bnew; double y0b = y0bnew; if(fDisplayOn) { char goOnChar; display->cd(1); TEllipse *ell1 = new TEllipse(x0a, y0a, a, b, 0, 360, -beta); ell1->SetFillStyle(0); ell1->SetLineColor(4); ell1->Draw("SAME"); TEllipse *ell2 = new TEllipse(x0b, y0b, a, b, 0, 360, -beta); ell2->SetFillStyle(0); ell2->SetLineColor(6); ell2->Draw("SAME"); TMarker *mrkinta = new TMarker(intxa, intya, 20); mrkinta->SetMarkerColor(4); mrkinta->Draw("SAME"); TMarker *mrkintb = new TMarker(intxb, intyb, 20); mrkintb->SetMarkerColor(6); mrkintb->Draw("SAME"); // cin >> goOnChar; } // 5. calculate z coordinate for each intersection // calculate z0a, z0b of the center of the ellipse Double_t t = ((x0a + y0a) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y())); Double_t z0a = first.Z() + (second.Z() - first.Z()) * t; // cout << "0 : calculate t, z0a " << t << " " << z0a << endl; t = ((x0b + y0b) - (first.X() + first.Y())) / ((second.X() - first.X()) + (second.Y() - first.Y())); Double_t z0b = first.Z() + (second.Z() - first.Z()) * t; TVector3 center1(x0a, y0a, z0a); TVector3 center2(x0b, y0b, z0b); if(fDisplayOn) { char goOnChar; display->cd(3); // cout << "COMPUTE Z COORDINATE" << endl; RefreshZ(); DrawZGeometry(); TLine *linezx = new TLine(first.X(), first.Z(), second.X(), second.Z()); linezx->Draw("SAME"); TMarker *mrkza = new TMarker(x0a, z0a, 20); mrkza->SetMarkerColor(4); mrkza->Draw("SAME"); TMarker *mrkzb = new TMarker(x0b, z0b, 20); mrkzb->SetMarkerColor(6); mrkzb->Draw("SAME"); // cin >> goOnChar; } // calculate the z of the intersection ON the ellipse (CHECK this step calculations!) double dx = intxa - x0a; double dy = intya - y0a; TVector3 dxdy(dx, dy, 0.0); TVector3 tfirst = first + dxdy; TVector3 tsecond = second + dxdy; t = ((intxa + intya) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y())); double intza = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t; if(fDisplayOn) { char goOnChar; display->cd(3); TLine *linezx1 = new TLine(tfirst.X(), tfirst.Z(), tsecond.X(), tsecond.Z()); linezx1->SetLineStyle(1); linezx1->Draw("SAME"); TMarker *mrkza1 = new TMarker(intxa, intza, 20); mrkza1->SetMarkerColor(kBlue - 9); mrkza1->Draw("SAME"); // cin >> goOnChar; } tfirst = first - dxdy; tsecond = second - dxdy; t = ((intxb + intyb) - (tfirst.X() + tfirst.Y())) / ((tsecond.X() - tfirst.X()) + (tsecond.Y() - tfirst.Y())); double intzb = tfirst.Z() + (tsecond.Z() - tfirst.Z()) * t; TVector3 fin_intersection1(intxa, intya, intza); TVector3 fin_intersection2(intxb, intyb, intzb); if(fDisplayOn) { char goOnChar; display->cd(3); TLine *linezx2 = new TLine(tfirst.X(), tfirst.Z(), tsecond.X(), tsecond.Z()); linezx2->SetLineStyle(1); linezx2->Draw("SAME"); TMarker *mrkzb1 = new TMarker(intxb, intzb, 20); mrkzb1->SetMarkerColor(kMagenta - 7); mrkzb1->Draw("SAME"); // cin >> goOnChar; } // int1.SetXYZ(intxa, intya, intza); // int2.SetXYZ(intxb, intyb, intzb); // errz = fabs(intza - intzb)/2. ; // CHECK to be changed int trackID = 1; double phi1 = track->ComputePhi(fin_intersection1); double phi2 = track->ComputePhi(fin_intersection2); PndTrkSkewHit *skewhit = new PndTrkSkewHit(*hit, trackID, center1, fin_intersection1, phi1, center2, fin_intersection2, phi2, a, b, -1, beta); // skewhit->Print(); skewhitlist.AddHit(skewhit); } // cout << "PHI: " << phimin << " " << phimax << endl; // cout << "Z : " << zmin << " " << zmax << endl; return skewhitlist; } // Bool_t PndTrkLegendreTask3::CheckPairOfHits(PndTrkHit *hit1, PndTrkHit *hit2) { // // cout << "hit1 " << hit1 << endl; // // cout << "hit2 " << hit2 << endl; // // both mvd // if(hit1->IsMvd() && hit2->IsMvd()) { // cout << "BOTH MVD " << hit1->GetDistance(hit2) << endl; // if(hit1->GetDistance(hit2) < 7) return kTRUE; // } // // one mvd - one stt // else if((hit1->IsMvd() && hit2->IsStt()) || (hit1->IsStt() && hit2->IsMvd())) { // cout << "ONE MVD/ONE STT " << hit1->GetXYDistance(hit2) << endl; // // if the stt one is skew -> WRONG for sure! // if((hit2->IsStt() && hit2->IsSttSkew()) || (hit1->IsStt() && hit1->IsSttSkew())) return kFALSE; // // if parallel: // else if(hit1->GetXYDistance(hit2) < 8) return kTRUE; // } // // both stt // else if(hit1->IsStt() && hit2->IsStt()) { // // one parallel one skew // if((hit1->IsSttParallel() && hit2->IsSttSkew()) || (hit2->IsSttParallel() && hit1->IsSttSkew())) { // if(hit1->GetXYDistance(hit2) < 3) return kTRUE; // } // // both parallel // else if(hit1->IsSttParallel() && hit2->IsSttParallel()) { // PndSttTube *tube1 = (PndSttTube*) fTubeArray->At(hit1->GetTubeID()); // int layerID1 = tube1->GetLayerID(); // PndSttTube *tube2 = (PndSttTube*) fTubeArray->At(hit2->GetTubeID()); // int layerID2 = tube2->GetLayerID(); // // // // cout << "layers " << layerID1 << " " << layerID2 << endl; // // same layer // if(layerID1 == tube2->GetLayerID()) return kTRUE; // // both inner/outer parallel layers // else if((layerID1 < 9 && layerID2 < 9) || (layerID1 > 8 && layerID2 > 8)) { // // adjacent layers // if(layerID1 == layerID2 - 1) return kTRUE; // else if(layerID1 == layerID2 + 1) return kTRUE; // } // // one inner(outer) one outer(inner) // else { // // must be on boundaries // CHECK isboundary // if(layerID1 == 7 && layerID2 == 16) return kTRUE; // else if(layerID1 == 16 && layerID2 == 7) return kTRUE; // } // } // // both skew // else if(hit1->IsSttSkew() && hit2->IsSttSkew()) { // PndSttTube *tube1 = (PndSttTube*) fTubeArray->At(hit1->GetTubeID()); // int layerID1 = layerID1; // PndSttTube *tube2 = (PndSttTube*) fTubeArray->At(hit2->GetTubeID()); // int layerID2 = layerID2; // // same layer // if(layerID1 == layerID2) return kTRUE; // // adjacent layers // else if(layerID1 == layerID2 - 1) return kTRUE; // else if(layerID1 == layerID2 + 1) return kTRUE; // } // } // return kFALSE; // } // void PndTrkLegendreTask3::Cleanup2(PndTrkCluster *cluster) { // // get most distance point // PndTrkHit *disthit = NULL; // double tmpdistance = -1; // for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { // PndTrkHit *hit = cluster->GetHit(ihit); // double dist = hit->GetXYDistance(TVector3(0., 0., 0.)); // if(dist > tmpdistance) { // disthit = hit; // tmpdistance = dist; // } // } // // sort // for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { // PndTrkHit *hit = cluster->GetHit(ihit); // double dist = hit->GetXYDistance(disthit); // hit->SetSortVariable(dist); // } // cluster->Sort(); // for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { // PndTrkHit *hit = cluster->GetHit(ihit); // // cout << "sorted hit " << ihit << " " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << hit->GetSortVariable() << " " << hit->IsSortable() << endl; // } // std::vector< int > failedhits, breakpoints; // for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { // PndTrkHit *hit = cluster->GetHit(ihit); // PndTrkHit *follow = cluster->GetNextHit(ihit); // PndTrkHit *previous = cluster->GetPreviousHit(ihit); // // cout << "pre " << previous << " fol " << follow << endl; // Bool_t pre = kTRUE, fol = kTRUE; // if(previous != NULL) pre = CheckPairOfHits(previous, hit); // if(follow != NULL) fol = CheckPairOfHits(follow, hit); // cout << "pre " << pre << " fol " << fol << endl; // // isolated hit wrong // if(pre == kFALSE && fol == kFALSE) { // cout << "isolated " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << pre << " " << fol << endl; // failedhits.push_back(ihit); // } // // first hit wrong // else if(previous == NULL && fol == kFALSE) { // PndTrkHit *follow2 = cluster->GetNextHit(ihit + 1); // if(CheckPairOfHits(follow, follow2) == kFALSE) { // failedhits.push_back(ihit); // cout << "first " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << pre << " " << fol << endl; // } // } // // last hit wrong // else if(follow == NULL && pre == kFALSE) { // PndTrkHit *previous2 = cluster->GetPreviousHit(ihit - 1); // if(CheckPairOfHits(previous, previous2) == kFALSE) { // cout << "last " << hit->GetDetectorID() << " " << hit->GetHitID() << " " << pre << " " << fol << endl; // failedhits.push_back(ihit); // } // } // // breakpoints // else if(pre == kTRUE && fol == kFALSE) { // if(find(failedhits.begin(), failedhits.end(), ihit - 1) == failedhits.end()) breakpoints.push_back(ihit); // } // } // cout << "CLEANUP PROCEDURE " << endl; // cout << "failedhits " << failedhits.size() << " breakpoints " << breakpoints.size() << endl; // if(fDisplayOn) { // char goOnChar; // display->cd(1); // for(int ihit = 0; ihit < failedhits.size(); ihit++) { // int hitno = failedhits[ihit]; // PndTrkHit *hit = cluster->GetHit(hitno); // hit->Draw(kMagenta); // } // for(int ihit = 0; ihit < breakpoints.size(); ihit++) { // int hitno = breakpoints[ihit]; // PndTrkHit *hit = cluster->GetHit(hitno); // hit->Draw(kGreen); // } // display->Update(); // display->Modified(); // } // } PndTrkCluster PndTrkLegendreTask3::CleanUpSkewHitList(PndTrkCluster *skewhitlist) { // 6. FIRST CLEANING of the track skewed associations // find most probable z - CHECK what happens if a track is very fwd peaked? // TRY WITH DELTA z = COST // cout << "##################### FIND MOST PROBABLE Z" << endl; double phimin = 400, phimax = -1, zmin = 1000, zmax = -1; for(int ihit = 0; ihit < skewhitlist->GetNofHits(); ihit++) { PndTrkSkewHit *skewhit = (PndTrkSkewHit*) skewhitlist->GetHit(ihit); if(!skewhit) continue; TVector3 fin_intersection1 = skewhit->GetIntersection1(); TVector3 fin_intersection2 = skewhit->GetIntersection2(); double phi1 = skewhit->GetPhi1(); double phi2 = skewhit->GetPhi2(); if(phi1 < phimin) phimin = phi1; if(phi2 < phimin) phimin = phi2; if(fin_intersection1.Z() < zmin) zmin = fin_intersection1.Z(); if(fin_intersection2.Z() < zmin) zmin = fin_intersection2.Z(); if(phi1 > phimax) phimax = phi1; if(phi2 > phimax) phimax = phi2; if(fin_intersection1.Z() > zmax) zmax = fin_intersection1.Z(); if(fin_intersection2.Z() > zmax) zmax = fin_intersection2.Z(); // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ CHECK } // cout << "PHI: " << phimin << " " << phimax << endl; // cout << "Z : " << zmin << " " << zmax << endl; // DISPLAY ---------------------------- if(fDisplayOn) { char goOnChar; display->cd(4); phimin -= 30; phimax += 30; zmin -= 13; zmax += 13; DrawZGeometry(2, phimin, phimax, zmin, zmax); hzphi->SetXTitle("#phi"); hzphi->SetYTitle("#z"); hzphi->Draw(); display->Update(); display->Modified(); for(int ihit = 0; ihit < skewhitlist->GetNofHits(); ihit++) { PndTrkSkewHit *skewhit = (PndTrkSkewHit*) skewhitlist->GetHit(ihit); if(!skewhit) continue; TVector3 fin_intersection1 = skewhit->GetIntersection1(); TVector3 fin_intersection2 = skewhit->GetIntersection2(); double phi1 = skewhit->GetPhi1(); double phi2 = skewhit->GetPhi2(); TLine *linezphi = new TLine(phi1, fin_intersection1.Z(), phi2, fin_intersection2.Z()); // TLine *linezphi = new TLine(fin_intersection1.Z(), phi1, fin_intersection2.Z(), phi2); linezphi->SetLineStyle(1); linezphi->Draw("SAME"); TMarker *mrkzphi1 = new TMarker(phi1, fin_intersection1.Z(), 20); // TMarker *mrkzphi1 = new TMarker(fin_intersection1.Z(), phi1, 20); mrkzphi1->SetMarkerColor(kBlue - 9); mrkzphi1->Draw("SAME"); TMarker *mrkzphi2 = new TMarker(phi2, fin_intersection2.Z(), 20); // TMarker *mrkzphi2 = new TMarker(fin_intersection2.Z(), phi2, 20); mrkzphi2->SetMarkerColor(kMagenta - 7); mrkzphi2->Draw("SAME"); } display->Update(); display->Modified(); // cin >> goOnChar; } // DISPLAY ---------------------------- TH1F hz("hz", "z", (zmax - zmin)/10., zmin, zmax); // CHECK for(int ihit = 0; ihit < skewhitlist->GetNofHits(); ihit++) { PndTrkSkewHit *skewhit = (PndTrkSkewHit*) skewhitlist->GetHit(ihit); if(!skewhit) continue; TVector3 fin_intersection1 = skewhit->GetIntersection1(); TVector3 fin_intersection2 = skewhit->GetIntersection2(); hz.Fill(fin_intersection1.Z()); hz.Fill(fin_intersection2.Z()); } int maxbinz = hz.GetMaximumBin(); double mostprobZ = hz.GetBinCenter(maxbinz); // cout << mostprobZ << endl; // delete hits too far away .......... // cout << "##################### DELETE HITS" << endl; PndTrkCluster tmpskewhitlist; for(int ihit = 0; ihit < skewhitlist->GetNofHits(); ihit++) { PndTrkSkewHit *skewhit = (PndTrkSkewHit*) skewhitlist->GetHit(ihit); if(!skewhit) continue; TVector3 fin_intersection1 = skewhit->GetIntersection1(); TVector3 fin_intersection2 = skewhit->GetIntersection2(); double phi1 = skewhit->GetPhi1(); double phi2 = skewhit->GetPhi2(); if(fabs(fin_intersection1.Z() - mostprobZ) > 30. && fabs(fin_intersection2.Z() - mostprobZ) > 30.) { // cout << "THROW AWAY " << ihit << " " << fabs(fin_intersection1.Z() - mostprobZ) << " " << fabs(fin_intersection2.Z() - mostprobZ) << endl; continue; } skewhit->SetSortVariable((phi1 + phi2)/2.); tmpskewhitlist.AddHit(skewhit); } tmpskewhitlist.Sort(); return tmpskewhitlist; } void PndTrkLegendreTask3::RegisterTrack(PndTrkTrack *track) { PndTrack *finaltrack = track->ConvertToPndTrack(); TClonesArray& clref1 = *fTrackArray; Int_t size = clref1.GetEntriesFast(); PndTrack *outputtrack = new(clref1[size]) PndTrack(finaltrack->GetParamFirst(), finaltrack->GetParamLast(), finaltrack->GetTrackCand()); TClonesArray& clref2 = *fTrackCandArray; size = clref2.GetEntriesFast(); PndTrackCand *outputtrackcand = new(clref2[size]) PndTrackCand(finaltrack->GetTrackCand()); if(fVerbose > 1) { cout << "MOM FIRST: TOT, PT, PL " << outputtrack->GetParamFirst().GetMomentum().Mag() << " " << outputtrack->GetParamFirst().GetMomentum().Perp() << " " << outputtrack->GetParamFirst().GetMomentum().Z() << endl; cout << "MOM LAST: TOT, PT, PL " << outputtrack->GetParamLast().GetMomentum().Mag() << " " << outputtrack->GetParamLast().GetMomentum().Perp() << " " << outputtrack->GetParamLast().GetMomentum().Z() << endl; } if(fDisplayOn) { char goOnChar; display->cd(1); Refresh(); track->Draw(kRed); display->Update(); display->Modified(); cout << "X, Y, R " << track->GetCenter().X() << " " << track->GetCenter().Y() << " " << track->GetRadius() << endl; cout << "Z0, TANL " << track->GetZ0() << " " << track->GetTanL() << endl; cout << "CHARGE " << track->GetCharge() << endl; // cin >> goOnChar; } } // ----------------------- FOR SECONDARIES PndTrkHit *PndTrkLegendreTask3::FindSttReferenceHit() { if(stthitlist->GetNofHits() == 0) return NULL; // loop on stt hits int tmphitid = -1; Double_t tmpiso = 1.; for(int jhit = 0; jhit < stthitlist->GetNofHits(); jhit++) { PndTrkHit *hit = stthitlist->GetHit(jhit); if(hit->IsUsed()) { if(fVerbose > 1) cout << "STT hit " << jhit << "already used " << endl; continue; } if(hit->IsSttSkew()) continue; if(hit->GetIsochrone() < tmpiso) { tmphitid = jhit; tmpiso = hit->GetIsochrone(); } } if(tmphitid == -1) return NULL; PndTrkHit *refhit = stthitlist->GetHit(tmphitid); if(fVerbose > 1) cout << "STT REFERENCE HIT " << tmphitid << " " << refhit->GetIsochrone() << endl; return refhit; } PndTrkHit *PndTrkLegendreTask3::FindMvdPixelReferenceHit() { if(mvdpixhitlist->GetNofHits() == 0) return NULL; // loop on mvd pix hits int tmphitid = -1; PndTrkHit *refhit = NULL; for(int jhit = 0; jhit < mvdpixhitlist->GetNofHits(); jhit++) { PndTrkHit *hit = mvdpixhitlist->GetHit(jhit); if(hit->IsUsed()) { if(fVerbose > 1) cout << "already used V" << endl; continue; } tmphitid = jhit; break; } if(tmphitid == -1) return NULL; refhit = mvdpixhitlist->GetHit(tmphitid); if(fVerbose > 1) cout << "MVD PIXEL REFERENCE HIT " << refhit->GetHitID() << endl; return refhit; } PndTrkHit *PndTrkLegendreTask3::FindMvdStripReferenceHit() { if(mvdstrhitlist->GetNofHits() == 0) return NULL; // loop on mvd str hits int tmphitid = -1; PndTrkHit *refhit = NULL; for(int jhit = 0; jhit < mvdstrhitlist->GetNofHits(); jhit++) { PndTrkHit *hit = mvdstrhitlist->GetHit(jhit); if(hit->IsUsed()) { if(fVerbose > 1) cout << "already used V" << endl; continue; } tmphitid = jhit; break; } if(tmphitid == -1) return NULL; refhit = mvdstrhitlist->GetHit(tmphitid); if(fVerbose > 1) cout << "MVD STRIP REFERENCE HIT " << refhit->GetHitID() << endl; return refhit; } PndTrkHit *PndTrkLegendreTask3::FindMvdReferenceHit() { PndTrkHit *refhit = NULL; refhit = FindMvdStripReferenceHit(); // refhit = FindMvdPixelReferenceHit(); if(refhit != NULL) return refhit; // refhit = FindMvdStripReferenceHit(); FindMvdPixelReferenceHit(); return refhit; } PndTrkHit *PndTrkLegendreTask3::FindReferenceHit() { PndTrkHit *refhit = NULL; // refhit = FindMvdReferenceHit(); refhit = FindSttReferenceHit(); if(refhit != NULL) return refhit; // refhit = FindSttReferenceHit(); FindMvdReferenceHit(); return refhit; } void PndTrkLegendreTask3::ComputeTraAndRot(PndTrkHit *hit, Double_t &delta, Double_t trasl[2]) { trasl[0] = hit->GetPosition().X(); trasl[1] = hit->GetPosition().Y(); delta = 0.; // TMath::ATan2(hit->GetPosition().Y() - 0., hit->GetPosition().X() - 0.); // CHECK } void PndTrkLegendreTask3::RePrepareLegendre(PndTrkCluster *cluster) { // cout << "RESETTING LEGENDRE HISTO" << endl; legendre->ResetLegendre(); if(fDisplayOn) { RefreshConf(); if(fSecondary) DrawGeometryConf(-1., 1., -1., 1.); else DrawGeometryConf(-0.07, 0.07, -0.07, 0.07); } // cout << "%%%%%%%%%%%%%%%%%%%% XY FINDE %%%%%%%%%%%%%%%%%%%%%%%%%%" << endl; FillLegendreHisto(cluster); } void PndTrkLegendreTask3::PrepareLegendre() { // cout << "RESETTING LEGENDRE HISTO" << endl; legendre->ResetLegendre(); if(fDisplayOn) { RefreshConf(); if(fSecondary) DrawGeometryConf(-1., 1., -1., 1.); else DrawGeometryConf(-0.07, 0.07, -0.07, 0.07); } if(fVerbose > 1) cout << "%%%%%%%%%%%%%%%%%%%% XY FINDER %%%%%%%%%%%%%%%%%%%%%%%%%%" << endl; FillLegendreHisto(0); } Int_t PndTrkLegendreTask3::ApplyLegendre(int icounter, double &theta_max, double &r_max) { PrepareLegendre(); return ExtractLegendre(icounter, 0, theta_max, r_max); } Int_t PndTrkLegendreTask3::ApplyLegendre(int icounter, PndTrkCluster *cluster, double &theta_max, double &r_max) { RePrepareLegendre(cluster); return ExtractLegendre(icounter, 1, theta_max, r_max); } Int_t PndTrkLegendreTask3::ExtractLegendre(int icounter, Int_t mode, double &theta_max, double &r_max) { if(fDisplayOn) { char goOnChar; // cin >> goOnChar; DrawLegendreHisto(); // cout << "LEGENDRE (nof conf hits = " << conformalhitlist->GetNofHits() << ")" << endl; display->cd(); // cin >> goOnChar; display->Update(); display->Modified(); } // FIND MAXIMUM IN LEGENDRE HISTO // legendre->ApplyThresholdLegendreHisto(0.3); cout << "icounter " << icounter << endl; PndTrkLegendreCluster thislist = legendre->ExtractLegendreMaximum(icounter); theta_max = thislist.GetTheta(); r_max = thislist.GetR(); int maxpeak = thislist.GetNofHits(); cout << "MAXPEAK " << maxpeak << endl; bool alreadythere = false; if(mode == 0) { // if(maxpeak <= 3) { // if(fVerbose > 1) cout << "MAXPEAK " << maxpeak << ", BREAK NOW! " << endl; // return maxpeak; // } for(int ialready = 0; ialready < fFoundPeaks.size(); ialready++) { std::pair foundthetar = fFoundPeaks.at(ialready); double foundtheta = foundthetar.first; double foundr = foundthetar.second; // IF THIS PEAK WAS ALREADY FOUND, DELETE THE PEAK AND GO ON (TO AVOID INFINITE LOOPS) if(theta_max == foundtheta && r_max == foundr) { // legendre->DeleteZoneAroundXYLegendre(theta_max, r_max); // maxpeak = legendre->ExtractLegendreMaximum(theta_max, r_max); alreadythere = true; if(fVerbose > 0) cout << "OH NO! THIS PEAK IS ALREADY THERE" << endl; return -1; } } if(alreadythere == false) { std::pair tr(theta_max, r_max); fFoundPeaks.push_back(tr); } } return maxpeak; } PndTrkCluster PndTrkLegendreTask3::Cleanup(PndTrkCluster cluster) { // CLEANUP // ------------------------------------------------------- // CleanCluster(cluster); // cluster->SortMvdPixelByLayerID(); // cluster->SortByLayerID(); int leftvotes = 0; int rightvotes = 0; TVector2 trasl = conform->GetTranslation(); Double_t angle = conform->GetRotation(); // SET SORT VARIABLE // MVD PIX PART OF CLUSTER PndTrkCluster mvdpixcluster = cluster.GetMvdPixelHitList(); for(int ihit = 0; ihit < mvdpixcluster.GetNofHits(); ihit++) { PndTrkHit *hit = mvdpixcluster.GetHit(ihit); if(fSecondary) hit->SetSortVariable(hit->GetDistance(TVector3(trasl.X(), trasl.Y(), angle))); else hit->SetSortVariable(hit->GetDistance(TVector3(0., 0., 0))); int sensorID = hit->GetSensorID(); int layerID = fClean->FindMvdLayer(sensorID); if(layerID%2 == 0) rightvotes++; else leftvotes++; } // MVD STR PART OF CLUSTER PndTrkCluster mvdstrcluster = cluster.GetMvdStripHitList(); for(int ihit = 0; ihit < mvdstrcluster.GetNofHits(); ihit++) { PndTrkHit *hit = mvdstrcluster.GetHit(ihit); if(fSecondary) hit->SetSortVariable(hit->GetDistance(TVector3(trasl.X(), trasl.Y(), angle))); else hit->SetSortVariable(hit->GetDistance(TVector3(0., 0., 0))); int sensorID = hit->GetSensorID(); int layerID = fClean->FindMvdLayer(sensorID); if(layerID%2 == 0) rightvotes++; else leftvotes++; } // STT PARALLEL PART OF CLUSTER PndTrkCluster sttcluster = cluster.GetSttParallelHitList(); for(int ihit = 0; ihit < sttcluster.GetNofHits(); ihit++) { PndTrkHit *hit = sttcluster.GetHit(ihit); int layerID = -999, tubeID = hit->GetTubeID(), sectorID = -999; PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID); layerID = tube->GetLayerID(); sectorID = tube->GetSectorID(); if(fSecondary) hit->SetSortVariable(hit->GetDistance(TVector3(trasl.X(), trasl.Y(), angle))); else hit->SetSortVariable(1000+layerID); if(sectorID < 3) rightvotes++; else leftvotes++; } // cout << "VOTES: LEFT " << leftvotes << " RIGHT " << rightvotes << endl; for(int ihit = 0; ihit < cluster.GetNofHits(); ihit++) { PndTrkHit *hit = cluster.GetHit(ihit); // cout << "for sorting " << hit->GetDetectorID () << " " << hit->GetHitID() << " " << hit->GetSortVariable() << endl; } cluster.Sort(); for(int ihit = 0; ihit < cluster.GetNofHits(); ihit++) { PndTrkHit *hit = cluster.GetHit(ihit); // cout << "sorted " << hit->GetDetectorID () << " " << hit->GetHitID() << " " << hit->GetSortVariable() << endl; } // cluster.Print(); // DELETE HITS PndTrkCluster deletioncluster; // LEFT/RIGHT CRITERION ------------------------------------------ NOW OFF // PIXEL mvdpixcluster = cluster.GetMvdPixelHitList(); for(int ihit = 0; ihit < mvdpixcluster.GetNofHits(); ihit++) { PndTrkHit *hit = mvdpixcluster.GetHit(ihit); int sensorID = hit->GetSensorID(); int layerID = fClean->FindMvdLayer(sensorID); // cout << "pixel sorted " << hit->GetHitID() << " " << hit->GetSortVariable() << " " << sensorID << " " << layerID << endl; if((leftvotes > rightvotes) && (layerID%2 == 0)) { // deletioncluster.AddHit(hit); // cout << "pixel DELETED " << hit->GetHitID() << " " << ihit << " " << mvdpixcluster.GetNofHits() << endl; } else if((leftvotes < rightvotes) && (layerID%2 != 0)) { // deletioncluster.AddHit(hit); // cout << "pixel DELETED " << hit->GetHitID() << " " << ihit << " " << mvdpixcluster.GetNofHits() << endl; } if(fDisplayOn) { char goOnChar; display->cd(1); hit->Draw(kGreen); display->Update(); display->Modified(); cout << "WAITING" << endl; cin >> goOnChar; } } // STRIP mvdstrcluster = cluster.GetMvdStripHitList(); for(int ihit = 0; ihit < mvdstrcluster.GetNofHits(); ihit++) { PndTrkHit *hit = mvdstrcluster.GetHit(ihit); int sensorID = hit->GetSensorID(); // cout << "strip sorted " << hit->GetHitID() << " " << hit->GetSortVariable() << " " << sensorID << " " << clean.FindMvdLayer(sensorID) << endl; int layerID = fClean->FindMvdLayer(sensorID); if((leftvotes > rightvotes) && (layerID%2 == 0)) { // deletioncluster.AddHit(hit); // cout << "strip DELETED " << hit->GetHitID() << endl; } else if((leftvotes < rightvotes) && (layerID%2 != 0)) { // deletioncluster.AddHit(hit); // cout << "strip DELETED " << hit->GetHitID() << endl; } if(fDisplayOn) { char goOnChar; display->cd(1); hit->Draw(kBlue); display->Update(); display->Modified(); // cin >> goOnChar; } } // STT PARALLEL sttcluster = cluster.GetSttParallelHitList(); for(int ihit = 0; ihit < sttcluster.GetNofHits(); ihit++) { PndTrkHit *hit = sttcluster.GetHit(ihit); int layerID = -999, tubeID = hit->GetTubeID(); PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID); layerID = tube->GetLayerID(); // cout << "stt sorted " << hit->GetHitID() << " " << hit->GetSortVariable() << " " << tubeID << " " << layerID << endl; int sectorID = tube->GetSectorID(); if((leftvotes > rightvotes) && (sectorID < 3)) { // deletioncluster.AddHit(hit); // cout << "stt DELETED " << hit->GetHitID() << " " << sectorID << endl; } else if((leftvotes < rightvotes) && (sectorID >= 3)) { // deletioncluster.AddHit(hit); // cout << "stt DELETED " << hit->GetHitID() << " " << sectorID << endl; } // ~~~~~~~~~~~~~~~~~~ check ******************* else { // NEIGHBORING CRITERION ------------------------------------------------ TArrayI neighboring = tube->GetNeighborings(); // cout << "NEIGHBORING TO TUBE " << tubeID << endl; // for(int nhit = 0; nhit < neighboring.GetSize(); nhit++) { // cout << " " << neighboring.At(nhit); // } // cout << endl; int neighboringcounter = 0; for(int ohit = 0; ohit < sttcluster.GetNofHits(); ohit++) { if(ihit == ohit) continue; PndTrkHit *otherhit = sttcluster.GetHit(ohit); int otherhitID = otherhit->GetTubeID(); for(int nhit = 0; nhit < neighboring.GetSize(); nhit++) { if(otherhitID == neighboring.At(nhit)) neighboringcounter++; } // if(isthere == false) // { // deletioncluster.AddHit(hit); // cout << "stt DELETEDbis " << hit->GetHitID() << " " << sectorID << endl; // } } // cout << "hit " << ihit << " has " << neighboringcounter << " NEIGHBORINS" << endl; if(layerID != 0 && layerID != 7 && layerID != 8 && ihit != sttcluster.GetNofHits() - 1 && neighboringcounter == 1) { // cout << "I SHOULD DELETE THIS" << endl; // deletioncluster.AddHit(hit); // OFF FOR NOW: delete if there is just 1 neighboring // cout << "stt DELETED " << hit->GetHitID() << " " << sectorID << endl; } // ON, delete if there is no neighboring (isolated hit) if(neighboringcounter == 0) { // cout << "I WILL DELETE THIS" << endl; deletioncluster.AddHit(hit); // cout << "stt DELETED " << hit->GetHitID() << " " << sectorID << endl; } if(fDisplayOn) { char goOnChar; display->cd(1); hit->Draw(kCyan); display->Update(); display->Modified(); // cin >> goOnChar; } } } if(fDisplayOn) { char goOnChar; // cin >> goOnChar; } // ~~~~~~~~~~~~~~~~~~ // ACTUAL DELETION OF HITS PndTrkCluster finalcluster; for(int ihit = 0; ihit < cluster.GetNofHits(); ihit++) { PndTrkHit *hit = cluster.GetHit(ihit); bool stop = false; for(int jhit = 0; jhit < deletioncluster.GetNofHits(); jhit++) { PndTrkHit *dhit = deletioncluster.GetHit(jhit); if(hit == dhit) { // ********** CHECK *********** // cout << "trhrow in again" << endl; // dhit->SetUsedFlag(0); stop = true; break; } } if(!stop) finalcluster.AddHit(hit); } return finalcluster; } PndTrkCluster* PndTrkLegendreTask3::CleanupZPhiFit(PndTrkCluster *cluster, double fitm, double fitp) { PndTrkCluster *cleancluster = new PndTrkCluster(); for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { PndTrkHit *hit = cluster->GetHit(ihit); if(hit->IsSttParallel()) { cleancluster->AddHit(hit); continue; // CHECK IsSttParrallel will be changed } TVector3 position = hit->GetPosition(); double phi = hit->GetPhi(); // d = | m x - y + p | / sqrt(m**2 + 1) double distance = fabs(fitm * phi - position.Z() + fitp )/sqrt(fitm * fitm + 1); // cout << "distance " << hit->GetHitID() << " " << distance << endl; if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch) || hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) { if(distance < 1) cleancluster->AddHit(hit); } else if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) { if(distance < 5) cleancluster->AddHit(hit); } } // cout << "cleancluster " << cleancluster->GetNofHits () << endl; return cleancluster; } Bool_t PndTrkLegendreTask3::ZPhiFit(int iter, PndTrkCluster *cluster, double &fitm, double &fitp) { if(iter != 0) { if(fDisplayOn) { char goOnChar; display->cd(3); DrawZGeometry(); hzphi->SetXTitle("#phi"); hzphi->SetYTitle("#z"); hzphi->Draw(); display->Update(); display->Modified(); } } fFitter->Reset(); PndTrkHit *hit = NULL; PndTrkHit *refhit = NULL; double refiso = 1000; for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { hit = cluster->GetHit(ihit); if(hit->IsSttParallel()) continue; // CHECK IsSttParrallel will be changed TVector3 position = hit->GetPosition(); double phi = hit->GetPhi(); // fFitter->SetPointToFit(phi, position.Z(), 0.1); // CHECK refhit stuff if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch) || hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) { if(iter != 0) refhit = hit; fFitter->SetPointToFit(phi, position.Z(), 0.001); } else if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) { if(iter != 0 && refhit == NULL && hit->GetIsochrone() < refiso) refhit = hit; fFitter->SetPointToFit(phi, position.Z(), 0.1); } if(fDisplayOn) { char goOnChar; if(iter == 0) display->cd(4); TMarker *mrkfoundzphi = NULL; if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdPixelBranch) || hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fMvdStripBranch)) { mrkfoundzphi = new TMarker(phi, position.Z(), 21); if(iter == 0) mrkfoundzphi->SetMarkerColor(kOrange); else mrkfoundzphi->SetMarkerColor(4); } else if(hit->GetDetectorID() == FairRootManager::Instance()->GetBranchId(fSttBranch)) { mrkfoundzphi = new TMarker(phi, position.Z(), 20); if(iter == 0) mrkfoundzphi->SetMarkerColor(kGreen); else mrkfoundzphi->SetMarkerColor(4); } mrkfoundzphi->Draw("SAME"); // display->cd(3); // mrkfoundzphi->Draw("SAME"); } } // =================================== bool fit = kFALSE; if(iter == 0) fit = fFitter->StraightLineFit(fitm, fitp); else { if(refhit != NULL) fit = fFitter->ConstrainedStraightLineFit(refhit->GetPhi(), refhit->GetPosition().Z(), fitm, fitp); else fit = fFitter->StraightLineFit(fitm, fitp); } if(fit) { if(fDisplayOn) { char goOnChar; if(iter != 0) display->cd(4); TLine *l22 = new TLine(-1000, -1000 * fitm + fitp, 1000, 1000 * fitm + fitp); if(iter == 0) l22->SetLineColor(3); else if(iter == 1) l22->SetLineColor(4); l22->Draw("SAME"); display->cd(3); l22->Draw("SAME"); display->Update(); display->Modified(); cin >> goOnChar; } } // ++++++++++++++++++++++++++++++++++ return fit; } double PndTrkLegendreTask3::ComputeZRediduals(PndTrkCluster *cluster, double fitm, double fitp) { TH1F hresiduals("hresiduals", "residuals in z", 20, -5, 5); for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { PndTrkHit *hit = cluster->GetHit(ihit); if(hit->GetDetectorID() != FairRootManager::Instance()->GetBranchId(fSttBranch)) continue; if(hit->IsSttParallel()) continue; TVector3 position = hit->GetPosition(); double phi = hit->GetPhi(); double distance = fitm * phi + fitp - position.Z(); hresiduals.Fill(distance); } if(fDisplayOn) { display->cd(4); char goOnChar; cout << "residuals"; cin >> goOnChar; hresiduals.Draw(); display->Update(); display->Modified(); cin >> goOnChar; } return hresiduals.GetMean(); } double PndTrkLegendreTask3::CorrectZ(PndTrkCluster *cluster, double deltaz, double fitm, double fitp) { for(int ihit = 0; ihit < cluster->GetNofHits(); ihit++) { PndTrkHit *hit = cluster->GetHit(ihit); if(hit->GetDetectorID() != FairRootManager::Instance()->GetBranchId(fSttBranch)) continue; if(hit->IsSttParallel()) continue; TVector3 position = hit->GetPosition(); double newz = position.Z() - deltaz; int tubeID = hit->GetTubeID(); PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID); TVector3 wireDirection = tube->GetWireDirection(); Double_t halflength = tube->GetHalfLength(); TVector3 first = tube->GetPosition() + wireDirection * halflength; // CHECK TVector3 second = tube->GetPosition() - wireDirection * halflength; // CHECK double m, p; double newx = position.X(); double newy = position.Y(); if((second.X() - first.X()) != 0) { m = (second.Z() - first.Z())/(second.X() - first.X()); p = position.Z() - m * position.X(); newx = (newz - p)/m; } if((second.Y() - first.Y()) != 0) { m = (second.Z() - first.Z())/(second.Y() - first.Y()); p = position.Z() - m * position.Y(); newy = (newz - p)/m; } hit->SetPosition(TVector3(newx, newy, newz)); if(fDisplayOn) { display->cd(1); char goOnChar; cout << endl; cout << "SEE IT" ; TMarker *mrk = new TMarker(newx, newy, 6); mrk->SetMarkerColor(kOrange); mrk->Draw("SAME"); display->Update(); display->Modified(); cout << "old " << position.X() << " " << position.Y() << " " << position.Z() << endl; cout << "new " << newx << " " << newy << " " << newz << endl; cin >> goOnChar; } } } void PndTrkLegendreTask3::HandleKey(Event_t *event) { // // Handle keyboard events. // char input[10]; // Int_t n; // UInt_t keysym; // if (event->fType == kGKeyPress) { // gVirtualX->LookupString(event, input, sizeof(input), keysym); // n = strlen(input); // if (event->fState & kKeyControlMask) { // Cntrl key modifier pressed // switch ((EKeySym)keysym & ~0x20) { // treat upper and lower the same // case kKey_C: // exit(0); // // fMenuFile->Activated(kQuitRoot); // default: // break; // } // } // } } ClassImp(PndTrkLegendreTask3)