Changeset 341014c in git for modules/TrackSmearing.cc
- Timestamp:
- Feb 12, 2019, 9:29:17 PM (6 years ago)
- Branches:
- ImprovedOutputFile, Timing, llp, master
- Children:
- 6455202
- Parents:
- 45e58be
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/TrackSmearing.cc
r45e58be r341014c 13 13 #include "classes/DelphesFormula.h" 14 14 15 #include "ExRootAnalysis/ExRootClassifier.h" 16 #include "ExRootAnalysis/ExRootFilter.h" 15 17 #include "ExRootAnalysis/ExRootResult.h" 16 #include "ExRootAnalysis/ExRootFilter.h" 17 #include "ExRootAnalysis/ExRootClassifier.h" 18 18 19 #include "TDatabasePDG.h" 20 #include "TFile.h" 21 #include "TFormula.h" 22 #include "TLorentzVector.h" 19 23 #include "TMath.h" 24 #include "TObjArray.h" 25 #include "TProfile2D.h" 26 #include "TRandom3.h" 20 27 #include "TString.h" 21 #include "TFormula.h" 22 #include "TRandom3.h" 23 #include "TObjArray.h" 24 #include "TDatabasePDG.h" 25 #include "TLorentzVector.h" 26 #include "TFile.h" 27 #include "TProfile2D.h" 28 29 #include <algorithm> 30 #include <stdexcept> 28 29 #include <algorithm> 31 30 #include <iostream> 32 31 #include <sstream> 32 #include <stdexcept> 33 33 34 34 using namespace std; … … 66 66 67 67 // !!! IF WE WANT TO KEEP ROOT INPUT !!! 68 if (string(GetString("D0ResolutionFormula", "0.0")) != "0.0")69 70 71 72 73 else 74 75 76 77 78 79 if (string(GetString("DZResolutionFormula", "0.0")) != "0.0")80 81 82 83 84 else 85 86 87 88 89 90 if (string(GetString("PResolutionFormula", "0.0")) != "0.0")91 92 93 94 95 else 96 97 98 99 100 101 if (string(GetString("CtgThetaResolutionFormula", "0.0")) != "0.0")102 103 104 105 106 else 107 108 109 110 111 112 if (string(GetString("PhiResolutionFormula", "0.0")) != "0.0")113 114 115 116 117 else 118 119 120 121 122 68 if(string(GetString("D0ResolutionFormula", "0.0")) != "0.0") 69 { 70 fD0Formula->Compile(GetString("D0ResolutionFormula", "0.0")); 71 fUseD0Formula = true; 72 } 73 else 74 { 75 fD0ResolutionFile = GetString("D0ResolutionFile", "errors.root"); 76 fD0ResolutionHist = GetString("D0ResolutionHist", "d0"); 77 fUseD0Formula = false; 78 } 79 if(string(GetString("DZResolutionFormula", "0.0")) != "0.0") 80 { 81 fDZFormula->Compile(GetString("DZResolutionFormula", "0.0")); 82 fUseDZFormula = true; 83 } 84 else 85 { 86 fDZResolutionFile = GetString("DZResolutionFile", "errors.root"); 87 fDZResolutionHist = GetString("DZResolutionHist", "dz"); 88 fUseDZFormula = false; 89 } 90 if(string(GetString("PResolutionFormula", "0.0")) != "0.0") 91 { 92 fPFormula->Compile(GetString("PResolutionFormula", "0.0")); 93 fUsePFormula = true; 94 } 95 else 96 { 97 fPResolutionFile = GetString("PResolutionFile", "errors.root"); 98 fPResolutionHist = GetString("PResolutionHist", "p"); 99 fUsePFormula = false; 100 } 101 if(string(GetString("CtgThetaResolutionFormula", "0.0")) != "0.0") 102 { 103 fCtgThetaFormula->Compile(GetString("CtgThetaResolutionFormula", "0.0")); 104 fUseCtgThetaFormula = true; 105 } 106 else 107 { 108 fCtgThetaResolutionFile = GetString("CtgThetaResolutionFile", "errors.root"); 109 fCtgThetaResolutionHist = GetString("CtgThetaResolutionHist", "ctgTheta"); 110 fUseCtgThetaFormula = false; 111 } 112 if(string(GetString("PhiResolutionFormula", "0.0")) != "0.0") 113 { 114 fPhiFormula->Compile(GetString("PhiResolutionFormula", "0.0")); 115 fUsePhiFormula = true; 116 } 117 else 118 { 119 fPhiResolutionFile = GetString("PhiResolutionFile", "errors.root"); 120 fPhiResolutionHist = GetString("PhiResolutionHist", "phi"); 121 fUsePhiFormula = false; 122 } 123 123 124 124 fApplyToPileUp = GetBool("ApplyToPileUp", true); … … 128 128 fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/stableParticles")); 129 129 fItInputArray = fInputArray->MakeIterator(); 130 130 131 131 // import beamspot 132 132 try … … 137 137 { 138 138 fBeamSpotInputArray = 0; 139 } 140 139 } 140 141 141 // create output array 142 142 … … 170 170 *phiErrorHist = NULL; 171 171 172 if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0)172 if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0) 173 173 beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0); 174 174 else 175 175 { 176 Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At(0));176 Candidate &beamSpotCandidate = *((Candidate *)fBeamSpotInputArray->At(0)); 177 177 beamSpotPosition = beamSpotCandidate.Position; 178 178 } 179 179 180 if 181 { 182 TFile *fin = TFile::Open (fD0ResolutionFile.c_str());183 d0ErrorHist = (TProfile2D *) fin->Get (fD0ResolutionHist.c_str());184 d0ErrorHist->SetDirectory(0);185 fin->Close();186 } 187 if 188 { 189 TFile *fin = TFile::Open (fDZResolutionFile.c_str());190 dzErrorHist = (TProfile2D *) fin->Get (fDZResolutionHist.c_str());191 dzErrorHist->SetDirectory(0);192 fin->Close();193 } 194 if 195 { 196 TFile *fin = TFile::Open (fPResolutionFile.c_str());197 pErrorHist = (TProfile2D *) fin->Get (fPResolutionHist.c_str());198 pErrorHist->SetDirectory(0);199 fin->Close();200 } 201 if 202 { 203 TFile *fin = TFile::Open (fCtgThetaResolutionFile.c_str());204 ctgThetaErrorHist = (TProfile2D *) fin->Get (fCtgThetaResolutionHist.c_str());205 ctgThetaErrorHist->SetDirectory(0);206 fin->Close();207 } 208 if 209 { 210 TFile *fin = TFile::Open (fPhiResolutionFile.c_str());211 phiErrorHist = (TProfile2D *) fin->Get (fPhiResolutionHist.c_str());212 phiErrorHist->SetDirectory(0);213 fin->Close();180 if(!fUseD0Formula) 181 { 182 TFile *fin = TFile::Open(fD0ResolutionFile.c_str()); 183 d0ErrorHist = (TProfile2D *)fin->Get(fD0ResolutionHist.c_str()); 184 d0ErrorHist->SetDirectory(0); 185 fin->Close(); 186 } 187 if(!fUseDZFormula) 188 { 189 TFile *fin = TFile::Open(fDZResolutionFile.c_str()); 190 dzErrorHist = (TProfile2D *)fin->Get(fDZResolutionHist.c_str()); 191 dzErrorHist->SetDirectory(0); 192 fin->Close(); 193 } 194 if(!fUsePFormula) 195 { 196 TFile *fin = TFile::Open(fPResolutionFile.c_str()); 197 pErrorHist = (TProfile2D *)fin->Get(fPResolutionHist.c_str()); 198 pErrorHist->SetDirectory(0); 199 fin->Close(); 200 } 201 if(!fUseCtgThetaFormula) 202 { 203 TFile *fin = TFile::Open(fCtgThetaResolutionFile.c_str()); 204 ctgThetaErrorHist = (TProfile2D *)fin->Get(fCtgThetaResolutionHist.c_str()); 205 ctgThetaErrorHist->SetDirectory(0); 206 fin->Close(); 207 } 208 if(!fUsePhiFormula) 209 { 210 TFile *fin = TFile::Open(fPhiResolutionFile.c_str()); 211 phiErrorHist = (TProfile2D *)fin->Get(fPhiResolutionHist.c_str()); 212 phiErrorHist->SetDirectory(0); 213 fin->Close(); 214 214 } 215 215 216 216 fItInputArray->Reset(); 217 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))218 { 219 217 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) 218 { 219 220 220 const TLorentzVector &momentum = candidate->Momentum; 221 221 const TLorentzVector &position = candidate->InitialPosition; 222 222 223 223 pt = momentum.Pt(); 224 224 eta = momentum.Eta(); … … 226 226 d0 = trueD0 = candidate->D0; 227 227 dz = trueDZ = candidate->DZ; 228 228 229 229 p = trueP = candidate->P; 230 230 ctgTheta = trueCtgTheta = candidate->CtgTheta; 231 231 phi = truePhi = candidate->Phi; 232 232 233 if 233 if(fUseD0Formula) 234 234 d0Error = fD0Formula->Eval(pt, eta); 235 235 else 236 236 { 237 Int_t xbin, ybin; 238 239 xbin = pt < d0ErrorHist->GetXaxis ()->GetXmax () ? d0ErrorHist->GetXaxis ()->FindBin (pt) 240 : d0ErrorHist->GetXaxis ()->GetBinCenter (d0ErrorHist->GetXaxis ()->GetNbins ()); 241 ybin = d0ErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta)); 242 d0Error = d0ErrorHist->GetBinContent (xbin, ybin); 243 if (!d0Error) 244 d0Error = -1.0; 245 } 246 if (d0Error < 0.0) 237 Int_t xbin, ybin; 238 239 xbin = pt < d0ErrorHist->GetXaxis()->GetXmax() ? d0ErrorHist->GetXaxis()->FindBin(pt) : d0ErrorHist->GetXaxis()->GetBinCenter(d0ErrorHist->GetXaxis()->GetNbins()); 240 ybin = d0ErrorHist->GetYaxis()->FindBin(TMath::Abs(eta)); 241 d0Error = d0ErrorHist->GetBinContent(xbin, ybin); 242 if(!d0Error) 243 d0Error = -1.0; 244 } 245 if(d0Error < 0.0) 247 246 continue; 248 247 249 if 248 if(fUseDZFormula) 250 249 dzError = fDZFormula->Eval(pt, eta); 251 250 else 252 251 { 253 Int_t xbin, ybin; 254 255 xbin = pt < dzErrorHist->GetXaxis ()->GetXmax () ? dzErrorHist->GetXaxis ()->FindBin (pt) 256 : dzErrorHist->GetXaxis ()->GetBinCenter (dzErrorHist->GetXaxis ()->GetNbins ()); 257 ybin = dzErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta)); 258 dzError = dzErrorHist->GetBinContent (xbin, ybin); 259 if (!dzError) 260 dzError = -1.0; 261 } 262 if (dzError < 0.0) 252 Int_t xbin, ybin; 253 254 xbin = pt < dzErrorHist->GetXaxis()->GetXmax() ? dzErrorHist->GetXaxis()->FindBin(pt) : dzErrorHist->GetXaxis()->GetBinCenter(dzErrorHist->GetXaxis()->GetNbins()); 255 ybin = dzErrorHist->GetYaxis()->FindBin(TMath::Abs(eta)); 256 dzError = dzErrorHist->GetBinContent(xbin, ybin); 257 if(!dzError) 258 dzError = -1.0; 259 } 260 if(dzError < 0.0) 263 261 continue; 264 262 265 if 263 if(fUsePFormula) 266 264 pError = fPFormula->Eval(pt, eta) * p; 267 265 else 268 266 { 269 Int_t xbin, ybin; 270 271 xbin = pt < pErrorHist->GetXaxis ()->GetXmax () ? pErrorHist->GetXaxis ()->FindBin (pt) 272 : pErrorHist->GetXaxis ()->GetBinCenter (pErrorHist->GetXaxis ()->GetNbins ()); 273 ybin = pErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta)); 274 pError = pErrorHist->GetBinContent (xbin, ybin) * p; 275 if (!pError) 276 pError = -1.0; 277 } 278 if (pError < 0.0) 267 Int_t xbin, ybin; 268 269 xbin = pt < pErrorHist->GetXaxis()->GetXmax() ? pErrorHist->GetXaxis()->FindBin(pt) : pErrorHist->GetXaxis()->GetBinCenter(pErrorHist->GetXaxis()->GetNbins()); 270 ybin = pErrorHist->GetYaxis()->FindBin(TMath::Abs(eta)); 271 pError = pErrorHist->GetBinContent(xbin, ybin) * p; 272 if(!pError) 273 pError = -1.0; 274 } 275 if(pError < 0.0) 279 276 continue; 280 277 281 if 278 if(fUseCtgThetaFormula) 282 279 ctgThetaError = fCtgThetaFormula->Eval(pt, eta); 283 280 else 284 281 { 285 Int_t xbin, ybin; 286 287 xbin = pt < ctgThetaErrorHist->GetXaxis ()->GetXmax () ? ctgThetaErrorHist->GetXaxis ()->FindBin (pt) 288 : ctgThetaErrorHist->GetXaxis ()->GetBinCenter (ctgThetaErrorHist->GetXaxis ()->GetNbins ()); 289 ybin = ctgThetaErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta)); 290 ctgThetaError = ctgThetaErrorHist->GetBinContent (xbin, ybin); 291 if (!ctgThetaError) 292 ctgThetaError = -1.0; 293 } 294 if (ctgThetaError < 0.0) 282 Int_t xbin, ybin; 283 284 xbin = pt < ctgThetaErrorHist->GetXaxis()->GetXmax() ? ctgThetaErrorHist->GetXaxis()->FindBin(pt) : ctgThetaErrorHist->GetXaxis()->GetBinCenter(ctgThetaErrorHist->GetXaxis()->GetNbins()); 285 ybin = ctgThetaErrorHist->GetYaxis()->FindBin(TMath::Abs(eta)); 286 ctgThetaError = ctgThetaErrorHist->GetBinContent(xbin, ybin); 287 if(!ctgThetaError) 288 ctgThetaError = -1.0; 289 } 290 if(ctgThetaError < 0.0) 295 291 continue; 296 292 297 if 293 if(fUsePhiFormula) 298 294 phiError = fPhiFormula->Eval(pt, eta); 299 295 else 300 296 { 301 Int_t xbin, ybin; 302 303 xbin = pt < phiErrorHist->GetXaxis ()->GetXmax () ? phiErrorHist->GetXaxis ()->FindBin (pt) 304 : phiErrorHist->GetXaxis ()->GetBinCenter (phiErrorHist->GetXaxis ()->GetNbins ()); 305 ybin = phiErrorHist->GetYaxis ()->FindBin (TMath::Abs (eta)); 306 phiError = phiErrorHist->GetBinContent (xbin, ybin); 307 if (!phiError) 308 phiError = -1.0; 309 } 310 if (phiError < 0.0) 297 Int_t xbin, ybin; 298 299 xbin = pt < phiErrorHist->GetXaxis()->GetXmax() ? phiErrorHist->GetXaxis()->FindBin(pt) : phiErrorHist->GetXaxis()->GetBinCenter(phiErrorHist->GetXaxis()->GetNbins()); 300 ybin = phiErrorHist->GetYaxis()->FindBin(TMath::Abs(eta)); 301 phiError = phiErrorHist->GetBinContent(xbin, ybin); 302 if(!phiError) 303 phiError = -1.0; 304 } 305 if(phiError < 0.0) 311 306 continue; 312 307 313 if 314 { 315 316 317 318 319 308 if(fApplyToPileUp || !candidate->IsPU) 309 { 310 d0 = gRandom->Gaus(d0, d0Error); 311 dz = gRandom->Gaus(dz, dzError); 312 p = gRandom->Gaus(p, pError); 313 ctgTheta = gRandom->Gaus(ctgTheta, ctgThetaError); 314 phi = gRandom->Gaus(phi, phiError); 320 315 } 321 316 322 317 if(p < 0.0) continue; 323 while (phi > TMath::Pi ()) phi -= TMath::TwoPi();324 while (phi <= -TMath::Pi ()) phi += TMath::TwoPi();318 while(phi > TMath::Pi()) phi -= TMath::TwoPi(); 319 while(phi <= -TMath::Pi()) phi += TMath::TwoPi(); 325 320 326 321 mother = candidate; 327 candidate = static_cast<Candidate *>(candidate->Clone());322 candidate = static_cast<Candidate *>(candidate->Clone()); 328 323 candidate->D0 = d0; 329 324 candidate->DZ = dz; … … 332 327 candidate->Phi = phi; 333 328 334 theta = TMath::ACos(ctgTheta / TMath::Sqrt 335 candidate->Momentum.SetPx (p * TMath::Cos (phi) * TMath::Sin(theta));336 candidate->Momentum.SetPy (p * TMath::Sin (phi) * TMath::Sin(theta));337 candidate->Momentum.SetPz (p * TMath::Cos(theta));338 candidate->Momentum.SetE (candidate->Momentum.Pt () * TMath::CosH(eta));339 candidate->PT = candidate->Momentum.Pt 340 341 x = position.X 342 y = position.Y 343 z = position.Z 344 t = position.T 345 px = candidate->Momentum.Px 346 py = candidate->Momentum.Py 347 pz = candidate->Momentum.Pz 348 pt = candidate->Momentum.Pt 349 329 theta = TMath::ACos(ctgTheta / TMath::Sqrt(1.0 + ctgTheta * ctgTheta)); 330 candidate->Momentum.SetPx(p * TMath::Cos(phi) * TMath::Sin(theta)); 331 candidate->Momentum.SetPy(p * TMath::Sin(phi) * TMath::Sin(theta)); 332 candidate->Momentum.SetPz(p * TMath::Cos(theta)); 333 candidate->Momentum.SetE(candidate->Momentum.Pt() * TMath::CosH(eta)); 334 candidate->PT = candidate->Momentum.Pt(); 335 336 x = position.X(); 337 y = position.Y(); 338 z = position.Z(); 339 t = position.T(); 340 px = candidate->Momentum.Px(); 341 py = candidate->Momentum.Py(); 342 pz = candidate->Momentum.Pz(); 343 pt = candidate->Momentum.Pt(); 344 350 345 // -- solve for delta: d0' = ( (x+delta)*py' - (y+delta)*px' )/pt' 351 352 candidate->InitialPosition.SetX 353 candidate->InitialPosition.SetY 354 x = candidate->InitialPosition.X 355 y = candidate->InitialPosition.Y 356 candidate->InitialPosition.SetZ (z + ((pz * (px * (x - beamSpotPosition.X ()) + py * (y - beamSpotPosition.Y())) + pt * pt * (dz - z)) / (pt * pt)));357 z = candidate->InitialPosition.Z 346 347 candidate->InitialPosition.SetX(x + ((px * y - py * x + d0 * pt) / (py - px))); 348 candidate->InitialPosition.SetY(y + ((px * y - py * x + d0 * pt) / (py - px))); 349 x = candidate->InitialPosition.X(); 350 y = candidate->InitialPosition.Y(); 351 candidate->InitialPosition.SetZ(z + ((pz * (px * (x - beamSpotPosition.X()) + py * (y - beamSpotPosition.Y())) + pt * pt * (dz - z)) / (pt * pt))); 352 z = candidate->InitialPosition.Z(); 358 353 359 354 candidate->InitialPosition.SetT(t); … … 366 361 q = candidate->Charge; 367 362 368 r = pt / (q * fBz) * 1.0E9 /c_light;// in [m]363 r = pt / (q * fBz) * 1.0E9 / c_light; // in [m] 369 364 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi] 370 365 371 366 // 2. helix axis coordinates 372 x_c = x + r *TMath::Sin(phi_0);373 y_c = y - r *TMath::Cos(phi_0);367 x_c = x + r * TMath::Sin(phi_0); 368 y_c = y - r * TMath::Cos(phi_0); 374 369 r_c = TMath::Hypot(x_c, y_c); 375 370 376 371 rcu = TMath::Abs(r); 377 rc2 = r_c *r_c;372 rc2 = r_c * r_c; 378 373 379 374 // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd 380 xd = x_c *x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;375 xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c; 381 376 xd = (rc2 > 0.0) ? xd / rc2 : -999; 382 yd = y_c *(-rcu*r_c + rc2);377 yd = y_c * (-rcu * r_c + rc2); 383 378 yd = (rc2 > 0.0) ? yd / rc2 : -999; 384 zd = z + (TMath::Sqrt(xd *xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;385 386 candidate->Xd = xd *1.0E3;387 candidate->Yd = yd *1.0E3;388 candidate->Zd = zd *1.0E3;389 390 if 391 { 392 393 394 395 396 397 candidate->ErrorPT = ptError(p, ctgTheta, pError, ctgThetaError);398 candidate->TrackResolution = pError/p;399 } 400 379 zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt; 380 381 candidate->Xd = xd * 1.0E3; 382 candidate->Yd = yd * 1.0E3; 383 candidate->Zd = zd * 1.0E3; 384 385 if(fApplyToPileUp || !candidate->IsPU) 386 { 387 candidate->ErrorD0 = d0Error; 388 candidate->ErrorDZ = dzError; 389 candidate->ErrorP = pError; 390 candidate->ErrorCtgTheta = ctgThetaError; 391 candidate->ErrorPhi = phiError; 392 candidate->ErrorPT = ptError(p, ctgTheta, pError, ctgThetaError); 393 candidate->TrackResolution = pError / p; 394 } 395 401 396 candidate->AddCandidate(mother); 402 397 fOutputArray->Add(candidate); … … 406 401 } 407 402 408 Double_t TrackSmearing::ptError 403 Double_t TrackSmearing::ptError(const Double_t p, const Double_t ctgTheta, const Double_t dP, const Double_t dCtgTheta) 409 404 { 410 405 Double_t a, b; 411 406 a = (p * p * ctgTheta * ctgTheta * dCtgTheta * dCtgTheta) / ((ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1) * (ctgTheta * ctgTheta + 1)); 412 407 b = (dP * dP) / (ctgTheta * ctgTheta + 1); 413 return sqrt 414 } 415 416 //------------------------------------------------------------------------------ 408 return sqrt(a + b); 409 } 410 411 //------------------------------------------------------------------------------
Note:
See TracChangeset
for help on using the changeset viewer.