Fork me on GitHub

Changeset bd33fce in git


Ignore:
Timestamp:
Mar 8, 2021, 9:52:15 AM (4 years ago)
Author:
GitHub <noreply@…>
Branches:
master
Children:
a10293e, b1cb322
Parents:
2eaa9fd (diff), 91ef0b8 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Michele Selvaggi <michele.selvaggi@…> (03/08/21 09:52:15)
git-committer:
GitHub <noreply@…> (03/08/21 09:52:15)
Message:

Merge pull request #88 from fbedesch/master

Fixed conflicts and Reginv3

Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • examples/ExamplePVtxFit.C

    r2eaa9fd rbd33fce  
    4040        //
    4141        // Loop over all events
    42         Int_t Nev = TMath::Min(Nevent, (Int_t) numberOfEntries);
     42        Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
    4343        for (Int_t entry = 0; entry < Nev; ++entry)
    4444        {
    … …  
    6161                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    6262                                //
    63                                 // Position of origin in meters
    64                                 Double_t x = 1.0e-3 * gp->X;
    65                                 Double_t y = 1.0e-3 * gp->Y;
    66                                 Double_t z = 1.0e-3 * gp->Z;
     63                                // Position of origin in mm
     64                                Double_t x = gp->X;
     65                                Double_t y = gp->Y;
     66                                Double_t z = gp->Z;
    6767                                //
    6868                                // group tracks originating from the primary vertex
    … …  
    7878                                        Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
    7979                                        TVectorD obsPar(5, oPar);       // Fill observed parameters
    80                                         TVector3 xv(x, y, z);
    8180                                        //
    8281                                        pr[Ntr] = new TVectorD(obsPar);
    83                                         cv[Ntr] = new TMatrixDSym(TrkUtil::CovToMm(trk->CovarianceMatrix()));
     82                                        //std::cout<<"Cov Matrix:"<<std::endl;
     83                                        //trk->CovarianceMatrix().Print();
     84                                        cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
    8485                                        Ntr++;
    8586                                }
    8687                        }               // End loop on tracks
    87                         //std::cout << "Total of " << Ntr << " primary tracks out of " << NtrG << " tracks" << std::endl;
    8888                }
    8989                //
    … …  
    118118        // Show resulting histograms
    119119        //
    120         TCanvas* Cnv = new TCanvas("Cnv", "Delphes generated track plots", 50, 50, 900, 500);
     120        TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
    121121        Cnv->Divide(2, 2);
    122122        Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);

  • external/TrackCovariance/VertexFit.cc

    r2eaa9fd rbd33fce  
    1515        fVtxCst = kFALSE;
    1616        fxCst.ResizeTo(3);
    17         fCovCst.ResizeTo(3,3);
     17        fCovCst.ResizeTo(3, 3);
    1818        fXv.ResizeTo(3);
    19         fcovXv.ResizeTo(3,3);
     19        fcovXv.ResizeTo(3, 3);
    2020}
    2121// Parameters and covariances
    … …  
    2626        fVtxCst = kFALSE;
    2727        fxCst.ResizeTo(3);
    28         fCovCst.ResizeTo(3,3);
     28        fCovCst.ResizeTo(3, 3);
    2929        fXv.ResizeTo(3);
    30         fcovXv.ResizeTo(3,3);
     30        fcovXv.ResizeTo(3, 3);
    3131        //
    3232        fPar = trkPar;
    … …  
    3535        //
    3636        ffi = new Double_t[Ntr];                                // Fit phases
    37         fx0i = new TVectorD* [Ntr];                     // Track expansion points
     37        fx0i = new TVectorD * [Ntr];                    // Track expansion points
    3838        for (Int_t i = 0; i < Ntr; i++) fx0i[i] = new TVectorD(3);
    3939        fai = new TVectorD * [Ntr];                     // dx/dphi
    4040        for (Int_t i = 0; i < Ntr; i++) fai[i] = new TVectorD(3);
    4141        fa2i = new Double_t[Ntr];                               // a'Wa
    42         fDi = new TMatrixDSym*[Ntr];            // W-WBW
     42        fDi = new TMatrixDSym * [Ntr];          // W-WBW
    4343        for (Int_t i = 0; i < Ntr; i++) fDi[i] = new TMatrixDSym(3);
    4444        fWi = new TMatrixDSym * [Ntr];  // (ACA')^-1
    … …  
    4646        fWinvi = new TMatrixDSym * [Ntr];       // ACA'
    4747        for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3);
    48 
    4948}
    5049// ObsTrk list
    … …  
    5554        fVtxCst = kFALSE;
    5655        fxCst.ResizeTo(3);
    57         fCovCst.ResizeTo(3,3);
     56        fCovCst.ResizeTo(3, 3);
    5857        fXv.ResizeTo(3);
    59         fcovXv.ResizeTo(3,3);
    60         //
    61         fPar = new TVectorD*[Ntr];
    62         fCov = new TMatrixDSym*[Ntr];
     58        fcovXv.ResizeTo(3, 3);
     59        //
     60        fPar = new TVectorD * [Ntr];
     61        fCov = new TMatrixDSym * [Ntr];
    6362        fChi2List.ResizeTo(Ntr);
    6463        for (Int_t i = 0; i < Ntr; i++)
    6564        {
    6665                fPar[i] = new TVectorD(track[i]->GetObsPar());
    67                 //std::cout << "fPar = " << std::endl; fPar[i]->Print();
    6866                fCov[i] = new TMatrixDSym(track[i]->GetCov());
    69                 //std::cout << "fCov = " << std::endl; fCov[i]->Print();
    7067        }
    7168        //
    … …  
    8279        fWinvi = new TMatrixDSym * [Ntr];       // ACA'
    8380        for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3);
    84 
    85         //std::cout << "VertexFit constructor executed" << std::endl;
    8681}
    8782//
    … …  
    9590        fChi2List.Clear();
    9691        //
    97         for (Int_t i= 0; i < fNtr; i++)
     92        for (Int_t i = 0; i < fNtr; i++)
    9893        {
    9994                fPar[i]->Clear();
    … …  
    148143        H(1, 1) = nn;
    149144        TVectorD c(2);
    150         c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) +=  n(i) * (y0(i) - x0(i));
     145        c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i));
    151146        c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i));
    152147        TVectorD smin = (H * c);
    … …  
    158153        Double_t R2 = TMath::Sqrt(Y(0) * Y(0) + Y(1) * Y(1));
    159154        //
    160         return 0.5*(R1+R2);
     155        return 0.5 * (R1 + R2);
    161156}
    162157Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2)
    … …  
    170165        TMatrixDSym H(2);
    171166        H(0, 0) = -TMath::Cos(p2(1));
    172         H(0, 1) =  TMath::Cos(p1(1));
     167        H(0, 1) = TMath::Cos(p1(1));
    173168        H(1, 0) = -TMath::Sin(p2(1));
    174         H(1, 1) =  TMath::Sin(p1(1));
     169        H(1, 1) = TMath::Sin(p1(1));
    175170        Double_t Det = TMath::Sin(p2(1) - p1(1));
    176171        H *= 1.0 / Det;
    … …  
    208203}
    209204
    210 TMatrixDSym VertexFit::RegInv3(TMatrixDSym &Smat0)
     205TMatrixDSym VertexFit::RegInv3(TMatrixDSym& Smat0)
    211206{
    212207        //
    … …  
    232227                        for (Int_t j = 0; j < 2; j++)Q(i, j) = RegMat(i, j);
    233228                }
    234                 Double_t Det = 1 - Q(0, 1)*Q(1, 0);
     229                Double_t Det = 1 - Q(0, 1) * Q(1, 0);
    235230                TMatrixDSym H(2);
    236231                H = Q;
    … …  
    241236                p(1) = RegMat(1, 2);
    242237                Double_t pHp = H.Similarity(p);
    243                 Double_t h = pHp-Det;
     238                Double_t h = pHp - Det;
    244239                //
    245240                TMatrixDSym pp(2); pp.Rank1Update(p);
    246                 TMatrixDSym F = (h*H) - pp.Similarity(H);
     241                TMatrixDSym F = (h * H) - pp.Similarity(H);
    247242                F *= 1.0 / Det;
    248                 TVectorD b = H*p;
     243                TVectorD b = H * p;
    249244                TMatrixDSym InvReg(3);
    250245                for (Int_t i = 0; i < 2; i++)
    … …  
    263258        else
    264259        {
    265                 std::cout << "RegInv3: found negative elements in diagonal. Return standard inversion." << std::endl;
    266                 return Smat.Invert();
     260                D.Zero();
     261                for (Int_t i = 0; i < N; i++) D(i, i) = 1.0 / TMath::Sqrt(TMath::Abs(Smat(i, i)));
     262                TMatrixDSym RegMat = Smat.Similarity(D);
     263                RegMat.Invert();
     264                return RegMat.Similarity(D);
    267265        }
    268266}
    … …  
    273271{
    274272        //
    275         // Derivative of track 3D position vector with respect to track parameters at constant phase
     273        // Derivative of track 3D position vector with respect to track parameters at constant phase 
    276274        //
    277275        // par = vector of track parameters
    … …  
    294292        A(2, 0) = 0.0;
    295293        // phi0
    296         A(0, 1) = -D*TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
    297         A(1, 1) = -D*TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
     294        A(0, 1) = -D * TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
     295        A(1, 1) = -D * TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
    298296        A(2, 1) = 0.0;
    299297        // C
    300         A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C*C);
    301         A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C*C);
    302         A(2, 2) = -ct*phi / (2 * C*C);
     298        A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C * C);
     299        A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C * C);
     300        A(2, 2) = -ct * phi / (2 * C * C);
    303301        // z0
    304302        A(0, 3) = 0.0;
    … …  
    353351        Double_t ct = par(4);
    354352        //
    355         x0(0) = -D *TMath::Sin(p0);
    356         x0(1) = D*TMath::Cos(p0);
     353        x0(0) = -D * TMath::Sin(p0);
     354        x0(1) = D * TMath::Cos(p0);
    357355        x0(2) = z0;
    358356        //
    … …  
    378376        x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
    379377        x(1) = x0(1) - (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
    380         x(2) = x0(2) + ct*phi / (2 * C);
     378        x(2) = x0(2) + ct * phi / (2 * C);
    381379        //
    382380        return x;
    … …  
    395393        //
    396394        // Stored quantities
    397         Double_t *fi = new Double_t[fNtr];                              // Phases
    398         TVectorD **x0i = new TVectorD*[fNtr];                   // Track expansion point
    399         TVectorD **ai = new TVectorD*[fNtr];                            // dx/dphi
    400         Double_t *a2i = new Double_t[fNtr];                             // a'Wa
    401         TMatrixDSym **Di = new TMatrixDSym*[fNtr];              // W-WBW
    402         TMatrixDSym **Wi = new TMatrixDSym*[fNtr];              // (ACA')^-1
    403         TMatrixDSym **Winvi = new TMatrixDSym*[fNtr];   // ACA'
     395        Double_t* fi = new Double_t[fNtr];                              // Phases
     396        TVectorD** x0i = new TVectorD * [fNtr];                 // Track expansion point
     397        TVectorD** ai = new TVectorD * [fNtr];                          // dx/dphi
     398        Double_t* a2i = new Double_t[fNtr];                             // a'Wa
     399        TMatrixDSym** Di = new TMatrixDSym * [fNtr];            // W-WBW
     400        TMatrixDSym** Wi = new TMatrixDSym * [fNtr];            // (ACA')^-1
     401        TMatrixDSym** Winvi = new TMatrixDSym * [fNtr]; // ACA'
    404402        //
    405403        // vertex radius approximation
    … …  
    415413        //
    416414        // Find track pair with largest phi difference
    417         Int_t isel=0; Int_t jsel=0; // selected track indices
     415        Int_t isel = 0; Int_t jsel = 0; // selected track indices
    418416        Double_t dphiMax = 0.0; // Max phi difference
    419         for (Int_t i = 0; i < fNtr-1; i++)
     417
     418        for (Int_t i = 0; i < fNtr - 1; i++)
    420419        {
    421420                //ObsTrk* ti = tracks[i];
    … …  
    423422                Double_t phi1 = pari(1);
    424423
    425                 for (Int_t j = i+1; j < fNtr; j++)
     424                for (Int_t j = i + 1; j < fNtr; j++)
    426425                {
    427426                        //ObsTrk* tj = tracks[j];
    … …  
    438437        }
    439438        //
    440         //
    441         //ObsTrk* t1 = tracks[isel];
    442439        TVectorD p1 = *fPar[isel];
    443         //ObsTrk* t2 = tracks[jsel];
    444440        TVectorD p2 = *fPar[jsel];
    445441        Double_t R = FastRv1(p1, p2);
    446         if (R > 1.0) R = Rd;
    447         R = 0.9*R + 0.1*Rd;
    448         //
    449         //std::cout << "VertexFinder: fast R = " << R << std::endl;
     442        if (R > 1000.0) R = Rd;
     443        R = 0.9 * R + 0.1 * Rd;
    450444        //
    451445        // Iteration properties
    … …  
    469463                {
    470464                        // Get track helix parameters and their covariance matrix
    471                         //ObsTrk *t = tracks[i];
    472465                        TVectorD par = *fPar[i];
    473466                        TMatrixDSym Cov = *fCov[i];
     467
    474468                        Double_t fs;
    475469                        if (Ntry <= 0)  // Initialize all phases on first pass
    … …  
    477471                                Double_t D = par(0);
    478472                                Double_t C = par(2);
    479                                 Double_t arg = TMath::Max(1.0e-6, (R*R - D*D) / (1 + 2 * C*D));
    480                                 fs = 2 * TMath::ASin(C*TMath::Sqrt(arg));
     473                                Double_t arg = TMath::Max(1.0e-6, (R * R - D * D) / (1 + 2 * C * D));
     474                                fs = 2 * TMath::ASin(C * TMath::Sqrt(arg));
    481475                                fi[i] = fs;
    482476                        }
    … …  
    484478                        // Starting values
    485479                        //
    486                         fs = fi[i];                                                             // Get phase
     480                        fs = fi[i];             // Get phase
    487481                        //std::cout << "VertexFinder: phase fs set" << std::endl;
    488482                        TVectorD xs = Fill_x(par, fs);
    … …  
    507501                        // Build D matrix
    508502                        TMatrixDSym B(3);
     503
    509504                        B.Rank1Update(a, 1.0);
    510505                        B *= -1. / a2;
    511506                        B.Similarity(W);
    512                         TMatrixDSym Ds = W+B;                                   // D matrix
     507                        TMatrixDSym Ds = W + B;                                 // D matrix
    513508                        Di[i] = new TMatrixDSym(Ds);                    // Store D matrix
    514509                        //std::cout << "VertexFinder: matrix Di stored" << std::endl;
    … …  
    523518                // update vertex position
    524519                TMatrixDSym H1 = RegInv3(H);
    525                 x = H1*cterm;
     520                x = H1 * cterm;
    526521                //std::cout << "VertexFinder: x vertex set" << std::endl;
    527522                // Update vertex covariance
    … …  
    532527                for (Int_t i = 0; i < fNtr; i++)
    533528                {
    534                         TVectorD lambda = (*Di[i])*(*x0i[i] - x);
     529                        TVectorD lambda = (*Di[i]) * (*x0i[i] - x);
    535530                        TMatrixDSym Wm1 = *Winvi[i];
    536531                        fChi2List(i) = Wm1.Similarity(lambda);
    537532                        Chi2 += fChi2List(i);
    538533                        TVectorD a = *ai[i];
    539                         TVectorD b = (*Wi[i])*(x - *x0i[i]);
    540                         for (Int_t j = 0; j < 3; j++)fi[i] += a(j)*b(j) / a2i[i];
     534                        TVectorD b = (*Wi[i]) * (x - *x0i[i]);
     535                        for (Int_t j = 0; j < 3; j++)fi[i] += a(j) * b(j) / a2i[i];
    541536                }
    542537
    543                 //std::cout << "VertexFinder: end chi2 calculation" << std::endl;
    544538                //
    545539                TVectorD dx = x - x0;
Note: See TracChangeset for help on using the changeset viewer.