Fork me on GitHub

Changeset 7dac4ea in git


Ignore:
Timestamp:
May 17, 2021, 6:05:45 PM (3 years ago)
Author:
GitHub <noreply@…>
Branches:
master
Children:
1c1c9c2, 33a6b3a, f8e61b2
Parents:
4afb18d (diff), 4acf2fd (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@…> (05/17/21 18:05:45)
git-committer:
GitHub <noreply@…> (05/17/21 18:05:45)
Message:

Merge pull request #95 from fbedesch/master

Fix for backward tracks and vertex fitting improvements

Files:
2 added
6 edited

Legend:

Unmodified
Added
Removed

  • examples/ExamplePVtxFind.C

    r4afb18d r7dac4ea  
    5151                // Load selected branches with data from specified event
    5252                treeReader->ReadEntry(entry);
    53                 Int_t Ntr = 0;  // # of starting tracks used in primary vertex
    5453                Int_t NtrG = branchTrack->GetEntries();
    55                 //std::cout << "Event opened containing " << NtrG << " tracks" << std::endl;
     54                if(entry%500 ==0)std::cout << "Event "<<entry<<" opened containing " << NtrG << " tracks" << std::endl;
    5655                TVectorD** pr = new TVectorD * [NtrG];
    5756                TMatrixDSym** cv = new TMatrixDSym * [NtrG];
    … …  
    6362                //
    6463                Double_t Nprim = 0.0;
     64                Double_t xpv = 0.0;             // Init true primary vertex position
     65                Double_t ypv = 0.0;
     66                Double_t zpv = 0.0;
    6567                if (branchTrack->GetEntries() > 0)
    6668                {
    … …  
    8082                                //std::cout << "Got track parameters for track " << it << std::endl;
    8183                                //
    82                                 // Load tracks for vertex fit if impact parameters is not ridiculous
    83                                 Double_t Dmax = 1.0;    // max is 1 mm
    84                                 if (TMath::Abs(obsD0) < Dmax) {
    85                                         Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
    86                                         TVectorD obsPar(5, oPar);       // Fill observed parameters
    87                                         pr[Ntr] = new TVectorD(obsPar);
    88                                         cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
    89                                         Ntr++;
    90                                         //std::cout << "Track loaded Ntr= " << Ntr << std::endl;
    91                                 }
    92                                 //
    93                                 // Get associated generated particle
     84                                // Load all tracks for vertex fit
     85                                Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
     86                                TVectorD obsPar(5, oPar);       // Fill observed parameters
     87                                pr[it] = new TVectorD(obsPar);
     88                                cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
     89                                //std::cout << "Track loaded it= " << it << std::endl;
     90                        //
     91                        // Find true primary vertex
    9492                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    9593                                //std::cout << "GenParticle pointer "<<gp << std::endl;
    … …  
    9997                                Double_t y = gp->Y;
    10098                                Double_t z = gp->Z;
    101                                 //std::cout << "Got position of origin " << std::endl;
    102                                 //
    103                                 // Count tracks originating from the primary vertex
    104                                 if (x == 0.0 && y == 0.0) Nprim++;
     99                                Bool_t prim = kTRUE;    // Is primary?
     100                                Int_t mp = gp->M1;      // Mother
     101                                while (mp > 0) {
     102                                        GenParticle* gm =
     103                                                (GenParticle*)branchGenPart->At(mp);
     104                                        Double_t xm = gm->X;
     105                                        Double_t ym = gm->Y;
     106                                        Double_t zm = gm->Z;
     107                                        if (x != xm || y != ym || z != zm) {
     108                                                prim = kFALSE;
     109                                                break;
     110                                        }
     111                                        else mp = gm->M1;
     112                                }
     113                                if (prim) {             // It's a primary track
     114                                        Nprim++;
     115                                        xpv = x;        // Store true primary
     116                                        ypv = y;
     117                                        zpv = z;
     118                                }
    105119
    106120                        }               // End loop on tracks
    107121                }
    108122                //
    109                 // Fit primary vertex
    110                 //
    111                 Int_t Nfound = Ntr;
    112                 //std::cout << "Found tracks "<<Nfound << std::endl;
    113                 Int_t MinTrk = 2;       // Minumum # tracks for vertex fit
    114                 Double_t MaxChi2 = 6.;
    115                 if (Ntr >= MinTrk) {
    116                         VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
     123                // Find primary vertex
     124                //
     125                //Beam constraint
     126                TVectorD xpvc(3);
     127                xpvc(0) = 1.0;
     128                xpvc(1) = -2.0;
     129                xpvc(2) = 10.0;
     130                TMatrixDSym covpvc(3); covpvc.Zero();
     131                covpvc(0, 0) = 0.0097 * 0.0097;
     132                covpvc(1, 1) = 2.55e-05 * 2.55e-05;
     133                covpvc(2, 2) = 0.64 * 0.64;
     134                //
     135                //
     136                // Skim tracks
     137                Int_t nSkim = 0;
     138                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.);
     142                for (Int_t n = 0; n < NtrG; n++) {
     143                        Vskim->AddTrk(pr[n], cv[n]);
     144                        Double_t Chi2One = Vskim->GetVtxChi2();
     145                        //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
     146                        if (Chi2One < MaxChi2) {
     147                                nSkimmed[nSkim] = n;
     148                                //std::cout << "nSkimmed[" << nSkim << "] = " << n << std::endl;
     149                                nSkim++;
     150                        }
     151                        Vskim->RemoveTrk(0);
     152                }
     153                // Load tracks for primary fit
     154                Int_t MinTrk = 1;       // Minumum # tracks for vertex fit
     155                TVectorD** PrFit = new TVectorD * [nSkim];
     156                TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
     157                if (nSkim >= MinTrk) {
     158                        for (Int_t n = 0; n < nSkim; n++) {
     159                                PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
     160                                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) {
     170                        VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
    117171                        //std::cout << "Vertex fit created " << std::endl;
     172                        Vtx->AddVtxConstraint(xpvc, covpvc);
    118173                        //
    119174                        // Remove tracks with large chi2
     175                        Double_t MaxChi2Fit = 9.;
    120176                        Bool_t Done = kFALSE;
    121                         while (!Done){
     177                        while (!Done) {
    122178                                //std::cout << "After while " << std::endl;
    123179                                // Find largest Chi2 contribution
    … …  
    130186                                Double_t Chi2Mx = Chi2L[iMax];
    131187                                //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
    132                                 if (Chi2Mx > MaxChi2 && Nfound > 2) {
    133                                         Vtx->RemoveTrk(iMax);
     188                                if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
     189                                        //std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
     190                                        Vtx->RemoveTrk(iMax);
     191                                        //std::cout << "After remove." << std::endl;
    134192                                        Nfound--;
    135193                                }
    … …  
    145203                        //
    146204                        // Require minimum number of tracks in vertex
    147                         Int_t Nmin = 4;
     205                        Int_t Nmin = 2;
    148206                        if (Nfound >= Nmin) {
    149207                                TVectorD xvtx = Vtx->GetVtx();
    … …  
    151209                                TMatrixDSym covX = Vtx->GetVtxCov();
    152210                                Double_t Chi2 = Vtx->GetVtxChi2();
    153                                 Double_t Ndof = 2 * (Double_t)Nfound - 3;
    154                                 Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
    155                                 Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
    156                                 Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
     211                                Double_t Ndof = 2 * (Double_t)Nfound;
     212                                Double_t PullX = (xvtx(0) - xpv) / TMath::Sqrt(covX(0, 0));
     213                                Double_t PullY = (xvtx(1) - ypv) / TMath::Sqrt(covX(1, 1));
     214                                Double_t PullZ = (xvtx(2) - zpv) / TMath::Sqrt(covX(2, 2));
    157215                                //
    158216                                // Fill histograms
    … …  
    173231                //
    174232                // Cleanup
    175                 for (Int_t i = 0; i < Ntr; i++) delete pr[i];
    176                 for (Int_t i = 0; i < Ntr; i++) delete cv[i];
     233                for (Int_t i = 0; i < NtrG; i++) delete pr[i];
     234                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];
    177237                delete[] pr;
    178238                delete[] cv;
     239                delete[] PrFit;
     240                delete[] CvFit;
    179241        }
    180242        //
    … …  
    195257        CnvN->cd(1);
    196258        hTrPrim->Draw();
    197         CnvN-> cd(2);
     259        CnvN->cd(2);
    198260        hTrFound->SetLineColor(kRed);
    199261        hTrFound->Draw();

  • examples/ExamplePVtxFit.C

    r4afb18d r7dac4ea  
    1414
    1515#endif
     16
    1617
    1718
    … …  
    4546                // Load selected branches with data from specified event
    4647                treeReader->ReadEntry(entry);
     48                //
    4749                Int_t Ntr = 0;  // # of tracks from primary vertex
    4850                Int_t NtrG = branchTrack->GetEntries();
    4951                TVectorD** pr = new TVectorD * [NtrG];
    5052                TMatrixDSym** cv = new TMatrixDSym * [NtrG];
     53                //
     54                // True vertex
     55                Double_t xpv, ypv, zpv;
    5156                // If event contains at least 1 track
    5257                //
    … …  
    6065                                // Get associated generated particle
    6166                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
     67                                TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
     68                                TVector3 pg(gp->Px,gp->Py,gp->Pz);
     69                                Double_t Q = (Double_t)gp->Charge;
     70                                Double_t Bz = 2.0;
     71                                TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
     72                                TVectorD genPar = TrkUtil::ParToMm(genParM);
    6273                                //
    6374                                // Position of origin in mm
    … …  
    6576                                Double_t y = gp->Y;
    6677                                Double_t z = gp->Z;
     78                                Bool_t prim = kTRUE;    // Is primary?
     79                                Int_t mp = gp->M1;      // Mother
     80                                while(mp>0){
     81                                        GenParticle* gm =
     82                                        (GenParticle*)branchGenPart->At(mp);
     83                                        Double_t xm = gm->X;
     84                                        Double_t ym = gm->Y;
     85                                        Double_t zm = gm->Z;
     86                                        if(x!=xm || y!=ym || z!=zm){
     87                                                prim = kFALSE;
     88                                                break;
     89                                        }else mp = gm->M1;
     90                                }
    6791                                //
    6892                                // group tracks originating from the primary vertex
    69                                 if (x == 0.0 && y == 0.0)
     93                                if (prim)
    7094                                {
     95                                        //std::cout<<"Event: "<<entry<<", track: "<<it;
     96                                        //std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl;
     97                                        xpv = x;
     98                                        ypv = y;
     99                                        zpv = z;
    71100                                        //
    72101                                        // Reconstructed track parameters
    … …  
    75104                                        Double_t obsC = trk->C;
    76105                                        Double_t obsZ0 = trk->DZ;
     106                                        //std::cout<<"Z0 track = "<< obsZ0
     107                                        //<<", gen Z0 = "<<genPar(3)
     108                                        //<<", gen cot = "<<genPar(4)<<std::endl;
    77109                                        Double_t obsCtg = trk->CtgTheta;
    78110                                        Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
    … …  
    96128                        Double_t Chi2 = Vtx->GetVtxChi2();
    97129                        Double_t Ndof = 2 * (Double_t)Ntr - 3;
    98                         Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
    99                         Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
    100                         Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
     130                        Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
     131                        Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
     132                        Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
     133                        //std::cout<<"**** True  vertex (x, y, z) "<<xpv<<", "
     134                        //<<ypv<<", "<<zpv<<std::endl;
     135                        //std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
     136                        //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
    101137                        //
    102138                        // Fill histograms

  • examples/Pythia8/ee_zh.cmnd

    r4afb18d r7dac4ea  
    11
     2! number of events to generate
    23Main:numberOfEvents = 1000         ! number of events to generate
    34
    … …  
    78
    89! Vertex smearing :
    9 !Beams:allowVertexSpread = on
    10 !Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
    11 !Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
    12 !Beams:sigmaVertexZ = 0.64      !  0.64 mm
    13 !Beams:sigmaTime    = 0.64      !  0.64 mm
    14 ï¿Œ
    1510
     11Beams:allowVertexSpread = on
     12Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
     13Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
     14Beams:sigmaVertexZ = 0.64      !  0.64 mm
     15Beams:sigmaTime    = 0.64      !  0.64 mm
    1616
    1717! Higgsstrahlung process
    … …  
    2020! 5) Force the Z decays to muons
    212123:onMode = off
    22 23:onIfAny = 11 -11
     2223:onIfAny = 13 -13

  • external/TrackCovariance/ObsTrk.cc

    r4afb18d r7dac4ea  
    1414ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
    1515{
    16         SetBfield(B);
     16        SetB(B);
    1717        fGC = GC;
    1818        fGenX = x;
    … …  
    5656ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)
    5757{
    58         SetBfield(B);
     58        SetB(B);
    5959        fGC = GC;
    6060        fGenX.SetXYZ(x[0],x[1],x[2]);

  • external/TrackCovariance/VertexFit.cc

    r4afb18d r7dac4ea  
    1818        fxCst.ResizeTo(3);
    1919        fCovCst.ResizeTo(3, 3);
     20        fCovCstInv.ResizeTo(3, 3);
    2021        fXv.ResizeTo(3);
    2122        fcovXv.ResizeTo(3, 3);
    … …  
    3132        fxCst.ResizeTo(3);
    3233        fCovCst.ResizeTo(3, 3);
     34        fCovCstInv.ResizeTo(3, 3);
    3335        fXv.ResizeTo(3);
    3436        fcovXv.ResizeTo(3, 3);
    … …  
    3840                TVectorD pr = *trkPar[i];
    3941                fPar.push_back(new TVectorD(pr));
     42                fParNew.push_back(new TVectorD(pr));
    4043                TMatrixDSym cv = *trkCov[i];
    4144                fCov.push_back(new TMatrixDSym(cv));
     45                fCovNew.push_back(new TMatrixDSym(cv));
    4246        }
    4347        fChi2List.ResizeTo(fNtr);
    … …  
    5458        fxCst.ResizeTo(3);
    5559        fCovCst.ResizeTo(3, 3);
     60        fCovCstInv.ResizeTo(3, 3);
    5661        fXv.ResizeTo(3);
    5762        fcovXv.ResizeTo(3, 3);
    … …  
    6166        {
    6267                fPar.push_back(new TVectorD(track[i]->GetObsPar()));
     68                fParNew.push_back(new TVectorD(track[i]->GetObsPar()));
    6369                fCov.push_back(new TMatrixDSym(track[i]->GetCov()));
     70                fCovNew.push_back(new TMatrixDSym(track[i]->GetCov()));
    6471        }
    6572}
    … …  
    7380        {
    7481                if (fx0i[i])  { fx0i[i]->Clear();               delete fx0i[i]; }
    75                 if (fai[i])   { fai[i]->Clear();                        delete fai[i]; }
    76                 if (fDi[i])   { fDi[i]->Clear();                        delete fDi[i]; }
    77                 if (fWi[i])   { fWi[i]->Clear();                        delete fWi[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]; }
    7887                if (fWinvi[i]){ fWinvi[i]->Clear();     delete fWinvi[i]; }
    7988        }
    … …  
    8190        fx0i.clear();
    8291        fai.clear();
     92        fdi.clear();
     93        fAti.clear();
    8394        fDi.clear();
    8495        fWi.clear();
    … …  
    8899{
    89100        fxCst.Clear();         
    90         fCovCst.Clear();               
     101        fCovCst.Clear();
     102        fCovCstInv.Clear();
    91103        fXv.Clear();           
    92104        fcovXv.Clear();         
    … …  
    95107        for (Int_t i = 0; i < fNtr; i++)
    96108        {
    97                 fPar[i]->Clear();       delete fPar[i];
    98                 fCov[i]->Clear();       delete fCov[i];
     109                fPar[i]->Clear();               delete fPar[i];
     110                fParNew[i]->Clear();    delete fParNew[i];
     111                fCov[i]->Clear();               delete fCov[i];
     112                fCovNew[i]->Clear();    delete fCovNew[i];
    99113        }
    100114        fPar.clear();
    101         fCov.clear();   
     115        fParNew.clear();
     116        fCov.clear();
     117        fCovNew.clear();               
    102118        //
    103119        ResetWrkArrays();
    … …  
    424440        //
    425441        // Get track parameters and their covariance
    426         TVectorD par = *fPar[i];
     442        TVectorD par = *fParNew[i];
    427443        TMatrixDSym Cov = *fCov[i];
    428444        //
    … …  
    433449        //
    434450        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
    435455        TMatrixDSym Winv = Cov.Similarity(A);           // W^-1 = A*C*A'
    436         fWinvi.push_back(new TMatrixDSym(Winv));                // Store W^-1 matrix
     456        fWinvi.push_back(new TMatrixDSym(Winv));        // Store W^-1 matrix
    437457        //
    438458        TMatrixDSym W = RegInv(Winv);                           // W = (A*C*A')^-1
    … …  
    456476{
    457477        //std::cout << "VertexFitter: just in" << std::endl;
    458         if (fNtr < 2)
    459         {
     478        if (fNtr < 2 && !fVtxCst){
    460479                std::cout << "VertexFit::VertexFitter - Method called with less than 2 tracks - Aborting " << std::endl;
    461480                std::exit(1);
    … …  
    469488        Double_t Chi2 = 0;
    470489        TVectorD x0 = fXv;      // If previous fit done
    471         if (fRold < 0.0)for (Int_t i = 0; i < 3; i++)x0(i) = 1000.;     // Set to arbitrary large value if not
     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        }
    472496        //
    473497        // Starting vertex radius approximation
    … …  
    475499        Double_t R = fRold;                                             // Use previous fit if available
    476500        if (R < 0.0) R = StartRadius();                 // Rough vertex estimate
     501        //std::cout << "Start radius: " << R << std::endl;
    477502        //
    478503        // Iteration properties
    … …  
    529554                        H += Ds;
    530555                        // update constant term
    531                         TVectorD xs = *fx0i[i];
     556                        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                        */
    532563                        cterm += Ds * xs;
    533564                }                               // End loop on tracks
     565                // Some additions in case of external constraints
     566                if (fVtxCst) {
     567                        H += fCovCstInv;
     568                        cterm += fCovCstInv * fxCst;
     569                        DW1D += fCovCstInv;
     570                }
    534571                //
    535572                // update vertex position
    … …  
    544581                for (Int_t i = 0; i < fNtr; i++)
    545582                {
    546                         TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x);
     583                        TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x - *fdi[i]);
    547584                        TMatrixDSym Wm1 = *fWinvi[i];
    548585                        fChi2List(i) = Wm1.Similarity(lambda);
    … …  
    551588                        TVectorD b = (*fWi[i]) * (x - (*fx0i[i]));
    552589                        for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i];
     590                        TVectorD newPar = *fPar[i] - ((*fCov[i]) * (*fAti[i])) * lambda;
     591                        fParNew[i] = new TVectorD(newPar);
    553592                }
     593                // Add external constraint to Chi2
     594                if (fVtxCst) Chi2 += fCovCstInv.Similarity(x - fxCst);
    554595                //
    555596                TVectorD dx = x - x0;
    … …  
    562603                // Store result
    563604                //
    564                 fXv = x;                                // Vertex position
     605                fXv = x;                        // Vertex position
    565606                fcovXv = covX;          // Vertex covariance
    566607                fChi2 = Chi2;           // Vertex fit Chi2
     608                //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) << std::endl;
    567609        }               // end of iteration loop
    568610        //
    … …  
    603645void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov)  // Add gaussian vertex constraint
    604646{
    605         std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
     647        //std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
     648        fVtxCst = kTRUE;                                // Vertex constraint flag
     649        fxCst = xv;                                             // Constraint value
     650        fCovCst = cov;                                  // Constraint covariance
     651        fCovCstInv = cov;
     652        fCovCstInv.Invert();                            // Its inverse
     653        //
     654        // Set starting vertex as external constraint
     655        fXv = fxCst;
     656        fcovXv = fCovCst;
    606657}
    607658//
    … …  
    613664        fPar.push_back(par);                    // add new track
    614665        fCov.push_back(Cov);
     666        fParNew.push_back(par);                 // add new track
     667        fCovNew.push_back(Cov);
    615668        //
    616669        // Reset previous vertex temp arrays
    … …  
    628681        fPar.erase(fPar.begin() + iTrk);                // Remove track
    629682        fCov.erase(fCov.begin() + iTrk);
     683        fParNew.erase(fParNew.begin() + iTrk);          // Remove track
     684        fCovNew.erase(fCovNew.begin() + iTrk);
    630685        //
    631686        // Reset previous vertex temp arrays

  • external/TrackCovariance/VertexFit.h

    r4afb18d r7dac4ea  
    2323        Int_t fNtr;                                     // Number of tracks
    2424        std::vector<TVectorD*> fPar;            // Input parameter array
    25         std::vector<TMatrixDSym*> fCov;// Input parameter covariances
     25        std::vector<TVectorD*> fParNew;         // Updated parameter array
     26        std::vector<TMatrixDSym*> fCov;         // Input parameter covariances
     27        std::vector<TMatrixDSym*> fCovNew;      // Updated parameter covariances
    2628        // Constraints
    2729        Bool_t fVtxCst;                         // Vertex constraint flag
    2830        TVectorD fxCst;                         // Constraint value
    29         TMatrixDSym fCovCst;                    // Constraint covariance
     31        TMatrixDSym fCovCst;                    // Constraint
     32        TMatrixDSym fCovCstInv;                 // Inverse of constraint covariance
    3033        //
    3134        // Results
    … …  
    3942        // Work arrays
    4043        std::vector<Double_t> ffi;                              // Fit phases
    41         std::vector<TVectorD*> fx0i;                            // Track expansion points
     44        std::vector<TVectorD*> fx0i;                    // Track expansion points
    4245        std::vector<TVectorD*> fai;                             // dx/dphi
     46        std::vector<TVectorD*> fdi;                             // x-shift
    4347        std::vector<Double_t> fa2i;                             // a'Wa
     48        std::vector<TMatrixD*> fAti;                    // A transposed
    4449        std::vector<TMatrixDSym*> fDi;                  // W-WBW
    4550        std::vector<TMatrixDSym*> fWi;                  // (ACA')^-1
Note: See TracChangeset for help on using the changeset viewer.