Fork me on GitHub

Changeset dd263e4 in git for examples/ExamplePVtxFitEC.C


Ignore:
Timestamp:
Jan 19, 2022, 12:08:01 PM (3 years ago)
Author:
GitHub <noreply@…>
Branches:
master
Children:
2c8865f, 56fb0be, 5c03893
Parents:
0c0c9af (diff), 78ce8d1 (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@…> (01/19/22 12:08:01)
git-committer:
GitHub <noreply@…> (01/19/22 12:08:01)
Message:

Merge pull request #104 from fbedesch/master

CovarianceMatrix scaling fir and vertexing improvements

File:
1 edited

Legend:

Unmodified
Added
Removed

  • examples/ExamplePVtxFitEC.C

    r0c0c9af rdd263e4  
    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.