Changeset a1a25d4 in git
- Timestamp:
- Mar 2, 2021, 3:38:05 PM (4 years ago)
- Branches:
- master
- Children:
- 4586f8f
- Parents:
- d3165fa (diff), 92b8d11 (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@…> (03/02/21 15:38:05)
- git-committer:
- GitHub <noreply@…> (03/02/21 15:38:05)
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
examples/Example6.C
rd3165fa ra1a25d4 51 51 // Track parameter pulls 52 52 TH1* histDpull = new TH1F("h_Dpull", "Pull impact parameter", 100, -10., 10.); 53 TH1* histP0pull = new TH1F("h_P0pull", "Pull phi_{0}", 100, -10., 10.);53 TH1* histP0pull = new TH1F("h_P0pull", "Pull #phi_{0}", 100, -10., 10.); 54 54 TH1* histCpull = new TH1F("h_Cpull", "Pull half curvature", 100, -10., 10.); 55 TH1* histZ0pull = new TH1F("h_Z0pull", "Pull Z 0", 100, -10., 10.);55 TH1* histZ0pull = new TH1F("h_Z0pull", "Pull Z_{0}", 100, -10., 10.); 56 56 TH1* histCtpull = new TH1F("h_Ctpull", "Pull cotg(#theta)", 100, -10., 10.); 57 57 // … … 64 64 // CEntral track over min Pt 65 65 TH1* histPtgenC = new TH1F("h_PtgenC", "Generated Pt - Central", 500, 0., 50.); // pt for central tracks; 66 TH1* histPtobsC = new TH1F("h_PtobsC", "Reconstructed Pt - Central", 500, 0., 50.); // pt for central 66 TH1* histPtobsC = new TH1F("h_PtobsC", "Reconstructed Pt - Central", 500, 0., 50.); // pt for central 67 67 TH1* histCtgenH = new TH1F("h_CtgenH", "Generateded Cotangent", 100, -10., 10.); // cot(theta) for high pt tracks; 68 68 TH1* histCtobsH = new TH1F("h_CtobsH", "Reconstructed Cotangent", 100, -10., 10.); // cot(theta) for high pt tracks … … 71 71 TH1* histAccPt = new TH1F("h_AccPt", "Pt acceptance", 500, 0., 50.); 72 72 // Tracks over min pt 73 TH1* histAccCtgH = new TH1F("h_AccCtgH", "Cotangent acceptance (high pt)", 500, -10., 10.);73 TH1* histAccCtgH = new TH1F("h_AccCtgH", "Cotangent acceptance (high pt)", 100, -10., 10.); 74 74 // tracks above min angle 75 75 TH1* histAccPtC = new TH1F("h_AccPtC", "Pt acceptance (central tracks)", 500, 0., 50.); 76 76 77 77 // 78 // Get magneti c field78 // Get magnetifc field 79 79 Double_t Bz = 2.0; 80 80 … … 106 106 histZ0obs->Fill(obsZ0); 107 107 histCtobs->Fill(obsCtg); 108 108 // 109 // Get associated generated particle 109 110 GenParticle* gp = (GenParticle*)trk->Particle.GetObject(); 110 111 // 111 112 // Position of origin in meters 112 113 Double_t x = 1.0e-3 * gp->X; … … 139 140 Double_t pullZ0 = (obsZ0 - genZ0) / trk->ErrorDZ; 140 141 Double_t pullCtg = (obsCtg - genCtg) / trk->ErrorCtgTheta; 141 142 // 143 histDpull->Fill(pullD0); 142 // 143 histDpull ->Fill(pullD0); 144 144 histP0pull->Fill(pullPhi); 145 histCpull ->Fill(pullC);145 histCpull ->Fill(pullC); 146 146 histZ0pull->Fill(pullZ0); 147 147 histCtpull->Fill(pullCtg); … … 149 149 // Acceptance plots 150 150 Double_t obsPt = trk->PT; 151 histPtobs->Fill(obsPt); 151 152 if (TMath::Abs(genCtg) < maxCtAcc)histPtobsC->Fill(obsPt); // pt for central tracks 152 if ( genPt > minPtAcc)histCtobsH->Fill(obsCtg); // cot(theta) for high pt tracks153 if (obsPt > minPtAcc)histCtobsH->Fill(obsCtg); // cot(theta) for high pt tracks 153 154 } // End loop on tracks 154 155 … … 161 162 GenParticle* gpart = (GenParticle*)branchGenPart->At(it); 162 163 // 163 // Plot charged particle parameters 164 // Plot charged particle parameters 164 165 // Only final state particles (Status = 1) 165 166 if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) { 166 167 // 167 // Position of origin 168 Double_t x = gpart->X;169 Double_t y = gpart->Y;170 Double_t z = gpart->Z;168 // Position of origin in meters 169 Double_t x = 1.0e-3 * gpart->X; 170 Double_t y = 1.0e-3 * gpart->Y; 171 Double_t z = 1.0e-3 * gpart->Z; 171 172 TVector3 xv(x, y, z); 172 173 // … … 179 180 // Original parameters 180 181 Double_t Q = gpart->Charge; 181 TVectorD genPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum 182 TVectorD gPar = TrkUtil::XPtoPar(xv, pv, Q, Bz); // Get parameters from position and momentum 183 TVectorD genPar = TrkUtil::ParToMm(gPar); // Convert to mm 182 184 // Unpack parameters 183 Double_t genD0 = genPar(0);185 Double_t genD0 = genPar(0); 184 186 Double_t genPhi = genPar(1); 185 Double_t genC = genPar(2);186 Double_t genZ0 = genPar(3);187 Double_t genC = genPar(2); 188 Double_t genZ0 = genPar(3); 187 189 Double_t genCtg = genPar(4); 190 Double_t genPt = pv.Pt(); 188 191 // Fill histograms 189 histDgen ->Fill(genD0);192 histDgen ->Fill(genD0); 190 193 histP0gen->Fill(genPhi); 191 histCgen ->Fill(genC);194 histCgen ->Fill(genC); 192 195 histZ0gen->Fill(genZ0); 193 196 histCtgen->Fill(genCtg); 194 197 // 195 198 // Acceptance plots 196 if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(gpart->PT); 197 if (gpart->PT > minPtAcc)histCtgenH->Fill(genCtg); 199 histPtgen->Fill(genPt); 200 if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(genPt); 201 if (genPt > minPtAcc)histCtgenH->Fill(genCtg); 198 202 } 199 203 } … … 222 226 TCanvas* Cnv2 = new TCanvas("Cnv2", "Delphes observed track pulls", 150, 150, 900, 500); 223 227 Cnv2->Divide(3, 2); 224 Cnv2->cd(1); histDpull->Draw(); 225 Cnv2->cd(2); histP0pull->Draw(); 226 Cnv2->cd(3); histCpull->Draw(); 227 Cnv2->cd(4); histZ0pull->Draw(); 228 Cnv2->cd(5); histCtpull->Draw(); 228 Cnv2->cd(1); gPad->SetLogy(1); 229 histDpull->Draw(); 230 Cnv2->cd(2); gPad->SetLogy(1); 231 histP0pull->Draw(); 232 Cnv2->cd(3); gPad->SetLogy(1); 233 histCpull->Draw(); 234 Cnv2->cd(4); gPad->SetLogy(1); 235 histZ0pull->Draw(); 236 Cnv2->cd(5); gPad->SetLogy(1); 237 histCtpull->Draw(); 229 238 // 230 239 // Acceptance plots 231 240 // 232 TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 350, 350, 500, 500);241 TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 200, 200, 500, 500); 233 242 Cnv4->Divide(2, 2); 234 243 Cnv4->cd(1); 235 244 histCtgen->Draw(); 245 histCtobs->SetLineColor(kRed); 236 246 histCtobs->Draw("SAME"); 237 247 Cnv4->cd(2); … … 246 256 } 247 257 histAccCtg->Draw(); 248 Cnv4->cd(3); 258 Cnv4->cd(3); gPad->SetLogy(1); 249 259 histPtgen->Draw(); 260 histPtobs->SetLineColor(kRed); 250 261 histPtobs->Draw("SAME"); 251 262 Cnv4->cd(4); … … 261 272 histAccPt->Draw(); 262 273 // 263 TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 400, 400, 500, 500);274 TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 250, 250, 500, 500); 264 275 Cnv5->Divide(2, 2); 265 276 Cnv5->cd(1); 266 277 histCtgenH->Draw(); 278 histCtobsH->SetLineColor(kRed); 267 279 histCtobsH->Draw("SAME"); 268 280 Cnv5->cd(2); … … 277 289 } 278 290 histAccCtgH->Draw(); 279 Cnv5->cd(3); 280 histPtgenC->Draw(); 291 Cnv5->cd(3); gPad->SetLogy(1); 292 histPtgenC->Draw(); 293 histPtobsC->SetLineColor(kRed); 281 294 histPtobsC->Draw("SAME"); 282 295 Cnv5->cd(4); -
external/TrackCovariance/TrkUtil.cc
rd3165fa ra1a25d4 1 1 #include "TrkUtil.h" 2 #include <iostream> 2 3 3 4 // Constructor … … 31 32 Double_t cross = x(0) * p(1) - x(1) * p(0); 32 33 Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2); 33 Double_t phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi034 Double_t phi0 = atan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi0 34 35 Double_t D; // Impact parameter D 35 36 if (pt < 10.0) D = (T - pt) / a; … … 41 42 //Longitudinal parameters 42 43 Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D)); 43 Double_t st = TMath::ASin(B) / C;44 Double_t st = asin(B) / C; 44 45 Double_t ct = p(2) / pt; 45 46 Double_t z0 = x(2) - ct * st; … … 68 69 // 69 70 TVector3 Xval; 70 Xval(0) = -D * TMath::Sin(phi0);71 Xval(1) = D * TMath::Cos(phi0);71 Xval(0) = -D * sin(phi0); 72 Xval(1) = D * cos(phi0); 72 73 Xval(2) = z0; 73 74 // … … 76 77 // 77 78 TVector3 TrkUtil::ParToP(TVectorD Par) 79 { 80 if (fBz == 0.0) 81 std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 82 // 83 return ParToP(Par,fBz); 84 } 85 // 86 TVector3 TrkUtil::ParToP(TVectorD Par, Double_t Bz) 78 87 { 79 88 Double_t C = Par(2); … … 82 91 // 83 92 TVector3 Pval; 84 Double_t pt = fBz * cSpeed() / TMath::Abs(2 * C);85 Pval(0) = pt * TMath::Cos(phi0);86 Pval(1) = pt * TMath::Sin(phi0);93 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C); 94 Pval(0) = pt * cos(phi0); 95 Pval(1) = pt * sin(phi0); 87 96 Pval(2) = pt * ct; 88 97 // … … 103 112 Double_t b = -cSpeed() * fBz / 2.; 104 113 pACTS(0) = 1000 * Par(0); // D from m to mm 105 pACTS(1) = 1000 * Par(3); // z0 from m to mm114 pACTS(1) = 1000 * Par(3); // z0 from m to mm 106 115 pACTS(2) = Par(1); // Phi0 is unchanged 107 pACTS(3) = TMath::ATan2(1.0, Par(4)); // Theta in [0, pi] range116 pACTS(3) = atan2(1.0, Par(4)); // Theta in [0, pi] range 108 117 pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4))); // q/p in GeV 109 118 pACTS(5) = 0.0; // Time: currently undefined -
external/TrackCovariance/TrkUtil.h
rd3165fa ra1a25d4 20 20 void SetBfield(Double_t Bz) { fBz = Bz; } 21 21 TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q); 22 TVector3 ParToX(TVectorD Par);23 22 TVector3 ParToP(TVectorD Par); 24 Double_t ParToQ(TVectorD Par);25 23 // 26 24 // Conversion to ACTS parametrization … … 46 44 static Double_t cSpeed() 47 45 { 48 return TMath::C()*1.0e-9; // Reduced speed of light 46 Double_t c = 2.99792458e8; // speed of light m/sec 47 //return TMath::C()*1.0e-9; // Incompatible with root5 48 return c*1.0e-9; // Reduced speed of light 49 49 } 50 50 // … … 52 52 // 53 53 static TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz); 54 static TVector3 ParToX(TVectorD Par); // position of minimum distance from z axis 55 static TVector3 ParToP(TVectorD Par, Double_t Bz); // Get Momentum from track parameters 56 static Double_t ParToQ(TVectorD Par); // Get track charge 54 57 // 55 58 // Conversion from meters to mm
Note:
See TracChangeset
for help on using the changeset viewer.