Fork me on GitHub

Changes in / [91ef0b8:363e269] in git


Ignore:
Files:
1 deleted
2 edited

Legend:

Unmodified
Added
Removed

  • examples/ExamplePVtxFit.C

    r91ef0b8 r363e269  
    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 mm
    64                                 Double_t x = gp->X;
    65                                 Double_t y = gp->Y;
    66                                 Double_t z = gp->Z;
     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;
    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);
    8081                                        //
    8182                                        pr[Ntr] = new TVectorD(obsPar);
    82                                         //std::cout<<"Cov Matrix:"<<std::endl;
    83                                         //trk->CovarianceMatrix().Print();
    84                                         cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
     83                                        cv[Ntr] = new TMatrixDSym(TrkUtil::CovToMm(trk->CovarianceMatrix()));
    8584                                        Ntr++;
    8685                                }
    8786                        }               // 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 primary vertex pulls", 50, 50, 900, 500);
     120        TCanvas* Cnv = new TCanvas("Cnv", "Delphes generated track plots", 50, 50, 900, 500);
    121121        Cnv->Divide(2, 2);
    122122        Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);

  • external/TrackCovariance/VertexFit.cc

    r91ef0b8 r363e269  
    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
    4849}
    4950// ObsTrk list
    … …  
    5455        fVtxCst = kFALSE;
    5556        fxCst.ResizeTo(3);
    56         fCovCst.ResizeTo(3, 3);
     57        fCovCst.ResizeTo(3,3);
    5758        fXv.ResizeTo(3);
    58         fcovXv.ResizeTo(3, 3);
    59         //
    60         fPar = new TVectorD * [Ntr];
    61         fCov = new TMatrixDSym * [Ntr];
     59        fcovXv.ResizeTo(3,3);
     60        //
     61        fPar = new TVectorD*[Ntr];
     62        fCov = new TMatrixDSym*[Ntr];
    6263        fChi2List.ResizeTo(Ntr);
    6364        for (Int_t i = 0; i < Ntr; i++)
    6465        {
    6566                fPar[i] = new TVectorD(track[i]->GetObsPar());
     67                //std::cout << "fPar = " << std::endl; fPar[i]->Print();
    6668                fCov[i] = new TMatrixDSym(track[i]->GetCov());
     69                //std::cout << "fCov = " << std::endl; fCov[i]->Print();
    6770        }
    6871        //
    … …  
    7982        fWinvi = new TMatrixDSym * [Ntr];       // ACA'
    8083        for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3);
     84
     85        //std::cout << "VertexFit constructor executed" << std::endl;
    8186}
    8287//
    … …  
    9095        fChi2List.Clear();
    9196        //
    92         for (Int_t i = 0; i < fNtr; i++)
     97        for (Int_t i= 0; i < fNtr; i++)
    9398        {
    9499                fPar[i]->Clear();
    … …  
    143148        H(1, 1) = nn;
    144149        TVectorD c(2);
    145         c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i));
     150        c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) +=  n(i) * (y0(i) - x0(i));
    146151        c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i));
    147152        TVectorD smin = (H * c);
    … …  
    153158        Double_t R2 = TMath::Sqrt(Y(0) * Y(0) + Y(1) * Y(1));
    154159        //
    155         return 0.5 * (R1 + R2);
     160        return 0.5*(R1+R2);
    156161}
    157162Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2)
    … …  
    165170        TMatrixDSym H(2);
    166171        H(0, 0) = -TMath::Cos(p2(1));
    167         H(0, 1) = TMath::Cos(p1(1));
     172        H(0, 1) =  TMath::Cos(p1(1));
    168173        H(1, 0) = -TMath::Sin(p2(1));
    169         H(1, 1) = TMath::Sin(p1(1));
     174        H(1, 1) =  TMath::Sin(p1(1));
    170175        Double_t Det = TMath::Sin(p2(1) - p1(1));
    171176        H *= 1.0 / Det;
    … …  
    203208}
    204209
    205 TMatrixDSym VertexFit::RegInv3(TMatrixDSym& Smat0)
     210TMatrixDSym VertexFit::RegInv3(TMatrixDSym &Smat0)
    206211{
    207212        //
    … …  
    227232                        for (Int_t j = 0; j < 2; j++)Q(i, j) = RegMat(i, j);
    228233                }
    229                 Double_t Det = 1 - Q(0, 1) * Q(1, 0);
     234                Double_t Det = 1 - Q(0, 1)*Q(1, 0);
    230235                TMatrixDSym H(2);
    231236                H = Q;
    … …  
    236241                p(1) = RegMat(1, 2);
    237242                Double_t pHp = H.Similarity(p);
    238                 Double_t h = pHp - Det;
     243                Double_t h = pHp-Det;
    239244                //
    240245                TMatrixDSym pp(2); pp.Rank1Update(p);
    241                 TMatrixDSym F = (h * H) - pp.Similarity(H);
     246                TMatrixDSym F = (h*H) - pp.Similarity(H);
    242247                F *= 1.0 / Det;
    243                 TVectorD b = H * p;
     248                TVectorD b = H*p;
    244249                TMatrixDSym InvReg(3);
    245250                for (Int_t i = 0; i < 2; i++)
    … …  
    258263        else
    259264        {
    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);
     265                std::cout << "RegInv3: found negative elements in diagonal. Return standard inversion." << std::endl;
     266                return Smat.Invert();
    265267        }
    266268}
    … …  
    271273{
    272274        //
    273         // Derivative of track 3D position vector with respect to track parameters at constant phase 
     275        // Derivative of track 3D position vector with respect to track parameters at constant phase
    274276        //
    275277        // par = vector of track parameters
    … …  
    292294        A(2, 0) = 0.0;
    293295        // phi0
    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);
     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);
    296298        A(2, 1) = 0.0;
    297299        // 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);
     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);
    301303        // z0
    302304        A(0, 3) = 0.0;
    … …  
    351353        Double_t ct = par(4);
    352354        //
    353         x0(0) = -D * TMath::Sin(p0);
    354         x0(1) = D * TMath::Cos(p0);
     355        x0(0) = -D *TMath::Sin(p0);
     356        x0(1) = D*TMath::Cos(p0);
    355357        x0(2) = z0;
    356358        //
    … …  
    376378        x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
    377379        x(1) = x0(1) - (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
    378         x(2) = x0(2) + ct * phi / (2 * C);
     380        x(2) = x0(2) + ct*phi / (2 * C);
    379381        //
    380382        return x;
    … …  
    393395        //
    394396        // Stored quantities
    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'
     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'
    402404        //
    403405        // vertex radius approximation
    … …  
    413415        //
    414416        // Find track pair with largest phi difference
    415         Int_t isel = 0; Int_t jsel = 0; // selected track indices
     417        Int_t isel=0; Int_t jsel=0; // selected track indices
    416418        Double_t dphiMax = 0.0; // Max phi difference
    417 
    418         for (Int_t i = 0; i < fNtr - 1; i++)
     419        for (Int_t i = 0; i < fNtr-1; i++)
    419420        {
    420421                //ObsTrk* ti = tracks[i];
    … …  
    422423                Double_t phi1 = pari(1);
    423424
    424                 for (Int_t j = i + 1; j < fNtr; j++)
     425                for (Int_t j = i+1; j < fNtr; j++)
    425426                {
    426427                        //ObsTrk* tj = tracks[j];
    … …  
    437438        }
    438439        //
     440        //
     441        //ObsTrk* t1 = tracks[isel];
    439442        TVectorD p1 = *fPar[isel];
     443        //ObsTrk* t2 = tracks[jsel];
    440444        TVectorD p2 = *fPar[jsel];
    441445        Double_t R = FastRv1(p1, p2);
    442         if (R > 1000.0) R = Rd;
    443         R = 0.9 * R + 0.1 * Rd;
     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;
    444450        //
    445451        // Iteration properties
    … …  
    463469                {
    464470                        // Get track helix parameters and their covariance matrix
     471                        //ObsTrk *t = tracks[i];
    465472                        TVectorD par = *fPar[i];
    466473                        TMatrixDSym Cov = *fCov[i];
    467 
    468474                        Double_t fs;
    469475                        if (Ntry <= 0)  // Initialize all phases on first pass
    … …  
    471477                                Double_t D = par(0);
    472478                                Double_t C = par(2);
    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));
     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));
    475481                                fi[i] = fs;
    476482                        }
    … …  
    478484                        // Starting values
    479485                        //
    480                         fs = fi[i];             // Get phase
     486                        fs = fi[i];                                                             // Get phase
    481487                        //std::cout << "VertexFinder: phase fs set" << std::endl;
    482488                        TVectorD xs = Fill_x(par, fs);
    … …  
    501507                        // Build D matrix
    502508                        TMatrixDSym B(3);
    503 
    504509                        B.Rank1Update(a, 1.0);
    505510                        B *= -1. / a2;
    506511                        B.Similarity(W);
    507                         TMatrixDSym Ds = W + B;                                 // D matrix
     512                        TMatrixDSym Ds = W+B;                                   // D matrix
    508513                        Di[i] = new TMatrixDSym(Ds);                    // Store D matrix
    509514                        //std::cout << "VertexFinder: matrix Di stored" << std::endl;
    … …  
    518523                // update vertex position
    519524                TMatrixDSym H1 = RegInv3(H);
    520                 x = H1 * cterm;
     525                x = H1*cterm;
    521526                //std::cout << "VertexFinder: x vertex set" << std::endl;
    522527                // Update vertex covariance
    … …  
    527532                for (Int_t i = 0; i < fNtr; i++)
    528533                {
    529                         TVectorD lambda = (*Di[i]) * (*x0i[i] - x);
     534                        TVectorD lambda = (*Di[i])*(*x0i[i] - x);
    530535                        TMatrixDSym Wm1 = *Winvi[i];
    531536                        fChi2List(i) = Wm1.Similarity(lambda);
    532537                        Chi2 += fChi2List(i);
    533538                        TVectorD a = *ai[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];
     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];
    536541                }
    537542
     543                //std::cout << "VertexFinder: end chi2 calculation" << std::endl;
    538544                //
    539545                TVectorD dx = x - x0;
Note: See TracChangeset for help on using the changeset viewer.