Fork me on GitHub

Changeset b750b0a in git for examples


Ignore:
Timestamp:
Jan 19, 2022, 10:11:39 AM (3 years ago)
Author:
Franco BEDESCHI <bed@…>
Branches:
master
Children:
78ce8d1, 7bca620
Parents:
bd376e3
Message:

Fix CovarianceMatrix scaling error. Vertexing improvements.

Location:
examples
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • examples/ExamplePVtxFind.C

    rbd376e3 rb750b0a  
    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

    rbd376e3 rb750b0a  
    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}
Note: See TracChangeset for help on using the changeset viewer.