Changeset 4acf2fd in git
- Timestamp:
- May 17, 2021, 6:05:38 PM (3 years ago)
- Branches:
- master
- Children:
- 7dac4ea
- Parents:
- 46d3442 (diff), 4afb18d (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:38)
- git-committer:
- GitHub <noreply@…> (05/17/21 18:05:38)
- Files:
-
- 4 added
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
Makefile
r46d3442 r4acf2fd 412 412 modules/ImpactParameterSmearing.h \ 413 413 modules/TimeSmearing.h \ 414 modules/TimeOfFlight.h \ 414 415 modules/SimpleCalorimeter.h \ 415 416 modules/DenseTrackFilter.h \ … … 448 449 modules/DecayFilter.h \ 449 450 modules/ParticleDensity.h \ 451 modules/TruthVertexFinder.h \ 450 452 modules/ExampleModule.h 451 453 tmp/modules/ModulesDict$(PcmSuf): \ … … 1012 1014 classes/DelphesFactory.h \ 1013 1015 classes/DelphesFormula.h 1016 tmp/modules/TimeOfFlight.$(ObjSuf): \ 1017 modules/TimeOfFlight.$(SrcSuf) \ 1018 modules/TimeOfFlight.h \ 1019 classes/DelphesClasses.h \ 1020 classes/DelphesFactory.h \ 1021 classes/DelphesFormula.h \ 1022 external/ExRootAnalysis/ExRootClassifier.h \ 1023 external/ExRootAnalysis/ExRootFilter.h \ 1024 external/ExRootAnalysis/ExRootResult.h 1014 1025 tmp/modules/TimeSmearing.$(ObjSuf): \ 1015 1026 modules/TimeSmearing.$(SrcSuf) \ … … 1071 1082 external/ExRootAnalysis/ExRootResult.h \ 1072 1083 external/ExRootAnalysis/ExRootTreeBranch.h 1084 tmp/modules/TruthVertexFinder.$(ObjSuf): \ 1085 modules/TruthVertexFinder.$(SrcSuf) \ 1086 modules/TruthVertexFinder.h \ 1087 classes/DelphesClasses.h \ 1088 classes/DelphesFactory.h \ 1089 classes/DelphesPileUpReader.h \ 1090 classes/DelphesTF2.h \ 1091 external/ExRootAnalysis/ExRootClassifier.h \ 1092 external/ExRootAnalysis/ExRootFilter.h \ 1093 external/ExRootAnalysis/ExRootResult.h 1073 1094 tmp/modules/UniqueObjectFinder.$(ObjSuf): \ 1074 1095 modules/UniqueObjectFinder.$(SrcSuf) \ … … 1215 1236 tmp/modules/TaggingParticlesSkimmer.$(ObjSuf) \ 1216 1237 tmp/modules/TauTagging.$(ObjSuf) \ 1238 tmp/modules/TimeOfFlight.$(ObjSuf) \ 1217 1239 tmp/modules/TimeSmearing.$(ObjSuf) \ 1218 1240 tmp/modules/TrackCountingBTagging.$(ObjSuf) \ … … 1222 1244 tmp/modules/TrackSmearing.$(ObjSuf) \ 1223 1245 tmp/modules/TreeWriter.$(ObjSuf) \ 1246 tmp/modules/TruthVertexFinder.$(ObjSuf) \ 1224 1247 tmp/modules/UniqueObjectFinder.$(ObjSuf) \ 1225 1248 tmp/modules/VertexFinder.$(ObjSuf) \ … … 2067 2090 @touch $@ 2068 2091 2092 modules/TimeOfFlight.h: \ 2093 classes/DelphesModule.h 2094 @touch $@ 2095 2069 2096 external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh: \ 2070 2097 external/fastjet/ClusterSequence.hh \ … … 2185 2212 @touch $@ 2186 2213 2214 modules/TruthVertexFinder.h: \ 2215 classes/DelphesModule.h 2216 @touch $@ 2217 2187 2218 classes/DelphesSTDHEPReader.h: \ 2188 2219 classes/DelphesXDRReader.h … … 2300 2331 @touch $@ 2301 2332 2333 modules/TreeWriter.h: \ 2334 classes/DelphesModule.h 2335 @touch $@ 2336 2302 2337 modules/TimeSmearing.h: \ 2303 classes/DelphesModule.h2304 @touch $@2305 2306 modules/TreeWriter.h: \2307 2338 classes/DelphesModule.h 2308 2339 @touch $@ -
cards/delphes_card_IDEA.tcl
r46d3442 r4acf2fd 11 11 set B 2.0 12 12 13 ## Drift chamber coordinates 14 set DCHZMIN -2.125 15 set DCHZMAX 2.125 16 set DCHRMIN 0.345 17 set DCHRMAX 2.02 18 19 13 20 ####################################### 14 21 # Order of execution of various modules … … 16 23 17 24 set ExecutionPath { 25 26 TruthVertexFinder 18 27 ParticlePropagator 19 28 … … 25 34 TrackSmearing 26 35 ClusterCounting 36 TimeSmearing 37 TimeOfFlight 27 38 28 39 TrackMerger … … 61 72 ScalarHT 62 73 TreeWriter 74 } 75 76 ################################# 77 # Truth Vertex Finder 78 ################################# 79 80 module TruthVertexFinder TruthVertexFinder { 81 82 ## below this distance two vertices are assumed to be the same 83 set Resolution 1E-06 84 85 set InputArray Delphes/stableParticles 86 set VertexOutputArray vertices 63 87 } 64 88 … … 101 125 } 102 126 } 103 104 # (pt <= 0.1) * (0.00) +105 # (abs(eta) <= 3.0) * (pt > 0.1) * (1.00) +106 # (abs(eta) > 3) * (0.00)107 127 108 128 … … 212 232 2 VTXDSK 0.141 0.3 0.92 0.00028 0.0937 2 0 1.5708 7e-006 7e-006 1 213 233 214 1 DCHCANI -2.125 2.125 0.3450.0002 0.237223 0 0 0 0 0 0234 1 DCHCANI $DCHZMIN $DCHZMAX $DCHRMIN 0.0002 0.237223 0 0 0 0 0 0 215 235 1 DCH -2 2 0.36 0.0147748 1400 1 0.0203738 0 0.0001 0 1 216 236 1 DCH -2 2 0.374775 0.0147748 1400 1 -0.0212097 0 0.0001 0 1 … … 325 345 1 DCH -2 2 1.98523 0.0147748 1400 1 0.111898 0 0.0001 0 1 326 346 1 DCH -2 2 2 0.0147748 1400 1 -0.112723 0 0.0001 0 1 327 1 DCHCANO -2.125 2.125 2.020.02 1.667 0 0 0 0 0 0347 1 DCHCANO $DCHZMIN $DCHZMAX $DCHRMAX $DCHRMAX 0.02 1.667 0 0 0 0 0 0 328 348 1 BSILWRP -2.35 2.35 2.04 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1 329 349 1 BSILWRP -2.35 2.35 2.06 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1 330 350 1 MAG -2.5 2.5 2.25 0.05 0.0658 0 0 0 0 0 0 331 351 1 BPRESH -2.55 2.55 2.45 0.02 1 2 0 1.5708 7e-005 0.01 1 332 2 DCHWALL 0.345 2.02 2.1250.25 5.55 0 0 0 0 0 0333 2 DCHWALL 0.345 2.02 -2.1250.25 5.55 0 0 0 0 0 0352 2 DCHWALL $DCHRMIN $DCHRMAX $DCHZMAX 0.25 5.55 0 0 0 0 0 0 353 2 DCHWALL $DCHRMIN $DCHRMAX $DCHZMIN 0.25 5.55 0 0 0 0 0 0 334 354 2 FSILWRP 0.354 2.02 -2.32 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1 335 355 2 FSILWRP 0.35 2.02 -2.3 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1 … … 355 375 set Bz $B 356 376 357 set Rmin 0.35 358 set Rmax 2.0 359 set Zmin -2.0 360 set Zmax 2.0 361 362 # gas mix option: 0 377 ## check that these are consistent with DCHCANI/DCHNANO parameters in TrackCovariance module 378 set Rmin $DCHRMIN 379 set Rmax $DCHRMAX 380 set Zmin $DCHZMIN 381 set Zmax $DCHZMAX 382 383 # gas mix option: 363 384 # 0: Helium 90% - Isobutane 10% 364 385 # 1: Helium 100% … … 371 392 372 393 394 ######################################## 395 # Time Smearing MIP 396 ######################################## 397 398 module TimeSmearing TimeSmearing { 399 set TrackInputArray ClusterCounting/tracks 400 set OutputArray tracks 401 402 # assume constant 30 ps resolution for now 403 set TimeResolution { 404 (abs(eta) > 0.0 && abs(eta) <= 3.0)* 30E-12 405 } 406 } 407 408 ######################################## 409 # Time Of Flight Measurement 410 ######################################## 411 412 module TimeOfFlight TimeOfFlight { 413 set TrackInputArray TimeSmearing/tracks 414 set VertexInputArray TruthVertexFinder/vertices 415 416 set OutputArray tracks 417 418 # 0: assume vertex time tV from MC Truth (ideal case) 419 # 1: assume vertex time tV=0 420 # 2: calculate vertex time as vertex TOF, assuming tPV=0 421 set VertexTimeMode 2 422 423 } 424 373 425 374 426 ############## … … 378 430 module Merger TrackMerger { 379 431 # add InputArray InputArray 380 add InputArray ClusterCounting/tracks432 add InputArray TimeOfFlight/tracks 381 433 set OutputArray tracks 382 434 } … … 779 831 780 832 module UniqueObjectFinder UniqueObjectFinder { 833 834 set UseUniqueID true 781 835 # earlier arrays take precedence over later ones 782 836 # add InputArray InputArray OutputArray … … 788 842 789 843 790 791 844 ################## 792 845 # ROOT tree writer -
classes/DelphesClasses.cc
r46d3442 r4acf2fd 231 231 Xd(0), Yd(0), Zd(0), 232 232 Nclusters(0.0), 233 dNdx(0.0), 233 234 TrackResolution(0), 234 235 NCharged(0), … … 408 409 object.Zd = Zd; 409 410 object.Nclusters = Nclusters; 411 object.dNdx = dNdx; 410 412 object.TrackResolution = TrackResolution; 411 413 object.NCharged = NCharged; … … 541 543 Zd = 0.0; 542 544 Nclusters = 0.0; 545 dNdx = 0.0; 543 546 TrackResolution = 0.0; 544 547 NCharged = 0; -
classes/DelphesClasses.h
r46d3442 r4acf2fd 462 462 Float_t DZ; // track longitudinal impact parameter 463 463 Float_t Nclusters; // Number of ionization clusters 464 Float_t dNdx; // Number of ionization clusters 464 465 465 466 Float_t ErrorP; // track momentum error … … 567 568 Float_t DZ; // track longitudinal impact parameter 568 569 Float_t Nclusters; // Number of ionization clusters 570 Float_t dNdx; // Number of ionization clusters 569 571 570 572 Float_t ErrorP; // track momentum error … … 697 699 698 700 Float_t Nclusters; // Number of ionization clusters 701 Float_t dNdx; // Number of ionization clusters per unit length 699 702 700 703 Float_t Xd; -
examples/Example5.C
r46d3442 r4acf2fd 37 37 // Load selected branches with data from specified event 38 38 treeReader->ReadEntry(entry); 39 39 40 40 if(entry>10) break; 41 41 42 42 cout<<"" <<endl; 43 43 cout<<"--------- New Event ---------" <<endl; 44 44 cout<<"" <<endl; 45 45 46 46 // loop over all input particles in the event 47 47 for(Int_t i=0; i < branchParticle->GetEntriesFast(); i++) 48 { 49 GenParticle *gen = (GenParticle*) branchParticle->At(i); 50 cout<<"N: "<<i<<", St: "<<gen->Status<<", PID: "<<gen->PID<<", E: "<<gen->E<<", Px: "<<gen->Px<<", Py: "<<gen->Py<<", Pz: "<<gen->Pz<<", M: "<<gen->Mass<<", M1: "<<gen->M1<<", M2: "<<gen->M2<<", D1: "<<gen->D1<<", D2: "<<gen->D2<<endl;48 { 49 GenParticle *gen = (GenParticle*) branchParticle->At(i); 50 cout<<"N: "<<i<<", St: "<<gen->Status<<", PID: "<<gen->PID<<", E: "<<gen->E<<", Px: "<<gen->Px<<", Py: "<<gen->Py<<", Pz: "<<gen->Pz<<", M: "<<gen->Mass<<", T: "<<gen->T<<", X: "<<gen->X<<", Y: "<<gen->Y<<", Z: "<<gen->Z<<", M1: "<<gen->M1<<", M2: "<<gen->M2<<", D1: "<<gen->D1<<", D2: "<<gen->D2<<endl; 51 51 } 52 52 } -
examples/Pythia8/ee_zh.cmnd
r46d3442 r4acf2fd 1 1 2 Main:numberOfEvents = 10 ! number of events to generate 2 ! number of events to generate 3 Main:numberOfEvents = 1000 ! number of events to generate 3 4 4 5 Beams:idA = 11 ! first beam, e- = -11 … … 7 8 8 9 ! Vertex smearing : 9 Beams:allowVertexSpread = off10 !Beams:sigmaVertexX = 9.70e-3 ! 13.7 mum / sqrt211 !Beams:sigmaVertexY = 25.5E-6 ! 36.1 nm / sqrt212 !Beams:sigmaVertexZ = 20.64 ! 0.64 mm13 !ï¿ŒBeams:sigmaVertexT = 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 -
modules/ClusterCounting.cc
r46d3442 r4acf2fd 1 /*1 /* 2 2 * Delphes: a framework for fast simulation of a generic collider experiment 3 3 * Copyright (C) 2020 Universite catholique de Louvain (UCLouvain), Belgium … … 28 28 */ 29 29 30 //FIXME add reference to Bedeschi-code31 //FIXME make sure about units of P, X32 //FIXME fix pt > 200 GeV issue and angle > 6.4133 34 30 #include "modules/ClusterCounting.h" 35 36 31 #include "classes/DelphesClasses.h" 37 32 #include "TrackCovariance/TrkUtil.h" … … 108 103 { 109 104 Candidate *candidate, *mother, *particle; 110 Double_t mass, Ncl;105 Double_t mass, trackLength, Ncl; 111 106 112 107 fItInputArray->Reset(); … … 124 119 mass = candidateMomentum.M(); 125 120 121 trackLength = fTrackUtil->TrkLen(Par); 122 126 123 mother = candidate; 127 124 candidate = static_cast<Candidate*>(candidate->Clone()); 128 125 129 Ncl = -999; 130 // computation of Nclusters is not supported for electrons 131 if (TMath::Abs(particle->PID) == 11) 126 Ncl = 0.; 127 if (fTrackUtil->IonClusters(Ncl, mass, Par)) 132 128 { 133 129 candidate->Nclusters = Ncl; 130 candidate->dNdx = (trackLength > 0.) ? Ncl/trackLength : -1; 134 131 } 135 else if (fTrackUtil->IonClusters(Ncl, mass, Par))136 {137 candidate->Nclusters = Ncl;138 }139 //cout<<candidate->PID<<", "<<mass<<", "<<candidate->Nclusters<<endl;140 132 141 133 candidate->AddCandidate(mother); -
modules/DualReadoutCalorimeter.cc
r46d3442 r4acf2fd 613 613 // create new photon tower 614 614 tower = static_cast<Candidate*>(fTower->Clone()); 615 pt = 615 pt = neutralEnergy / TMath::CosH(eta); 616 616 //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl; 617 617 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy); … … 623 623 tower->Ehad = 0.0; 624 624 tower->PID = 22; 625 }626 625 fEFlowPhotonOutputArray->Add(tower); 626 } 627 627 // if hadronic fraction > 0, use HCAL resolution 628 628 else … … 631 631 tower->Ehad = neutralEnergy; 632 632 tower->PID = 130; 633 } 634 635 fEFlowPhotonOutputArray->Add(tower); 636 633 fEFlowNeutralHadronOutputArray->Add(tower); 634 } 637 635 638 636 //clone tracks -
modules/Efficiency.cc
r46d3442 r4acf2fd 100 100 Double_t pt, eta, phi, e; 101 101 102 102 103 fItInputArray->Reset(); 103 104 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) … … 115 116 pt = candidateMomentum.Pt(); 116 117 e = candidateMomentum.E(); 117 118 118 119 // apply an efficency formula 119 120 if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e, candidate)) continue; -
modules/ModulesLinkDef.h
r46d3442 r4acf2fd 40 40 #include "modules/ImpactParameterSmearing.h" 41 41 #include "modules/TimeSmearing.h" 42 #include "modules/TimeOfFlight.h" 42 43 #include "modules/SimpleCalorimeter.h" 43 44 #include "modules/DenseTrackFilter.h" … … 76 77 #include "modules/DecayFilter.h" 77 78 #include "modules/ParticleDensity.h" 79 #include "modules/TruthVertexFinder.h" 78 80 #include "modules/ExampleModule.h" 79 81 … … 98 100 #pragma link C++ class ImpactParameterSmearing+; 99 101 #pragma link C++ class TimeSmearing+; 102 #pragma link C++ class TimeOfFlight+; 100 103 #pragma link C++ class SimpleCalorimeter+; 101 104 #pragma link C++ class DenseTrackFilter+; … … 134 137 #pragma link C++ class DecayFilter+; 135 138 #pragma link C++ class ParticleDensity+; 139 #pragma link C++ class TruthVertexFinder+; 136 140 #pragma link C++ class ExampleModule+; 137 141 -
modules/ParticlePropagator.cc
r46d3442 r4acf2fd 125 125 TLorentzVector particlePosition, particleMomentum, beamSpotPosition; 126 126 Double_t px, py, pz, pt, pt2, e, q; 127 Double_t x, y, z, t, r , phi;127 Double_t x, y, z, t, r; 128 128 Double_t x_c, y_c, r_c, phi_c, phi_0; 129 129 Double_t x_t, y_t, z_t, r_t; 130 Double_t t1, t2, t3, t4, t5, t6; 131 Double_t t_z, t_r, t_ra, t_rb; 132 Double_t tmp, discr, discr2; 133 Double_t delta, gammam, omega, asinrho; 134 Double_t rcu, rc2, xd, yd, zd; 135 Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha; 130 Double_t t_z, t_r; 131 Double_t discr; 132 Double_t gammam, omega; 133 Double_t xd, yd, zd; 134 Double_t l, d0, dz, ctgTheta, alpha; 136 135 Double_t bsx, bsy, bsz; 137 Double_t s0, s1, sd; 138 136 Double_t rxp, rdp, t_R; 137 Double_t td, pio, phid, sign_pz, vz; 138 139 139 const Double_t c_light = 2.99792458E8; 140 140 … … 161 161 particlePosition = particle->Position; 162 162 particleMomentum = particle->Momentum; 163 164 // Constants 165 163 166 x = particlePosition.X() * 1.0E-3; 164 167 y = particlePosition.Y() * 1.0E-3; … … 205 208 else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9) 206 209 { 207 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0 208 tmp = px * y - py * x; 209 discr2 = pt2 * fRadius2 - tmp * tmp; 210 211 if(discr2 < 0.0) 212 { 213 // no solutions 214 continue; 215 } 216 217 tmp = px * x + py * y; 218 discr = TMath::Sqrt(discr2); 219 t1 = (-tmp + discr) / pt2; 220 t2 = (-tmp - discr) / pt2; 221 t = (t1 < 0.0) ? t2 : t1; 222 223 z_t = z + pz * t; 224 if(TMath::Abs(z_t) > fHalfLength) 225 { 226 t3 = (+fHalfLength - z) / pz; 227 t4 = (-fHalfLength - z) / pz; 228 t = (t3 < 0.0) ? t4 : t3; 229 } 230 231 x_t = x + px * t; 232 y_t = y + py * t; 233 z_t = z + pz * t; 234 235 l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z)); 210 211 rxp = x*py - y*px; 212 rdp = x*px + y*py; 213 214 discr = fRadius*fRadius*pt*pt - rxp*rxp; 215 216 t_R = e * (sqrt(discr) - rdp) / (c_light * pt * pt); 217 t_z = e * (TMath::Sign(fHalfLengthMax, pz) - z) / ( c_light * pz); 218 219 t = TMath::Min(t_R, t_z); 220 221 x_t = x + px*t*c_light/e; 222 y_t = y + py*t*c_light/e; 223 z_t = z + pz*t*c_light/e; 224 r_t = TMath::Hypot(x_t, y_t); 225 226 l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z)); 236 227 237 228 mother = candidate; 238 candidate = static_cast<Candidate 229 candidate = static_cast<Candidate*>(candidate->Clone()); 239 230 240 231 candidate->InitialPosition = particlePosition; 241 candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e *1.0E3);242 candidate->L = l *1.0E3;232 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3); 233 candidate->L = l*1.0E3; 243 234 244 235 candidate->Momentum = particleMomentum; … … 246 237 247 238 fOutputArray->Add(candidate); 239 248 240 if(TMath::Abs(q) > 1.0E-9) 249 241 { … … 274 266 // helix radius r = p_{T0} / (omega gamma m) 275 267 276 gammam = e * 1.0E9 / (c_light * c_light);// gammam in [eV/c^2]277 omega = q * fBz / (gammam); // omega is here in [89875518/s]278 r = pt / (q * fBz) * 1.0E9 / c_light;// in [m]268 gammam = e*1.0E9 / (c_light*c_light); // gammam in [eV/c^2] 269 omega = q * fBz / (gammam); // omega is here in [89875518/s] 270 r = pt / (q * fBz) * 1.0E9/c_light; // in [m] 279 271 280 272 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi] 281 273 282 274 // 2. helix axis coordinates 283 x_c = x + r *TMath::Sin(phi_0);284 y_c = y - r *TMath::Cos(phi_0);275 x_c = x + r*TMath::Sin(phi_0); 276 y_c = y - r*TMath::Cos(phi_0); 285 277 r_c = TMath::Hypot(x_c, y_c); 286 phi_c = TMath::ATan2(y_c, x_c); 287 phi = phi_c; 288 if(x_c < 0.0) phi += TMath::Pi(); 289 290 rcu = TMath::Abs(r); 291 rc2 = r_c * r_c; 292 293 // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd 294 xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c; 295 xd = (rc2 > 0.0) ? xd / rc2 : -999; 296 yd = y_c * (-rcu * r_c + rc2); 297 yd = (rc2 > 0.0) ? yd / rc2 : -999; 298 zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt; 299 300 // proper calculation of the DCAz coordinate 301 // s0: track circle parameter at the track origin 302 // s1: track circle parameter at the closest approach to beam pipe 303 // sd: s1-s0 signed angular difference 304 s0 = atan2(y - y_c, x - x_c); 305 s1 = atan2(yd - y_c, xd - x_c); 306 sd = atan2(sin(s1 - s0), cos(s1 - s0)); 307 zd = z - r * pz / pt * sd; 308 309 // use perigee momentum rather than original particle 310 // momentum, since the orignal particle momentum isn't known 311 312 px = TMath::Sign(1.0, r) * pt * (-y_c / r_c); 313 py = TMath::Sign(1.0, r) * pt * (x_c / r_c); 314 etap = particleMomentum.Eta(); 315 phip = TMath::ATan2(py, px); 316 317 particleMomentum.SetPtEtaPhiE(pt, etap, phip, particleMomentum.E()); 278 phi_c = TMath::ATan(y_c/x_c); 279 if(x_c < 0.0) phi_c -= TMath::Sign(1., phi_c)*TMath::Pi(); 280 281 //Find the time of closest approach 282 td = (phi_0 - TMath::ATan(-x_c/y_c))/omega; 283 284 //Remove all the modulo pi that might have come from the atan 285 pio = fabs(TMath::Pi()/omega); 286 while(fabs(td) > 0.5*pio) 287 { 288 td -= TMath::Sign(1., td)*pio; 289 } 290 291 //Compute the coordinate of closed approach to z axis 292 //if wants wtr beamline need to be changedto re-center with a traslation of the z axis 293 phid = phi_0 - omega*td; 294 xd = x_c - r*TMath::Sin(phid); 295 yd = y_c + r*TMath::Cos(phid); 296 zd = z + c_light*(pz/e)*td; 297 298 //Compute momentum at closest approach (perigee??) 299 px = pt*TMath::Cos(phid); 300 py = pt*TMath::Sin(phid); 301 302 particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E()); 318 303 319 304 // calculate additional track parameters (correct for beamspot position) 320 321 d0 = ((x - bsx) * py - (y - bsy) * px) / pt; 322 dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt); 323 p = particleMomentum.P(); 324 ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta()); 305 d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt; 306 dz = zd - bsz; 307 ctgTheta = 1.0 / TMath::Tan (particleMomentum.Theta()); 325 308 326 309 // 3. time evaluation t = TMath::Min(t_r, t_z) 327 310 // t_r : time to exit from the sides 328 311 // t_z : time to exit from the front or the back 329 t_r = 0.0; // in [ns] 330 int sign_pz = (pz > 0.0) ? 1 : -1; 331 if(pz == 0.0) 332 t_z = 1.0E99; 312 t = 0; 313 t_z = 0; 314 sign_pz = (pz > 0.0) ? 1 : -1; 315 if(pz == 0.0) t_z = 1.0E99; 316 else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz); 317 318 if(r_c + TMath::Abs(r) < fRadius) // helix does not cross the cylinder sides 319 { 320 t = t_z; 321 } 333 322 else 334 t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz); 335 336 if(r_c + TMath::Abs(r) < fRadius) 337 { 338 // helix does not cross the cylinder sides 339 t = t_z; 340 } 341 else 342 { 343 asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c)); 344 delta = phi_0 - phi; 345 if(delta < -TMath::Pi()) delta += 2 * TMath::Pi(); 346 if(delta > TMath::Pi()) delta -= 2 * TMath::Pi(); 347 t1 = (delta + asinrho) / omega; 348 t2 = (delta + TMath::Pi() - asinrho) / omega; 349 t3 = (delta + TMath::Pi() + asinrho) / omega; 350 t4 = (delta - asinrho) / omega; 351 t5 = (delta - TMath::Pi() - asinrho) / omega; 352 t6 = (delta - TMath::Pi() + asinrho) / omega; 353 354 if(t1 < 0.0) t1 = 1.0E99; 355 if(t2 < 0.0) t2 = 1.0E99; 356 if(t3 < 0.0) t3 = 1.0E99; 357 if(t4 < 0.0) t4 = 1.0E99; 358 if(t5 < 0.0) t5 = 1.0E99; 359 if(t6 < 0.0) t6 = 1.0E99; 360 361 t_ra = TMath::Min(t1, TMath::Min(t2, t3)); 362 t_rb = TMath::Min(t4, TMath::Min(t5, t6)); 363 t_r = TMath::Min(t_ra, t_rb); 323 { 324 alpha = -(fRadius*fRadius - r*r - r_c*r_c)/(2*fabs(r)*r_c); 325 alpha = fabs(TMath::ACos(alpha)); 326 t_r = td + alpha/fabs(omega); 327 364 328 t = TMath::Min(t_r, t_z); 365 329 } 366 330 367 // 4. position in terms of x(t), y(t), z(t) 368 x_t = x_c + r * TMath::Sin(omega * t - phi_0); 369 y_t = y_c + r * TMath::Cos(omega * t - phi_0); 370 z_t = z + pz * 1.0E9 / c_light / gammam * t; 371 r_t = TMath::Hypot(x_t, y_t); 331 x_t = x_c - r*TMath::Sin(phi_0 - omega*t); 332 y_t = y_c + r*TMath::Cos(phi_0 - omega*t); 333 z_t = z + c_light*t*pz/e; 334 r_t = TMath::Hypot(x_t, y_t); 372 335 373 336 // compute path length for an helix 374 375 alpha = pz * 1.0E9 / c_light / gammam;376 l = t * TMath::Sqrt( alpha * alpha + r * r * omega *omega);337 vz = pz*1.0E9 / c_light / gammam; 338 //lenght of the path from production to tracker 339 l = t * TMath::Sqrt(vz*vz + r*r*omega*omega); 377 340 378 341 if(r_t > 0.0) 379 342 { 380 381 343 // store these variables before cloning 382 344 if(particle == candidate) … … 384 346 particle->D0 = d0 * 1.0E3; 385 347 particle->DZ = dz * 1.0E3; 386 particle->P = p ;348 particle->P = particleMomentum.P(); 387 349 particle->PT = pt; 388 350 particle->CtgTheta = ctgTheta; 389 particle->Phi = p hip;351 particle->Phi = particleMomentum.Phi(); 390 352 } 391 353 -
modules/TimeSmearing.cc
r46d3442 r4acf2fd 19 19 /** \class TimeSmearing 20 20 * 21 * Performs t ransverse momentum resolutionsmearing.21 * Performs time smearing. 22 22 * 23 * \author P. Demin - UCL, Louvain-la-Neuve23 * \author M. Selvaggi - CERN 24 24 * 25 25 */ … … 49 49 50 50 using namespace std; 51 52 51 //------------------------------------------------------------------------------ 53 52 54 53 TimeSmearing::TimeSmearing() : 55 fIt InputArray(0)54 fItTrackInputArray(0), fResolutionFormula(0) 56 55 { 56 fResolutionFormula = new DelphesFormula; 57 57 } 58 58 … … 61 61 TimeSmearing::~TimeSmearing() 62 62 { 63 if(fResolutionFormula) delete fResolutionFormula; 63 64 } 64 65 … … 69 70 // read resolution formula 70 71 71 fTimeResolution = GetDouble("TimeResolution", 1.0E-10);72 // import input array72 // read time resolution formula in seconds 73 fResolutionFormula->Compile(GetString("TimeResolution", "30e-12")); 73 74 74 fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons")); 75 fItInputArray = fInputArray->MakeIterator(); 75 // import track input array 76 fTrackInputArray = ImportArray(GetString("TrackInputArray", "MuonMomentumSmearing/muons")); 77 fItTrackInputArray = fTrackInputArray->MakeIterator(); 78 76 79 77 80 // create output array 78 79 fOutputArray = ExportArray(GetString("OutputArray", "muons")); 81 fOutputArray = ExportArray(GetString("OutputArray", "tracks")); 80 82 } 81 83 … … 84 86 void TimeSmearing::Finish() 85 87 { 86 if(fIt InputArray) delete fItInputArray;88 if(fItTrackInputArray) delete fItTrackInputArray; 87 89 } 88 90 … … 92 94 { 93 95 Candidate *candidate, *mother; 94 Double_t ti, tf_smeared, tf; 96 Double_t tf_smeared, tf; 97 Double_t eta, energy; 98 Double_t timeResolution; 99 95 100 const Double_t c_light = 2.99792458E8; 96 101 97 fIt InputArray->Reset();98 while((candidate = static_cast<Candidate *>(fIt InputArray->Next())))102 fItTrackInputArray->Reset(); 103 while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next()))) 99 104 { 100 const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;105 // converting to meters 101 106 const TLorentzVector &candidateFinalPosition = candidate->Position; 107 const TLorentzVector &candidateMomentum = candidate->Momentum; 102 108 103 ti = candidateInitialPosition.T() * 1.0E-3 / c_light;104 109 tf = candidateFinalPosition.T() * 1.0E-3 / c_light; 105 110 111 eta = candidateMomentum.Eta(); 112 energy = candidateMomentum.E(); 113 106 114 // apply smearing formula 107 tf_smeared = gRandom->Gaus(tf, fTimeResolution); 108 ti = ti + tf_smeared - tf; 109 tf = tf_smeared; 115 timeResolution = fResolutionFormula->Eval(0.0, eta, 0.0, energy); 116 tf_smeared = gRandom->Gaus(tf, timeResolution); 110 117 111 118 mother = candidate; 112 119 candidate = static_cast<Candidate *>(candidate->Clone()); 113 candidate->InitialPosition.SetT(ti * 1.0E3 * c_light);114 candidate->Position.SetT(tf * 1.0E3 * c_light);115 120 116 candidate->ErrorT = fTimeResolution * 1.0E3 * c_light; 121 candidate->Position.SetT(tf_smeared * 1.0E3 * c_light); 122 candidate->ErrorT = timeResolution * 1.0E3 * c_light; 117 123 118 124 candidate->AddCandidate(mother); 119 120 125 fOutputArray->Add(candidate); 121 126 } 122 127 } 123 124 //------------------------------------------------------------------------------ -
modules/TimeSmearing.h
r46d3442 r4acf2fd 22 22 /** \class TimeSmearing 23 23 * 24 * Performs t ransverse time smearing.24 * Performs time smearing. 25 25 * 26 * \author Michele Selvaggi - UCL, Louvain-la-Neuve26 * \author Michele Selvaggi - CERN 27 27 * 28 28 */ … … 32 32 class TIterator; 33 33 class TObjArray; 34 class DelphesFormula; 34 35 35 36 class TimeSmearing: public DelphesModule … … 44 45 45 46 private: 46 Double_t fTimeResolution;47 47 48 TIterator *fItInputArray; //! 48 DelphesFormula *fResolutionFormula; 49 Int_t fVertexTimeMode; 49 50 50 const TObjArray *fInputArray; //! 51 TIterator *fItTrackInputArray; //! 52 53 const TObjArray *fTrackInputArray; //! 51 54 52 55 TObjArray *fOutputArray; //! -
modules/TrackCovariance.cc
r46d3442 r4acf2fd 104 104 Double_t dd0, ddz, dphi, dct, dp, dpt, dC; 105 105 106 107 106 fItInputArray->Reset(); 108 107 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) -
modules/TreeWriter.cc
r46d3442 r4acf2fd 356 356 entry->DZ = candidate->DZ; 357 357 entry->Nclusters = candidate->Nclusters; 358 entry->dNdx = candidate->dNdx; 358 359 359 360 entry->ErrorP = candidate->ErrorP; … … 404 405 entry->Mass = m; 405 406 406 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 407 const TLorentzVector &initialPosition = particle->Position; 407 //particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 408 //const TLorentzVector &initialPosition = particle->Position; 409 const TLorentzVector &initialPosition = candidate->InitialPosition; 408 410 409 411 entry->X = initialPosition.X(); … … 411 413 entry->Z = initialPosition.Z(); 412 414 entry->T = initialPosition.T() * 1.0E-3 / c_light; 415 entry->ErrorT =candidate-> ErrorT * 1.0E-3 / c_light; 413 416 414 417 entry->Particle = particle; … … 509 512 entry->DZ = candidate->DZ; 510 513 entry->Nclusters = candidate->Nclusters; 514 entry->dNdx = candidate->dNdx; 511 515 512 516 entry->ErrorP = candidate->ErrorP; … … 557 561 entry->Mass = m; 558 562 559 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 560 const TLorentzVector &initialPosition = particle->Position; 563 //particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 564 //const TLorentzVector &initialPosition = particle->Position; 565 const TLorentzVector &initialPosition = candidate->InitialPosition; 561 566 562 567 entry->X = initialPosition.X(); … … 564 569 entry->Z = initialPosition.Z(); 565 570 entry->T = initialPosition.T() * 1.0E-3 / c_light; 571 entry->ErrorT = candidate-> ErrorT * 1.0E-3 / c_light; 566 572 567 573 entry->VertexIndex = candidate->ClusterIndex; … … 574 580 entry->Edges[3] = candidate->Edges[3]; 575 581 576 entry->T = position.T() * 1.0E-3 / c_light;582 //entry->T = position.T() * 1.0E-3 / c_light; 577 583 entry->NTimeHits = candidate->NTimeHits; 578 584
Note:
See TracChangeset
for help on using the changeset viewer.