Changes in examples/ExamplePVtxFitEC.C [b750b0a:46d3442] in git
- File:
-
- 1 edited
-
examples/ExamplePVtxFitEC.C (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
examples/ExamplePVtxFitEC.C
rb750b0a r46d3442 1 1 /* 2 2 Example of using vertex fitter class to fit primary vertex 3 assumed to be generated with Pythia8/ee_zh_smr-shf.cmn3 assumed to be generated in (0,0,0) 4 4 */ 5 5 … … 21 21 void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5) 22 22 { 23 //24 // Beam constraint25 // (consistent with generation with ee_zh_smr-shf.cmnd)26 //27 // Mean beam position28 TVectorD xpvc(3);29 xpvc(0) = 1.0;30 xpvc(1) = -2.0;31 xpvc(2) = 10.0;32 // Interaction region covariance33 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 //39 23 // Create chain of root trees 40 24 TChain chain("Delphes"); … … 51 35 // Book histograms 52 36 Int_t Nbin = 100; 53 // Vertex fit pulls54 37 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.); 55 38 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.); 56 39 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.); 57 40 TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.); 58 // Generation check59 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 63 41 // 64 42 // Loop over all events … … 87 65 // Get associated generated particle 88 66 GenParticle* gp = (GenParticle*)trk->Particle.GetObject(); 89 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); // mm -> meters67 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); 90 68 TVector3 pg(gp->Px,gp->Py,gp->Pz); 91 69 Double_t Q = (Double_t)gp->Charge; 92 70 Double_t Bz = 2.0; 93 71 TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz); 94 TVectorD genPar = TrkUtil::ParToMm(genParM); // -> back to mm 95 96 // 97 // Check if original track is from primary vertex 72 TVectorD genPar = TrkUtil::ParToMm(genParM); 98 73 // 99 74 // Position of origin in mm … … 104 79 Int_t mp = gp->M1; // Mother 105 80 while(mp>0){ 106 GenParticle* gm = (GenParticle*)branchGenPart->At(mp); 81 GenParticle* gm = 82 (GenParticle*)branchGenPart->At(mp); 107 83 Double_t xm = gm->X; 108 84 Double_t ym = gm->Y; 109 85 Double_t zm = gm->Z; 110 86 if(x!=xm || y!=ym || z!=zm){ 111 prim = kFALSE; // Not primary87 prim = kFALSE; 112 88 break; 113 89 }else mp = gm->M1; 114 90 } 115 116 91 // 117 92 // group tracks originating from the primary vertex 118 93 if (prim) 119 94 { 120 // 121 // Store true primary vertex for this event95 //std::cout<<"Event: "<<entry<<", track: "<<it; 96 //std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl; 122 97 xpv = x; 123 98 ypv = y; … … 129 104 Double_t obsC = trk->C; 130 105 Double_t obsZ0 = trk->DZ; 106 //std::cout<<"Z0 track = "<< obsZ0 107 //<<", gen Z0 = "<<genPar(3) 108 //<<", gen cot = "<<genPar(4)<<std::endl; 131 109 Double_t obsCtg = trk->CtgTheta; 132 110 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg }; … … 134 112 // 135 113 pr[Ntr] = new TVectorD(obsPar); 114 //std::cout<<"Cov Matrix:"<<std::endl; 115 //trk->CovarianceMatrix().Print(); 136 116 cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix()); 137 117 Ntr++; … … 142 122 // Fit primary vertex with beam constraint 143 123 // 144 Int_t MinTrk = 3; // Minumum # tracks for vertex fit 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 145 134 if (Ntr >= MinTrk) { 146 135 VertexFit* Vtx = new VertexFit(Ntr, pr, cv); … … 150 139 Double_t Chi2 = Vtx->GetVtxChi2(); 151 140 Double_t Ndof = 2 * (Double_t)Ntr; 152 delete Vtx;153 //154 141 Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0)); 155 142 Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1)); 156 143 Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2)); 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 // 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; 162 148 // 163 149 // Fill histograms … … 166 152 hZpull->Fill(PullZ); 167 153 hChi2->Fill(Chi2 / Ndof); 168 //169 hXvpull->Fill(PullXv);170 hYvpull->Fill(PullYv);171 hZvpull->Fill(PullZv);172 154 } 173 155 156 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl; 174 157 // 175 158 // Cleanup … … 179 162 delete[] cv; 180 163 } 181 182 164 // 183 165 // Show resulting histograms … … 192 174 hZpull->Fit("gaus"); hZpull->Draw(); 193 175 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();203 176 }
Note:
See TracChangeset
for help on using the changeset viewer.
