//////////////////////////////////////////////////////////// // // PndTrkCombiLegendreTransform // // Class to perform the Combinatorial Legendre Transformation: // for each couple ij of hits, the 4 intersections of the r = r(theta) // sinusoidal curves are computed and the accumulation histo is // filled with these points only. // 1. the nof entries in the histo and its binning are not dependent anymore // 2. the nof entries in the histo is 2 * N * (N - 1), with N = nof hits // 3. if the couple ij was considered, then the couple ji is skipped // 4. this method is equivalent to fill the accumulation plot with the 4 tangents // of each couple of hits (4 lines are tangent to both two circles) // // authors: Lia Lavezzi - INFN Pavia (2012) // //////////////////////////////////////////////////////////// #include "PndTrkCombiLegendreTransform.h" #include using namespace std; // ----- Default constructor ------------------------------------------- PndTrkCombiLegendreTransform::PndTrkCombiLegendreTransform() { fThetaNofBin = 180; fThetaMin = 0; fThetaMax = 180; fRNofBin = 1000; fRMin = -0.07; fRMax = 0.07; fhLegendre = NULL; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndTrkCombiLegendreTransform::~PndTrkCombiLegendreTransform() { } // ------------------------------------------------------------------------- // ------ HISTOGRAM -------------------------------------------------------- void PndTrkCombiLegendreTransform::SetUpLegendreHisto() { if(fhLegendre == NULL) fhLegendre = new TH2F("fhCombiLegendre", "", fThetaNofBin, fThetaMin, fThetaMax, fRNofBin, fRMin, fRMax); else fhLegendre->SetBins(fThetaNofBin, fThetaMin, fThetaMax, fRNofBin, fRMin, fRMax); } void PndTrkCombiLegendreTransform::SetUpLegendreHisto(double thetaNofBin, double thetaMin, double thetaMax, double rNofBin, double rMin, double rMax) { fThetaNofBin = thetaNofBin; fThetaMin = thetaMin; fThetaMax = thetaMax; fRNofBin = rNofBin; fRMin = rMin; fRMax = rMax; SetUpLegendreHisto(); } void PndTrkCombiLegendreTransform::ResetLegendreHisto() { fhLegendre->Reset(); fhLegendre->SetStats(kFALSE); } void PndTrkCombiLegendreTransform::FillHisto(TH2F *histo, double x1, double y1, double r1, double x2, double y2, double r2) { double theta, r; ComputeThetaR(x1, y1, r1, x2, y2, r2, theta, r); histo->Fill(theta * TMath::RadToDeg(), r); } void PndTrkCombiLegendreTransform::ComputeThetaR(double x1, double y1, double r1, double x2, double y2, double r2, double &theta, double &r) { double dx = x1 - x2; double dy = y1 - y2; if(r1 == 0 && r2 == 0) { // two points theta = TMath::ATan2(-dx, dy); r = x1 * TMath::Cos(theta) + y1 * TMath::Sin(theta); // cout << "2points " << theta * TMath::RadToDeg() << "/ r " << r << endl; } else if(r1 > 0 && r2 > 0) { // two circles double dr[4] = {r1 - r2, r1 + r2, -r1 - r2, -r1 + r2}; for(int ir = 0; ir < 4; ir++) { double alpha = TMath::ATan2(dy, dx); if(dy < 0) alpha += (2 * TMath::Pi()); double distance = TMath::Sqrt(dx * dx + dy * dy); theta = alpha + TMath::ACos(-dr[ir]/distance); while(theta < 0) theta += TMath::Pi(); while(theta > TMath::Pi()) theta -= TMath::Pi(); // try all the combinations (++ +- -+ --) and... double rA[2], rB[2]; rA[0] = x1 * TMath::Cos(theta) + y1 * TMath::Sin(theta) - r1; rA[1] = x1 * TMath::Cos(theta) + y1 * TMath::Sin(theta) + r1; rB[0] = x2 * TMath::Cos(theta) + y2 * TMath::Sin(theta) - r2; rB[1] = x2 * TMath::Cos(theta) + y2 * TMath::Sin(theta) + r2; // ...get the one where curve A and curve B have the same r --> they cross if(fabs(rA[0] - rB[0]) < 1e-10 || fabs(rA[0] - rB[1]) < 1e-10) r = rA[0]; else if(fabs(rA[1] - rB[0]) < 1e-10 || fabs(rA[1] - rB[1]) < 1e-10) r = rA[1]; // cout << "2circles " << theta * TMath::RadToDeg() << "/ r " << r << endl; // cout << "ir " << ir << " theta " << theta * TMath::RadToDeg() << "/ r " << r << endl; } } else { // one point/one circle double rj; if(r1 > 0) rj = r1; else rj = r2; double dr[2] = {-rj, rj}; for(int ir = 0; ir < 2; ir++) { double alpha = TMath::ATan2(dy, dx); if(dy < 0) alpha += (2 * TMath::Pi()); double distance = TMath::Sqrt(dx * dx + dy * dy); theta = alpha + TMath::ACos(dr[ir]/distance); while(theta < 0) theta += TMath::Pi(); while(theta > TMath::Pi()) theta -= TMath::Pi(); if(r1 > 0) r = x2 * TMath::Cos(theta) + y2 * TMath::Sin(theta); else r = x1 * TMath::Cos(theta) + y1 * TMath::Sin(theta); // cout << "1pt/1cir " << theta * TMath::RadToDeg() << "/ r " << r << endl; // cout << "ir " << ir << " theta " << theta * TMath::RadToDeg() << "/ r " << r << endl; } } } void PndTrkCombiLegendreTransform::FillLegendreHisto(double x1, double y1, double radius1, double x2, double y2, double radius2) { FillHisto(fhLegendre, x1, y1, radius1, x2, y2, radius2); } void PndTrkCombiLegendreTransform::Draw() { fhLegendre->Draw("colz"); } int PndTrkCombiLegendreTransform::ExtractLegendreMaximum(double &theta_max, double &r_max) { ExtractMaximumFromHisto(fhLegendre, theta_max, r_max); } int PndTrkCombiLegendreTransform::ExtractMaximumFromHisto(TH2F *histo, double &theta_max, double &r_max) { int bin = histo->GetMaximumBin(); int binx, biny, binz; histo->GetBinXYZ(bin, binx, biny, binz); theta_max = histo->GetXaxis()->GetBinCenter(binx); r_max = histo->GetYaxis()->GetBinCenter(biny); return histo->GetMaximum(); } void PndTrkCombiLegendreTransform::ExtractLegendreSingleLineParameters(double &slope, double &intercept) { double theta, r; ExtractLegendreMaximum(theta, r); ExtractLineParameters(theta, r, slope, intercept); } void PndTrkCombiLegendreTransform::ExtractLineParameters(double theta, double r, double &slope, double &intercept) { // r = x cost + y sint --> y = r/sint - x/tant if(theta == 0.) theta = 1.e-6; // CHECK slope = - 1./TMath::Tan(TMath::DegToRad() * theta); intercept = r/TMath::Sin(TMath::DegToRad() * theta); } ClassImp(PndTrkCombiLegendreTransform)