Changeset b750b0a in git for examples/ExamplePVtxFitEC.C
- Timestamp:
- Jan 19, 2022, 10:11:39 AM (3 years ago)
- Branches:
- master
- Children:
- 78ce8d1, 7bca620
- Parents:
- bd376e3
- File:
-
- 1 edited
-
examples/ExamplePVtxFitEC.C (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
examples/ExamplePVtxFitEC.C
rbd376e3 rb750b0a 1 1 /* 2 2 Example of using vertex fitter class to fit primary vertex 3 assumed to be generated in (0,0,0)3 assumed to be generated with Pythia8/ee_zh_smr-shf.cmn 4 4 */ 5 5 … … 21 21 void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5) 22 22 { 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 // 23 39 // Create chain of root trees 24 40 TChain chain("Delphes"); … … 35 51 // Book histograms 36 52 Int_t Nbin = 100; 53 // Vertex fit pulls 37 54 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.); 38 55 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.); 39 56 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.); 40 57 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 41 63 // 42 64 // Loop over all events … … 65 87 // Get associated generated particle 66 88 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 68 90 TVector3 pg(gp->Px,gp->Py,gp->Pz); 69 91 Double_t Q = (Double_t)gp->Charge; 70 92 Double_t Bz = 2.0; 71 93 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 73 98 // 74 99 // Position of origin in mm … … 79 104 Int_t mp = gp->M1; // Mother 80 105 while(mp>0){ 81 GenParticle* gm = 82 (GenParticle*)branchGenPart->At(mp); 106 GenParticle* gm = (GenParticle*)branchGenPart->At(mp); 83 107 Double_t xm = gm->X; 84 108 Double_t ym = gm->Y; 85 109 Double_t zm = gm->Z; 86 110 if(x!=xm || y!=ym || z!=zm){ 87 prim = kFALSE; 111 prim = kFALSE; // Not primary 88 112 break; 89 113 }else mp = gm->M1; 90 114 } 115 91 116 // 92 117 // group tracks originating from the primary vertex 93 118 if (prim) 94 119 { 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 97 122 xpv = x; 98 123 ypv = y; … … 104 129 Double_t obsC = trk->C; 105 130 Double_t obsZ0 = trk->DZ; 106 //std::cout<<"Z0 track = "<< obsZ0107 //<<", gen Z0 = "<<genPar(3)108 //<<", gen cot = "<<genPar(4)<<std::endl;109 131 Double_t obsCtg = trk->CtgTheta; 110 132 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg }; … … 112 134 // 113 135 pr[Ntr] = new TVectorD(obsPar); 114 //std::cout<<"Cov Matrix:"<<std::endl;115 //trk->CovarianceMatrix().Print();116 136 cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix()); 117 137 Ntr++; … … 122 142 // Fit primary vertex with beam constraint 123 143 // 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 134 145 if (Ntr >= MinTrk) { 135 146 VertexFit* Vtx = new VertexFit(Ntr, pr, cv); … … 139 150 Double_t Chi2 = Vtx->GetVtxChi2(); 140 151 Double_t Ndof = 2 * (Double_t)Ntr; 152 delete Vtx; 153 // 141 154 Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0)); 142 155 Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1)); 143 156 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 // 148 162 // 149 163 // Fill histograms … … 152 166 hZpull->Fill(PullZ); 153 167 hChi2->Fill(Chi2 / Ndof); 168 // 169 hXvpull->Fill(PullXv); 170 hYvpull->Fill(PullYv); 171 hZvpull->Fill(PullZv); 154 172 } 155 173 156 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;157 174 // 158 175 // Cleanup … … 162 179 delete[] cv; 163 180 } 181 164 182 // 165 183 // Show resulting histograms … … 174 192 hZpull->Fit("gaus"); hZpull->Draw(); 175 193 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(); 176 203 }
Note:
See TracChangeset
for help on using the changeset viewer.
