- Timestamp:
- May 17, 2021, 6:05:45 PM (3 years ago)
- 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)
- Location:
- examples
- Files:
-
- 2 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
examples/ExamplePVtxFind.C
r4afb18d r7dac4ea 51 51 // Load selected branches with data from specified event 52 52 treeReader->ReadEntry(entry); 53 Int_t Ntr = 0; // # of starting tracks used in primary vertex54 53 Int_t NtrG = branchTrack->GetEntries(); 55 //std::cout << "Eventopened containing " << NtrG << " tracks" << std::endl;54 if(entry%500 ==0)std::cout << "Event "<<entry<<" opened containing " << NtrG << " tracks" << std::endl; 56 55 TVectorD** pr = new TVectorD * [NtrG]; 57 56 TMatrixDSym** cv = new TMatrixDSym * [NtrG]; … … 63 62 // 64 63 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; 65 67 if (branchTrack->GetEntries() > 0) 66 68 { … … 80 82 //std::cout << "Got track parameters for track " << it << std::endl; 81 83 // 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 94 92 GenParticle* gp = (GenParticle*)trk->Particle.GetObject(); 95 93 //std::cout << "GenParticle pointer "<<gp << std::endl; … … 99 97 Double_t y = gp->Y; 100 98 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 } 105 119 106 120 } // End loop on tracks 107 121 } 108 122 // 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); 117 171 //std::cout << "Vertex fit created " << std::endl; 172 Vtx->AddVtxConstraint(xpvc, covpvc); 118 173 // 119 174 // Remove tracks with large chi2 175 Double_t MaxChi2Fit = 9.; 120 176 Bool_t Done = kFALSE; 121 while (!Done) {177 while (!Done) { 122 178 //std::cout << "After while " << std::endl; 123 179 // Find largest Chi2 contribution … … 130 186 Double_t Chi2Mx = Chi2L[iMax]; 131 187 //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; 134 192 Nfound--; 135 193 } … … 145 203 // 146 204 // Require minimum number of tracks in vertex 147 Int_t Nmin = 4;205 Int_t Nmin = 2; 148 206 if (Nfound >= Nmin) { 149 207 TVectorD xvtx = Vtx->GetVtx(); … … 151 209 TMatrixDSym covX = Vtx->GetVtxCov(); 152 210 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)); 157 215 // 158 216 // Fill histograms … … 173 231 // 174 232 // 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]; 177 237 delete[] pr; 178 238 delete[] cv; 239 delete[] PrFit; 240 delete[] CvFit; 179 241 } 180 242 // … … 195 257 CnvN->cd(1); 196 258 hTrPrim->Draw(); 197 CnvN-> 259 CnvN->cd(2); 198 260 hTrFound->SetLineColor(kRed); 199 261 hTrFound->Draw(); -
examples/ExamplePVtxFit.C
r4afb18d r7dac4ea 14 14 15 15 #endif 16 16 17 17 18 … … 45 46 // Load selected branches with data from specified event 46 47 treeReader->ReadEntry(entry); 48 // 47 49 Int_t Ntr = 0; // # of tracks from primary vertex 48 50 Int_t NtrG = branchTrack->GetEntries(); 49 51 TVectorD** pr = new TVectorD * [NtrG]; 50 52 TMatrixDSym** cv = new TMatrixDSym * [NtrG]; 53 // 54 // True vertex 55 Double_t xpv, ypv, zpv; 51 56 // If event contains at least 1 track 52 57 // … … 60 65 // Get associated generated particle 61 66 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); 62 73 // 63 74 // Position of origin in mm … … 65 76 Double_t y = gp->Y; 66 77 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 } 67 91 // 68 92 // group tracks originating from the primary vertex 69 if ( x == 0.0 && y == 0.0)93 if (prim) 70 94 { 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; 71 100 // 72 101 // Reconstructed track parameters … … 75 104 Double_t obsC = trk->C; 76 105 Double_t obsZ0 = trk->DZ; 106 //std::cout<<"Z0 track = "<< obsZ0 107 //<<", gen Z0 = "<<genPar(3) 108 //<<", gen cot = "<<genPar(4)<<std::endl; 77 109 Double_t obsCtg = trk->CtgTheta; 78 110 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg }; … … 96 128 Double_t Chi2 = Vtx->GetVtxChi2(); 97 129 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; 101 137 // 102 138 // Fill histograms -
examples/Pythia8/ee_zh.cmnd
r4afb18d r7dac4ea 1 1 2 ! number of events to generate 2 3 Main:numberOfEvents = 1000 ! number of events to generate 3 4 … … 7 8 8 9 ! Vertex smearing : 9 !Beams:allowVertexSpread = on10 !Beams:sigmaVertexX = 9.70e-3 ! 13.7 mum / sqrt211 !Beams:sigmaVertexY = 25.5E-6 ! 36.1 nm / sqrt212 !Beams:sigmaVertexZ = 0.64 ! 0.64 mm13 !Beams:sigmaTime = 0.64 ! 0.64 mm14 ï¿Œ15 10 11 Beams:allowVertexSpread = on 12 Beams:sigmaVertexX = 9.70e-3 ! 13.7 mum / sqrt2 13 Beams:sigmaVertexY = 25.5E-6 ! 36.1 nm / sqrt2 14 Beams:sigmaVertexZ = 0.64 ! 0.64 mm 15 Beams:sigmaTime = 0.64 ! 0.64 mm 16 16 17 17 ! Higgsstrahlung process … … 20 20 ! 5) Force the Z decays to muons 21 21 23:onMode = off 22 23:onIfAny = 1 1 -1122 23:onIfAny = 13 -13
Note:
See TracChangeset
for help on using the changeset viewer.