Fork me on GitHub

Changes in / [0c0c9af:78ce8d1] in git


Ignore:
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • examples/ExamplePVtxFind.C

    r0c0c9af r78ce8d1  
    11/*
    22Example of using vertex fitter class to fit primary vertex
    3 assumed to be generated in (0,0,0)
     3assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
    44*/
    55
     
    4444        TH1D* hTrFoundG = new TH1D("hTrFoundG", "Found GOOD primary tracks", 41, -0.5, 40.5);
    4545        TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD  primary tracks", 41, -0.5, 40.5);
     46        TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);
    4647        //
    4748        // Loop over all events
     
    8788                                pr[it] = new TVectorD(obsPar);
    8889                                cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
     90                                // Check covariance
     91                                /*
     92                                TMatrixDSymEigen Eign(*cv[it]);
     93                                TVectorD lambda = Eign.GetEigenValues();
     94                                Bool_t good = kTRUE;
     95                                for(Int_t i=0; i<5; i++) if(lambda(i)<0.0) good = kFALSE;
     96                                if(!good){
     97                                        std::cout<<"Cov["<<it<<"]"<<" eigenvalues: "
     98                                        <<lambda(0)<<", "<<lambda(1)<<", "<<lambda(2)<<", "
     99                                        <<lambda(3)<<", "<<lambda(4)<<std::endl;
     100                                        TMatrixDSym Ctmp = *cv[it];
     101                                        std::cout<<"Error D (dir - cov) "<<trk->ErrorD0
     102                                        <<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl;
     103                                }
     104                                */
    89105                                //std::cout << "Track loaded it= " << it << std::endl;
    90106                        //
     
    120136                        }               // End loop on tracks
    121137                }
     138                std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
    122139                //
    123140                // Find primary vertex
     
    132149                covpvc(1, 1) = 2.55e-05 * 2.55e-05;
    133150                covpvc(2, 2) = 0.64 * 0.64;
     151                if(Nprim == 0){
     152                        xpv = xpvc(0);
     153                        ypv = xpvc(1);
     154                        zpv = xpvc(2);
     155                }
    134156                //
    135157                //
     
    137159                Int_t nSkim = 0;
    138160                Int_t* nSkimmed = new Int_t[NtrG];
    139                 VertexFit* Vskim = new VertexFit();
    140                 Vskim->AddVtxConstraint(xpvc, covpvc);
    141                 Double_t MaxChi2 = 2.0+3*sqrt(2.);
     161                TVectorD** PrSk = new TVectorD * [1];
     162                TMatrixDSym** CvSk = new TMatrixDSym * [1];
     163                Double_t MaxChi2 = 2.0+3*2.;
    142164                for (Int_t n = 0; n < NtrG; n++) {
    143                         Vskim->AddTrk(pr[n], cv[n]);
     165                        PrSk[0] = new TVectorD(*pr[n]);
     166                        CvSk[0] = new TMatrixDSym(*cv[n]);
     167                        VertexFit* Vskim = new VertexFit(1,PrSk, CvSk);
     168                        Vskim->AddVtxConstraint(xpvc, covpvc);
    144169                        Double_t Chi2One = Vskim->GetVtxChi2();
    145170                        //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
     
    149174                                nSkim++;
    150175                        }
    151                         Vskim->RemoveTrk(0);
    152                 }
     176                        // Cleanup
     177                        delete Vskim;
     178                }
     179                delete PrSk[0];
     180                delete CvSk[0];
     181                delete[] PrSk;
     182                delete[] CvSk;
     183                //
    153184                // Load tracks for primary fit
    154185                Int_t MinTrk = 1;       // Minumum # tracks for vertex fit
    155                 TVectorD** PrFit = new TVectorD * [nSkim];
    156                 TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
    157186                if (nSkim >= MinTrk) {
     187                        TVectorD** PrFit = new TVectorD * [nSkim];
     188                        TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
    158189                        for (Int_t n = 0; n < nSkim; n++) {
    159190                                PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
    160191                                CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
    161                                 TVectorD test = *PrFit[n];
    162                                 //std::cout << "Test *PrFit[" << n << "](0) = " << test(0) << std::endl;
    163                         }
    164                 }
    165                 delete[] nSkimmed;
    166                 delete Vskim;
    167                 Int_t Nfound = nSkim;
    168                 if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
    169                 if (nSkim >= MinTrk) {
     192                        }
     193                        delete[] nSkimmed;
     194                        Int_t Nfound = nSkim;
     195                        if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
     196                        const Int_t MaxFound = 100; Double_t Chi2LL[MaxFound]; Double_t *Chi2L = Chi2LL;
    170197                        VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
    171198                        //std::cout << "Vertex fit created " << std::endl;
     
    173200                        //
    174201                        // Remove tracks with large chi2
    175                         Double_t MaxChi2Fit = 9.;
     202                        Double_t MaxChi2Fit = 4.5;
    176203                        Bool_t Done = kFALSE;
    177204                        while (!Done) {
     
    180207                                TVectorD Chi2List = Vtx->GetVtxChi2List();      // Get contributions to Chi2
    181208                                //std::cout << "After Chi2List.  " << std::endl; Chi2List.Print();
    182                                 Double_t* Chi2L = new Double_t[Nfound];
     209                                //Double_t* Chi2L = new Double_t[Nfound];
    183210                                Chi2L = Chi2List.GetMatrixArray();
    184211                                Int_t iMax = TMath::LocMax(Nfound, Chi2L);
     
    186213                                Double_t Chi2Mx = Chi2L[iMax];
    187214                                //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
    188                                 if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
     215                                if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
    189216                                        //std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
    190217                                        Vtx->RemoveTrk(iMax);
     
    195222                                        Done = kTRUE;
    196223                                }
    197                                 //std::cout << "Done =  " << Done << std::endl;
    198                                 //delete[] Chi2L;
    199                                 //std::cout << "Array Chi2L removed " << std::endl;
    200224                        }
    201225                        //
     
    203227                        //
    204228                        // Require minimum number of tracks in vertex
    205                         Int_t Nmin = 2;
     229                        Int_t Nmin = 1;
    206230                        if (Nfound >= Nmin) {
    207231                                TVectorD xvtx = Vtx->GetVtx();
    208                                 //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
     232                                std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
    209233                                TMatrixDSym covX = Vtx->GetVtxCov();
    210234                                Double_t Chi2 = Vtx->GetVtxChi2();
     
    222246                                hTrPrim->Fill(Nprim);
    223247                                hTrFound->Fill((Double_t)Nfound);
     248                                hTrDiff->Fill((Double_t)Nfound-Nprim);
    224249                                //std::cout << "Histograms filled " << std::endl;
    225250                        }
    226251                        //
     252                        // Clean
    227253                        delete Vtx;
     254                        for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
     255                        for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
     256                        delete[] PrFit;
     257                        delete[] CvFit;
    228258                }
    229259
     
    233263                for (Int_t i = 0; i < NtrG; i++) delete pr[i];
    234264                for (Int_t i = 0; i < NtrG; i++) delete cv[i];
    235                 for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
    236                 for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
    237265                delete[] pr;
    238266                delete[] cv;
    239                 delete[] PrFit;
    240                 delete[] CvFit;
    241267        }
    242268        //
     
    254280        //
    255281        TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500);
    256         CnvN->Divide(2, 1);
     282        CnvN->Divide(2, 2);
    257283        CnvN->cd(1);
    258284        hTrPrim->Draw();
     
    260286        hTrFound->SetLineColor(kRed);
    261287        hTrFound->Draw();
     288        CnvN->cd(3);
     289        hTrDiff->Draw();
    262290}
  • examples/ExamplePVtxFitEC.C

    r0c0c9af r78ce8d1  
    11/*
    22Example of using vertex fitter class to fit primary vertex
    3 assumed to be generated in (0,0,0)
     3assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
    44*/
    55
     
    2121void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
    2222{
     23        //
     24        // Beam constraint
     25        // (consistent with generation with ee_zh_smr-shf.cmnd)
     26        //
     27        // Mean beam position
     28        TVectorD xpvc(3);
     29        xpvc(0) = 1.0;
     30        xpvc(1) = -2.0;
     31        xpvc(2) = 10.0;
     32        // Interaction region covariance
     33        TMatrixDSym covpvc(3); covpvc.Zero();
     34        covpvc(0, 0) = 0.0097 * 0.0097;
     35        covpvc(1, 1) = 2.55e-05 * 2.55e-05;
     36        covpvc(2, 2) = 0.64 * 0.64;
     37
     38        //
    2339        // Create chain of root trees
    2440        TChain chain("Delphes");
     
    3551        // Book histograms
    3652        Int_t Nbin = 100;
     53        // Vertex fit pulls
    3754        TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
    3855        TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
    3956        TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
    4057        TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
     58        // Generation check
     59        TH1D* hXvpull = new TH1D("hXvpull", "Pull X generated vertex", Nbin, -10., 10.);
     60        TH1D* hYvpull = new TH1D("hYvpull", "Pull Y generated vertex", Nbin, -10., 10.);
     61        TH1D* hZvpull = new TH1D("hZvpull", "Pull Z generated vertex", Nbin, -10., 10.);
     62
    4163        //
    4264        // Loop over all events
     
    6587                                // Get associated generated particle
    6688                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    67                                 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
     89                                TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);       // mm -> meters
    6890                                TVector3 pg(gp->Px,gp->Py,gp->Pz);
    6991                                Double_t Q = (Double_t)gp->Charge;
    7092                                Double_t Bz = 2.0;
    7193                                TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
    72                                 TVectorD genPar = TrkUtil::ParToMm(genParM);
     94                                TVectorD genPar = TrkUtil::ParToMm(genParM);            // -> back to mm
     95
     96                                //
     97                                // Check if original track is from primary vertex
    7398                                //
    7499                                // Position of origin in mm
     
    79104                                Int_t mp = gp->M1;      // Mother
    80105                                while(mp>0){
    81                                         GenParticle* gm =
    82                                         (GenParticle*)branchGenPart->At(mp);
     106                                        GenParticle* gm = (GenParticle*)branchGenPart->At(mp);
    83107                                        Double_t xm = gm->X;
    84108                                        Double_t ym = gm->Y;
    85109                                        Double_t zm = gm->Z;
    86110                                        if(x!=xm || y!=ym || z!=zm){
    87                                                 prim = kFALSE;
     111                                                prim = kFALSE;  // Not primary
    88112                                                break;
    89113                                        }else mp = gm->M1;
    90114                                }
     115
    91116                                //
    92117                                // group tracks originating from the primary vertex
    93118                                if (prim)
    94119                                {
    95                                         //std::cout<<"Event: "<<entry<<", track: "<<it;
    96                                         //std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl;
     120                                        //
     121                                        // Store true primary vertex for this event
    97122                                        xpv = x;
    98123                                        ypv = y;
     
    104129                                        Double_t obsC = trk->C;
    105130                                        Double_t obsZ0 = trk->DZ;
    106                                         //std::cout<<"Z0 track = "<< obsZ0
    107                                         //<<", gen Z0 = "<<genPar(3)
    108                                         //<<", gen cot = "<<genPar(4)<<std::endl;
    109131                                        Double_t obsCtg = trk->CtgTheta;
    110132                                        Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
     
    112134                                        //
    113135                                        pr[Ntr] = new TVectorD(obsPar);
    114                                         //std::cout<<"Cov Matrix:"<<std::endl;
    115                                         //trk->CovarianceMatrix().Print();
    116136                                        cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
    117137                                        Ntr++;
     
    122142                // Fit primary vertex with beam constraint
    123143                //
    124                 //Beam constraint
    125                 TVectorD xpvc(3);
    126                 xpvc(0) = 1.0;
    127                 xpvc(1) = -2.0;
    128                 xpvc(2) = 10.0;
    129                 TMatrixDSym covpvc(3); covpvc.Zero();
    130                 covpvc(0, 0) = 0.0097 * 0.0097;
    131                 covpvc(1, 1) = 2.55e-05 * 2.55e-05;
    132                 covpvc(2, 2) = 0.64 * 0.64;
    133                 Int_t MinTrk = 2;       // Minumum # tracks for vertex fit
     144                Int_t MinTrk = 3;       // Minumum # tracks for vertex fit
    134145                if (Ntr >= MinTrk) {
    135146                        VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
     
    139150                        Double_t Chi2 = Vtx->GetVtxChi2();
    140151                        Double_t Ndof = 2 * (Double_t)Ntr;
     152                        delete Vtx;
     153                        //
    141154                        Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
    142155                        Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
    143156                        Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
    144                         //std::cout<<"**** True  vertex (x, y, z) "<<xpv<<", "
    145                         //<<ypv<<", "<<zpv<<std::endl;
    146                         //std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
    147                         //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
     157                        //
     158                        Double_t PullXv = (xpvc(0)-xpv) / TMath::Sqrt(covpvc(0, 0));
     159                        Double_t PullYv = (xpvc(1)-ypv) / TMath::Sqrt(covpvc(1, 1));
     160                        Double_t PullZv = (xpvc(2)-zpv) / TMath::Sqrt(covpvc(2, 2));
     161                        //
    148162                        //
    149163                        // Fill histograms
     
    152166                        hZpull->Fill(PullZ);
    153167                        hChi2->Fill(Chi2 / Ndof);
     168                        //
     169                        hXvpull->Fill(PullXv);
     170                        hYvpull->Fill(PullYv);
     171                        hZvpull->Fill(PullZv);
    154172                }
    155173
    156                 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
    157174                //
    158175                // Cleanup
     
    162179                delete[] cv;
    163180        }
     181
    164182        //
    165183        // Show resulting histograms
     
    174192        hZpull->Fit("gaus"); hZpull->Draw();
    175193        Cnv->cd(4); hChi2->Draw();
     194        //
     195        TCanvas* Cnv1 = new TCanvas("Cnv1", "Generated primary vertex pulls", 100, 100, 900, 500);
     196        Cnv1->Divide(3, 1);
     197        Cnv1->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     198        hXvpull->Fit("gaus"); hXvpull->Draw();
     199        Cnv1->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     200        hYvpull->Fit("gaus"); hYvpull->Draw();
     201        Cnv1->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     202        hZvpull->Fit("gaus"); hZvpull->Draw();
    176203}
  • external/TrackCovariance/TrkUtil.h

    r0c0c9af r78ce8d1  
    6868        //
    6969        static TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz);
    70         static TVector3 ParToX(TVectorD Par);                           // position of minimum distance from z axis
     70        static TVector3 ParToX(TVectorD Par);                   // position of minimum distance from z axis
    7171        static TVector3 ParToP(TVectorD Par, Double_t Bz);      // Get Momentum from track parameters
    72         static Double_t ParToQ(TVectorD Par);                           // Get track charge
     72        static Double_t ParToQ(TVectorD Par);                   // Get track charge
    7373        static void LineDistance(TVector3 x0, TVector3 y0, TVector3 dirx, TVector3 diry, Double_t &sx, Double_t &sy, Double_t &distance);
    7474        //
    7575        // Track trajectory
    7676        //
    77         static TVector3 Xtrack(TVectorD par, Double_t s);               // Parametric track trajectory
    78         TVectorD derRphi_R(TVectorD par, Double_t R);   // Derivatives of R-phi at constant R
     77        static TVector3 Xtrack(TVectorD par, Double_t s);       // Parametric track trajectory
     78        TVectorD derRphi_R(TVectorD par, Double_t R);           // Derivatives of R-phi at constant R
    7979        TVectorD derZ_R(TVectorD par, Double_t R);              // Derivatives of z at constant R
    80         TVectorD derRphi_Z(TVectorD par, Double_t z);   // Derivatives of R-phi at constant z
     80        TVectorD derRphi_Z(TVectorD par, Double_t z);           // Derivatives of R-phi at constant z
    8181        TVectorD derR_Z(TVectorD par, Double_t z);              // Derivatives of R at constant z
    8282        //
     
    8888        //
    8989        static TVectorD ParToMm(TVectorD Par);                  // Parameter conversion
    90         static TMatrixDSym CovToMm(TMatrixDSym Cov);    // Covariance conversion
     90        static TMatrixDSym CovToMm(TMatrixDSym Cov);            // Covariance conversion
    9191        //
    9292        // Inside cylindrical volume
  • external/TrackCovariance/VertexFit.cc

    r0c0c9af r78ce8d1  
    1 #include <TMath.h>
     1/*
     2Vertex fitting code
     3*/
     4#include <TMath.h>
    25#include <TVectorD.h>
    36#include <TVector3.h>
    47#include <TMatrixD.h>
    58#include <TMatrixDSym.h>
     9#include <TMatrixDSymEigen.h>
    610#include "VertexFit.h"
    711//
     
    3842        for (Int_t i = 0; i < fNtr; i++)
    3943        {
    40                 TVectorD pr = *trkPar[i];
    41                 fPar.push_back(new TVectorD(pr));
    42                 fParNew.push_back(new TVectorD(pr));
    43                 TMatrixDSym cv = *trkCov[i];
    44                 fCov.push_back(new TMatrixDSym(cv));
    45                 fCovNew.push_back(new TMatrixDSym(cv));
     44                fPar.push_back(new TVectorD(*trkPar[i]));
     45                fParNew.push_back(new TVectorD(*trkPar[i]));
     46                fCov.push_back(new TMatrixDSym(*trkCov[i]));
     47                fCovNew.push_back(new TMatrixDSym(*trkCov[i]));
    4648        }
    4749        fChi2List.ResizeTo(fNtr);
     
    7678void VertexFit::ResetWrkArrays()
    7779{
    78         Int_t N = (Int_t)ffi.size();
    79         for (Int_t i = 0; i < N; i++)
    80         {
    81                 if (fx0i[i])  { fx0i[i]->Clear();               delete fx0i[i]; }
    82                 if (fai[i])   { fai[i]->Clear();                delete fai[i]; }
    83                 if (fdi[i])   { fdi[i]->Clear();                delete fdi[i]; }
    84                 if (fAti[i])  { fAti[i]->Clear();               delete fAti[i]; }
    85                 if (fDi[i])   { fDi[i]->Clear();                delete fDi[i]; }
    86                 if (fWi[i])   { fWi[i]->Clear();                delete fWi[i]; }
    87                 if (fWinvi[i]){ fWinvi[i]->Clear();     delete fWinvi[i]; }
    88         }
    89         fa2i.clear();
    90         fx0i.clear();
    91         fai.clear();
    92         fdi.clear();
    93         fAti.clear();
    94         fDi.clear();
    95         fWi.clear();
    96         fWinvi.clear();
     80        Int_t N = (Int_t)fdi.size();
     81        if(N > 0){
     82                for (Int_t i = 0; i < N; i++)
     83                {
     84                        if (fx0i[i])  { fx0i[i]->Clear();       delete fx0i[i]; }
     85                        if (fai[i])   { fai[i]->Clear();        delete fai[i]; }
     86                        if (fdi[i])   { fdi[i]->Clear();        delete fdi[i]; }
     87                        if (fAti[i])  { fAti[i]->Clear();       delete fAti[i]; }
     88                        if (fDi[i])   { fDi[i]->Clear();        delete fDi[i]; }
     89                        if (fWi[i])   { fWi[i]->Clear();        delete fWi[i]; }
     90                        if (fWinvi[i]){ fWinvi[i]->Clear();     delete fWinvi[i]; }
     91                }
     92                fa2i.clear();
     93                fx0i.clear();
     94                fai.clear();
     95                fdi.clear();
     96                fAti.clear();
     97                fDi.clear();
     98                fWi.clear();
     99                fWinvi.clear();
     100        }
    97101}
    98102VertexFit::~VertexFit()
    99 {
     103{       
    100104        fxCst.Clear();         
    101105        fCovCst.Clear();
     
    103107        fXv.Clear();           
    104108        fcovXv.Clear();         
    105         fChi2List.Clear();     
     109        fChi2List.Clear();
    106110        //
    107111        for (Int_t i = 0; i < fNtr; i++)
     
    111115                fCov[i]->Clear();               delete fCov[i];
    112116                fCovNew[i]->Clear();    delete fCovNew[i];
    113         }
     117        }       
    114118        fPar.clear();
    115119        fParNew.clear();
     
    118122        //
    119123        ResetWrkArrays();
    120         ffi.clear();
     124        ffi.clear();   
    121125        fNtr = 0;
    122126}
    123127//
    124 Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2)
    125 {
    126         //
    127         // Find radius of minimum distance between two tracks
    128         //
    129         // p = (D,phi, C, z0, ct)
    130         //
    131         // Define arrays
    132         //
    133         Double_t C1 = p1(2);
    134         Double_t C2 = p2(2);
    135         Double_t ph1 = p1(1);
    136         Double_t ph2 = p2(1);
    137         TVectorD x0 = Fill_x0(p1);
    138         TVectorD y0 = Fill_x0(p2);
    139         TVectorD n = Fill_a(p1, 0.0);
    140         n *= (2*C1);
    141         TVectorD k = Fill_a(p2, 0.0);
    142         k *= (2*C2);
    143         //
    144         // Setup and solve linear system
    145         //
    146         Double_t nn = 0; for (Int_t i = 0; i < 3; i++)nn += n(i) * n(i);
    147         Double_t nk = 0; for (Int_t i = 0; i < 3; i++)nk += n(i) * k(i);
    148         Double_t kk = 0; for (Int_t i = 0; i < 3; i++)kk += k(i) * k(i);
    149         Double_t discr = nn * kk - nk * nk;
    150         TMatrixDSym H(2);
    151         H(0, 0) = kk;
    152         H(0, 1) = nk;
    153         H(1, 0) = H(0, 1);
    154         H(1, 1) = nn;
    155         TVectorD c(2);
    156         c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i));
    157         c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i));
    158         TVectorD smin = (H * c);
    159         smin *= 1.0 / discr;
    160         //
    161         TVectorD X = x0 + smin(0) * n;
    162         TVectorD Y = y0 + smin(1) * k;
    163         //
    164         // Higher order corrections
    165         X(0) += -C1 * smin(0) * smin(0) * TMath::Sin(ph1);
    166         X(1) +=  C1 * smin(0) * smin(0) * TMath::Cos(ph1);
    167         Y(0) += -C2 * smin(1) * smin(1) * TMath::Sin(ph2);
    168         Y(1) +=  C2 * smin(1) * smin(1) * TMath::Cos(ph2);
    169         //
    170         TVectorD Xavg = 0.5 * (X + Y);
    171         //
    172         //
    173         return TMath::Sqrt(Xavg(0)*Xavg(0)+Xavg(1)*Xavg(1));
    174 }
    175 //
    176 // Starting radius determination
    177 Double_t VertexFit::StartRadius()
    178 {
    179         //
    180         // Maximum impact parameter
    181         Double_t Rd = 0;
    182         for (Int_t i = 0; i < fNtr; i++)
    183         {
    184                 TVectorD par = *fPar[i];
    185                 Double_t Dabs = TMath::Abs(par(0));
    186                 if (Dabs > Rd)Rd = Dabs;
    187         }
    188         //-----------------------------
    189         //
    190         // Find track pair with phi difference closest to pi/2
    191         Int_t isel = 0; Int_t jsel = 0;         // selected track indices
    192         Double_t dSinMax = 0.0;                         // Max phi difference
    193         for (Int_t i = 0; i < fNtr - 1; i++)
    194         {
    195                 TVectorD pari = *fPar[i];
    196                 Double_t phi1 = pari(1);
    197 
    198                 for (Int_t j = i + 1; j < fNtr; j++)
    199                 {
    200                         TVectorD parj = *fPar[j];
    201                         Double_t phi2 = parj(1);
    202                         Double_t Sindphi = TMath::Abs(TMath::Sin(phi2 - phi1));
    203                         if (Sindphi > dSinMax)
    204                         {
    205                                 isel = i; jsel = j;
    206                                 dSinMax = Sindphi;
    207                         }
    208                 }
    209         }
    210         //
    211         //------------------------------------------
    212         //
    213         // Find radius of minimum distrance between tracks
    214         TVectorD p1 = *fPar[isel];
    215         TVectorD p2 = *fPar[jsel];
    216         Double_t R = FastRv(p1, p2);
    217         //
    218         R = 0.9 * R + 0.1 * Rd;         // Protect for overshoot
    219         //
    220         return R;
    221 }
    222 //
    223 // Regularized symmetric matrix inversion
    224 //
    225 TMatrixDSym VertexFit::RegInv(TMatrixDSym& Min)
    226 {
    227         TMatrixDSym M = Min;                            // Decouple from input
    228         Int_t N = M.GetNrows();                 // Matrix size
    229         TMatrixDSym D(N); D.Zero();             // Normaliztion matrix
    230         TMatrixDSym R(N);                               // Normarized matrix
    231         TMatrixDSym Rinv(N);                            // Inverse of R
    232         TMatrixDSym Minv(N);                            // Inverse of M
    233         //
    234         // Check for 0's and normalize
    235         for (Int_t i = 0; i < N; i++)
    236         {
    237                 if (M(i, i) != 0.0) D(i, i) = 1. / TMath::Sqrt(TMath::Abs(M(i, i)));
    238                 else D(i, i) = 1.0;
    239         }
    240         R = M.Similarity(D);
    241         //
    242         // Recursive algorithms stops when N = 2
    243         //
    244         //****************
    245         // case N = 2  ***
    246         //****************
    247         if (N == 2)
    248         {
    249                 Double_t det = R(0, 0) * R(1, 1) - R(0, 1) * R(1, 0);
    250                 if (det == 0)
    251                 {
    252                         std::cout << "VertexFit::RegInv: null determinant for N = 2" << std::endl;
    253                         Rinv.Zero();    // Return null matrix
    254                 }
    255                 else
    256                 {
    257                         // invert matrix
    258                         Rinv(0, 0) = R(1, 1);
    259                         Rinv(0, 1) = -R(0, 1);
    260                         Rinv(1, 0) = Rinv(0, 1);
    261                         Rinv(1, 1) = R(0, 0);
    262                         Rinv *= 1. / det;
    263                 }
    264         }
    265         //****************
    266         // case N > 2  ***
    267         //****************
    268         else
    269         {
    270                 // Break up matrix
    271                 TMatrixDSym Q = R.GetSub(0, N - 2, 0, N - 2);   // Upper left
    272                 TVectorD p(N - 1);
    273                 for (Int_t i = 0; i < N - 1; i++)p(i) = R(N - 1, i);
    274                 Double_t q = R(N - 1, N - 1);
    275                 //Invert pieces and re-assemble
    276                 TMatrixDSym Ainv(N - 1);
    277                 TMatrixDSym A(N - 1);
    278                 if (TMath::Abs(q) > 1.0e-15)
    279                 {
    280                         // Case |q| > 0
    281                         Ainv.Rank1Update(p, -1.0 / q);
    282                         Ainv += Q;
    283                         A = RegInv(Ainv);               // Recursive call
    284                         TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A;
    285                         //
    286                         TVectorD b = (-1.0 / q) * (A * p);
    287                         for (Int_t i = 0; i < N - 1; i++)
    288                         {
    289                                 Rinv(N - 1, i) = b(i);
    290                                 Rinv(i, N - 1) = b(i);
    291                         }
    292                         //
    293                         Double_t pdotb = 0.;
    294                         for (Int_t i = 0; i < N - 1; i++)pdotb += p(i) * b(i);
    295                         Double_t c = (1.0 - pdotb) / q;
    296                         Rinv(N - 1, N - 1) = c;
    297                 }
    298                 else
    299                 {
    300                         // case q = 0
    301                         TMatrixDSym Qinv = RegInv(Q);           // Recursive call
    302                         Double_t a = Qinv.Similarity(p);
    303                         Double_t c = -1.0 / a;
    304                         Rinv(N - 1, N - 1) = c;
    305                         //
    306                         TVectorD b = (1.0 / a) * (Qinv * p);
    307                         for (Int_t i = 0; i < N - 1; i++)
    308                         {
    309                                 Rinv(N - 1, i) = b(i);
    310                                 Rinv(i, N - 1) = b(i);
    311                         }
    312                         //
    313                         A.Rank1Update(p, -1 / a);
    314                         A += Q;
    315                         A.Similarity(Qinv);
    316                         TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A;
    317                 }
    318         }
    319         Minv = Rinv.Similarity(D);
    320         return Minv;
    321 }
    322 //
    323 //
    324 //
    325 TMatrixD VertexFit::Fill_A(TVectorD par, Double_t phi)
    326 {
    327         //
    328         // Derivative of track 3D position vector with respect to track parameters at constant phase
    329         //
    330         // par = vector of track parameters
    331         // phi = phase
    332         //
    333         TMatrixD A(3, 5);
     128//
     129//
     130TVectorD VertexFit::Fill_x0(TVectorD par)
     131{
     132        //
     133        // Calculate track 3D position at R = |D| (minimum approach to z-axis)
     134        //
     135        TVectorD x0(3);
    334136        //
    335137        // Decode input arrays
     
    341143        Double_t ct = par(4);
    342144        //
    343         // Fill derivative matrix dx/d alpha
    344         // D
    345         A(0, 0) = -TMath::Sin(p0);
    346         A(1, 0) = TMath::Cos(p0);
    347         A(2, 0) = 0.0;
    348         // phi0
    349         A(0, 1) = -D * TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
    350         A(1, 1) = -D * TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
    351         A(2, 1) = 0.0;
    352         // C
    353         A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C * C);
    354         A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C * C);
    355         A(2, 2) = -ct * phi / (2 * C * C);
    356         // z0
    357         A(0, 3) = 0.0;
    358         A(1, 3) = 0.0;
    359         A(2, 3) = 1.0;
    360         // ct = lambda
    361         A(0, 4) = 0.0;
    362         A(1, 4) = 0.0;
    363         A(2, 4) = phi / (2 * C);
    364         //
    365         return A;
    366 }
    367 //
    368 TVectorD VertexFit::Fill_a(TVectorD par, Double_t phi)
    369 {
    370         //
    371         // Derivative of track 3D position vector with respect to phase at constant track parameters
    372         //
    373         // par = vector of track parameters
    374         // phi = phase
    375         //
    376         TVectorD a(3);
     145        x0(0) = -D * TMath::Sin(p0);
     146        x0(1) = D * TMath::Cos(p0);
     147        x0(2) = z0;
     148        //
     149        return x0;
     150}
     151//
     152TVectorD VertexFit::Fill_x(TVectorD par, Double_t phi)
     153{
     154        //
     155        // Calculate track 3D position for a given phase, phi
     156        //
     157        TVectorD x(3);
    377158        //
    378159        // Decode input arrays
     
    384165        Double_t ct = par(4);
    385166        //
    386         a(0) = TMath::Cos(phi + p0) / (2 * C);
    387         a(1) = TMath::Sin(phi + p0) / (2 * C);
    388         a(2) = ct / (2 * C);
    389         //
    390         return a;
    391 }
    392 //
    393 TVectorD VertexFit::Fill_x0(TVectorD par)
    394 {
    395         //
    396         // Calculate track 3D position at R = |D| (minimum approach to z-axis)
    397         //
    398         TVectorD x0(3);
    399         //
    400         // Decode input arrays
    401         //
    402         Double_t D = par(0);
    403         Double_t p0 = par(1);
    404         Double_t C = par(2);
    405         Double_t z0 = par(3);
    406         Double_t ct = par(4);
    407         //
    408         x0(0) = -D * TMath::Sin(p0);
    409         x0(1) = D * TMath::Cos(p0);
    410         x0(2) = z0;
    411         //
    412         return x0;
    413 }
    414 //
    415 TVectorD VertexFit::Fill_x(TVectorD par, Double_t phi)
    416 {
    417         //
    418         // Calculate track 3D position for a given phase, phi
    419         //
    420         TVectorD x(3);
    421         //
    422         // Decode input arrays
    423         //
    424         Double_t D = par(0);
    425         Double_t p0 = par(1);
    426         Double_t C = par(2);
    427         Double_t z0 = par(3);
    428         Double_t ct = par(4);
    429         //
    430167        TVectorD x0 = Fill_x0(par);
    431168        x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
     
    439176{
    440177        //
    441         // Get track parameters and their covariance
     178        // Get track parameters, covariance and phase
     179        Double_t fs = ffi[i];                   // Get phase
    442180        TVectorD par = *fParNew[i];
    443181        TMatrixDSym Cov = *fCov[i];
    444182        //
    445183        // Fill all track related work arrays arrays
    446         Double_t fs = ffi[i];                                           // Get phase
     184        TMatrixD A = derXdPar(par, fs);         // A = dx/da = derivatives wrt track parameters
     185        TMatrixDSym Winv = Cov;                         
     186        Winv.Similarity(A);                     // W^-1 = A*C*A'
     187
     188        TMatrixD At(TMatrixD::kTransposed, A);          // A transposed
     189        fAti.push_back(new TMatrixD(At));               // Store A'
     190        fWinvi.push_back(new TMatrixDSym(Winv));        // Store W^-1 matrix
     191        //
    447192        TVectorD xs = Fill_x(par, fs);
    448193        fx0i.push_back(new TVectorD(xs));                       // Start helix position
    449         //
    450         TMatrixD A = Fill_A(par, fs);                           // A = dx/da = derivatives wrt track parameters
    451         TMatrixD At(TMatrixD::kTransposed, A);          // A transposed
    452         fAti.push_back(new TMatrixD(At));                       // Store A'
    453         TVectorD di = A * (par - *fPar[i]);             // x-shift
    454         fdi.push_back(new TVectorD(di));                        // Store x-shift
    455         TMatrixDSym Winv = Cov.Similarity(A);           // W^-1 = A*C*A'
    456         fWinvi.push_back(new TMatrixDSym(Winv));        // Store W^-1 matrix
     194        //
     195        TVectorD di = A * (par - *fPar[i]);                     // x-shift
     196        fdi.push_back(new TVectorD(di));                        // Store x-shift       
    457197        //
    458198        TMatrixDSym W = RegInv(Winv);                           // W = (A*C*A')^-1
    459199        fWi.push_back(new TMatrixDSym(W));                      // Store W matrix
    460200        //
    461         TVectorD a = Fill_a(par, fs);                           // a = dx/ds = derivatives wrt phase
     201        TVectorD a = derXds(par, fs);                           // a = dx/ds = derivatives wrt phase
    462202        fai.push_back(new TVectorD(a));                         // Store a
    463203        //
     
    473213}
    474214//
     215void VertexFit::VtxFitNoSteer()
     216{
     217        //
     218        // Initialize
     219        //
     220        std::vector<TVectorD*> x0i;                             // Tracks at ma
     221        std::vector<TVectorD*> ni;                              // Track derivative wrt phase
     222        std::vector<TMatrixDSym*> Ci;                   // Position error matrix at fixed phase
     223        std::vector<TVectorD*> wi;                              // Ci*ni
     224        //
     225        // Track loop
     226        for (Int_t i = 0; i < fNtr; i++)
     227        {
     228                Double_t s = 0.;
     229                TVectorD par = *fPar[i];
     230                TMatrixDSym Cov = *fCov[i];
     231                x0i.push_back(new TVectorD(Fill_x0(par)));
     232                ni.push_back(new TVectorD(derXds(par, s)));
     233                TMatrixD A = derXdPar(par, s);
     234                Ci.push_back(new TMatrixDSym(Cov.Similarity(A)));
     235                TMatrixDSym Cinv = RegInv(*Ci[i]);
     236                wi.push_back(new TVectorD(Cinv * (*ni[i])));
     237        }
     238        //std::cout << "Vtx init completed. fNtr = "<<fNtr << std::endl;
     239        //
     240        // Get fit vertex
     241        //
     242        TMatrixDSym D(3); D.Zero();
     243        TVectorD Dx(3); Dx.Zero();
     244        for (Int_t i = 0; i < fNtr; i++)
     245        {
     246                TMatrixDSym Cinv = RegInv(*Ci[i]);
     247                TMatrixDSym W(3);
     248                W.Rank1Update(*wi[i], 1. / Ci[i]->Similarity(*wi[i]));
     249                TMatrixDSym Dd = Cinv - W;
     250                D += Dd;
     251                Dx += Dd * (*x0i[i]);
     252        }
     253        if(fVtxCst){
     254                D  += fCovCstInv;
     255                Dx += fCovCstInv*fxCst;
     256        }
     257        fXv = RegInv(D) * Dx;
     258        //std::cout << "Fast vertex (x, y, z) = "<<fXv(0)<<", "<<fXv(1)<<", "<<fXv(2) << std::endl;
     259        //
     260        // Get fit phases
     261        //
     262        for (Int_t i = 0; i < fNtr; i++){
     263                Double_t si = Dot(*wi[i], fXv - (*x0i[i])) / Ci[i]->Similarity(*wi[i]);
     264                ffi.push_back(si);
     265                //TVectorD xvi = Fill_x(*fPar[i],si);
     266                //std::cout << "Fast vertex "<<i<<": xvi = "<<xvi(0)<<", "<<xvi(1)<<", "<<xvi(2)
     267                //      <<", si = "<<si<< std::endl;
     268        }
     269        //
     270        // Cleanup
     271        for (Int_t i = 0; i < fNtr; i++)
     272        {
     273                x0i[i]->Clear();        delete x0i[i];
     274                ni[i]->Clear();         delete ni[i];
     275                Ci[i]->Clear();         delete Ci[i];
     276                wi[i]->Clear();         delete wi[i];
     277        }
     278        x0i.clear();;
     279        ni.clear();
     280        Ci.clear();
     281        wi.clear();
     282}
     283//
    475284void  VertexFit::VertexFitter()
    476285{
     
    487296        TMatrixDSym covX(3);
    488297        Double_t Chi2 = 0;
    489         TVectorD x0 = fXv;      // If previous fit done
    490         if (fRold < 0.0) {
    491                 // External constraint
    492                 if (fVtxCst) fRold = TMath::Sqrt(fxCst(0) * fxCst(0) + fxCst(1) * fxCst(1));
    493                 // No constraint
    494                 else for (Int_t i = 0; i < 3; i++)x0(i) = 1000.;        // Set to arbitrary large value
    495         }
    496         //
    497         // Starting vertex radius approximation
    498         //
    499         Double_t R = fRold;                                             // Use previous fit if available
    500         if (R < 0.0) R = StartRadius();                 // Rough vertex estimate
    501         //std::cout << "Start radius: " << R << std::endl;
     298        //
     299        VtxFitNoSteer();        // Fast vertex finder on first pass (set ffi and fXv)
     300        TVectorD x0 = fXv;
    502301        //
    503302        // Iteration properties
     
    515314                TVectorD cterm(3); TMatrixDSym H(3); TMatrixDSym DW1D(3);
    516315                covX.Zero();            // Reset vertex covariance
    517                 cterm.Zero();   // Reset constant term
     316                cterm.Zero();           // Reset constant term
    518317                H.Zero();               // Reset H matrix
    519318                DW1D.Zero();
     
    530329                        TVectorD par = *fPar[i];
    531330                        TMatrixDSym Cov = *fCov[i];
    532                         //
    533                         // For first iteration only
    534                         Double_t fs;
    535                         if (Ntry <= 0)  // Initialize all phases on first pass
    536                         {
    537                                 Double_t D = par(0);
    538                                 Double_t C = par(2);
    539                                 Double_t arg = TMath::Max(1.0e-6, (R * R - D * D) / (1 + 2 * C * D));
    540                                 fs = 2 * TMath::ASin(C * TMath::Sqrt(arg));
    541                                 ffi.push_back(fs);
    542                         }
    543331                        //
    544332                        // Update track related arrays
     
    555343                        // update constant term
    556344                        TVectorD xs = *fx0i[i] - *fdi[i];
    557                         TVectorD xx0 = *fx0i[i];
    558                         /*
    559                         std::cout << "Iter. " << Ntry << ", trk " << i << ", x= "
    560                                 << xx0(0) << ", " << xx0(1) << ", " << xx0(2)<<
    561                                 ", ph0= "<<par(1)<< std::endl;
    562                         */
     345                        //TVectorD xx0 = *fx0i[i];
     346                       
     347                        //std::cout << "Iter. " << Ntry << ", trk " << i << ", xs= "
     348                        //      << xs(0) << ", " << xs(1) << ", " << xs(2)<<
     349                        //      ", ph0= "<<par(1)<< std::endl;
     350                       
    563351                        cterm += Ds * xs;
    564352                }                               // End loop on tracks
     
    584372                        TMatrixDSym Wm1 = *fWinvi[i];
    585373                        fChi2List(i) = Wm1.Similarity(lambda);
     374                        if(fChi2List(i) < 0.0){
     375                                //std::cout<<"# "<<i<<", Chi2= "<<fChi2List(i)<<", Wm1:"<<std::endl; Wm1.Print();
     376                                //std::cout<<"Lambda= "<<std::endl; lambda.Print();
     377                        }
    586378                        Chi2 += fChi2List(i);
    587379                        TVectorD a = *fai[i];
    588                         TVectorD b = (*fWi[i]) * (x - (*fx0i[i]));
    589                         for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i];
     380                        TVectorD b = (*fWi[i]) * (x - *fx0i[i] + *fdi[i]);
     381                        ffi[i] += Dot(a, b) / fa2i[i];
    590382                        TVectorD newPar = *fPar[i] - ((*fCov[i]) * (*fAti[i])) * lambda;
    591383                        fParNew[i] = new TVectorD(newPar);
     
    610402        //
    611403        fVtxDone = kTRUE;               // Set fit completion flag
    612         fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1));     // Store fit radius
     404        fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1));     // Store fit
     405        //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2)
     406        //      << ", after "<<Ntry<<" iterations"<<std::endl;
    613407        //
    614408}
  • external/TrackCovariance/VertexFit.h

    r0c0c9af r78ce8d1  
    66#include <TVectorD.h>
    77#include <TMatrixDSym.h>
     8#include "TrkUtil.h"
    89#include "ObsTrk.h"
    910#include <vector>
     
    1213// Class for vertex fitting
    1314
    14 class VertexFit {
     15class VertexFit: public TrkUtil
     16{
    1517        //
    1618        // Vertex fitting with track parameters steering
     
    2123        //
    2224        // Inputs
    23         Int_t fNtr;                                     // Number of tracks
     25        Int_t fNtr;                                                     // Number of tracks
    2426        std::vector<TVectorD*> fPar;            // Input parameter array
    2527        std::vector<TVectorD*> fParNew;         // Updated parameter array
     
    2729        std::vector<TMatrixDSym*> fCovNew;      // Updated parameter covariances
    2830        // Constraints
    29         Bool_t fVtxCst;                         // Vertex constraint flag
    30         TVectorD fxCst;                         // Constraint value
     31        Bool_t fVtxCst;                                 // Vertex constraint flag
     32        TVectorD fxCst;                                 // Constraint value
    3133        TMatrixDSym fCovCst;                    // Constraint
    3234        TMatrixDSym fCovCstInv;                 // Inverse of constraint covariance
    3335        //
    3436        // Results
    35         Bool_t fVtxDone;                                // Flag vertex fit completed
     37        Bool_t fVtxDone;                        // Flag vertex fit completed
    3638        Double_t fRold;                         // Current value of vertex radius
    3739        TVectorD fXv;                           // Found vertex
     
    5254        //
    5355        // Service routines
    54         //void InitWrkArrays();                                                 // Initializations
    5556        void ResetWrkArrays();                                                  // Clear work arrays
    56         Double_t StartRadius();                                                 // Starting vertex radius determination
    57         Double_t FastRv(TVectorD p1, TVectorD p2);              // Fast vertex radius determination
    58         TMatrixDSym RegInv(TMatrixDSym& Smat0);                 // Regularized 3D matrix inversion
    59         TMatrixD Fill_A(TVectorD par, Double_t phi);            // Derivative of track position wrt track parameters
    60         TVectorD Fill_a(TVectorD par, Double_t phi);            // Derivative of track position wrt track phase
     57        //Double_t StartRadius();                                                       // Starting vertex radius determination
     58        //Double_t FastRv(TVectorD p1, TVectorD p2);            // Fast vertex radius determination
     59        //TMatrixDSym RegInv(TMatrixDSym& Smat0);               // Regularized 3D matrix inversion
     60        //TMatrixD Fill_A(TVectorD par, Double_t phi);  // Derivative of track position wrt track parameters
     61        //TVectorD Fill_a(TVectorD par, Double_t phi);  // Derivative of track position wrt track phase
    6162        TVectorD Fill_x0(TVectorD par);                                 // Track position at dma to z-axis
    62         TVectorD Fill_x(TVectorD par, Double_t phi);            // Track position at given phase
     63        TVectorD Fill_x(TVectorD par, Double_t phi);    // Track position at given phase
    6364        void UpdateTrkArrays(Int_t i);                                  // Fill track realted arrays
    64         void VertexFitter();                                                            // Vertex finder routine
     65        void VtxFitNoSteer();                                                   // Vertex fitter routine w/o parameter steering
     66        void VertexFitter();                                                    // Vertex fitter routine w/  parameter steering
    6567public:
    6668        //
  • modules/TreeWriter.cc

    r0c0c9af r78ce8d1  
    369369
    370370    // add some offdiagonal covariance matrix elements
    371     entry->ErrorD0Phi          = candidate->TrackCovariance(0,1);
     371    entry->ErrorD0Phi          = candidate->TrackCovariance(0,1)*1.e3;
    372372    entry->ErrorD0C            = candidate->TrackCovariance(0,2);
    373     entry->ErrorD0DZ           = candidate->TrackCovariance(0,3);
    374     entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4);
    375     entry->ErrorPhiC           = candidate->TrackCovariance(1,2);
    376     entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3);
     373    entry->ErrorD0DZ           = candidate->TrackCovariance(0,3)*1.e6;
     374    entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4)*1.e3;
     375    entry->ErrorPhiC           = candidate->TrackCovariance(1,2)*1.e-3;
     376    entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3)*1.e3;
    377377    entry->ErrorPhiCtgTheta    = candidate->TrackCovariance(1,4);
    378378    entry->ErrorCDZ            = candidate->TrackCovariance(2,3);
    379     entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4);
    380     entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4);
     379    entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4)*1.e-3;
     380    entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4)*1.e3;
    381381
    382382    entry->Xd = candidate->Xd;
Note: See TracChangeset for help on using the changeset viewer.