Changeset 28027d5 in git
- Timestamp:
- Jun 26, 2015, 3:13:55 PM (9 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- f3c4047
- Parents:
- f53a4d2 (diff), fe0273c (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. - Files:
-
- 27 added
- 26 edited
Legend:
- Unmodified
- Added
- Removed
-
Makefile
rf53a4d2 r28027d5 257 257 classes/DelphesClasses.h \ 258 258 classes/DelphesFactory.h \ 259 classes/DelphesLHEFReader.h \ 259 260 external/ExRootAnalysis/ExRootTreeWriter.h \ 260 261 external/ExRootAnalysis/ExRootTreeBranch.h \ … … 337 338 modules/Weighter.h \ 338 339 modules/Hector.h \ 340 modules/RunPUPPI.h \ 341 modules/JetFlavorAssociation.h \ 339 342 modules/ExampleModule.h 340 343 ModulesDict$(PcmSuf): \ … … 521 524 tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \ 522 525 external/Hector/H_VerticalQuadrupole.$(SrcSuf) 526 tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf): \ 527 external/PUPPI/puppiCleanContainer.$(SrcSuf) \ 528 external/fastjet/Selector.hh 523 529 tmp/modules/AngularSmearing.$(ObjSuf): \ 524 530 modules/AngularSmearing.$(SrcSuf) \ … … 652 658 external/ExRootAnalysis/ExRootFilter.h \ 653 659 external/ExRootAnalysis/ExRootClassifier.h 660 tmp/modules/JetFlavorAssociation.$(ObjSuf): \ 661 modules/JetFlavorAssociation.$(SrcSuf) \ 662 modules/JetFlavorAssociation.h \ 663 classes/DelphesClasses.h \ 664 classes/DelphesFactory.h \ 665 classes/DelphesFormula.h \ 666 external/ExRootAnalysis/ExRootResult.h \ 667 external/ExRootAnalysis/ExRootFilter.h \ 668 external/ExRootAnalysis/ExRootClassifier.h 654 669 tmp/modules/JetPileUpSubtractor.$(ObjSuf): \ 655 670 modules/JetPileUpSubtractor.$(SrcSuf) \ … … 730 745 classes/DelphesClasses.h \ 731 746 classes/DelphesFactory.h \ 732 classes/Delphes Formula.h \747 classes/DelphesTF2.h \ 733 748 classes/DelphesPileUpReader.h \ 734 749 external/ExRootAnalysis/ExRootResult.h \ 735 750 external/ExRootAnalysis/ExRootFilter.h \ 736 751 external/ExRootAnalysis/ExRootClassifier.h 752 tmp/modules/RunPUPPI.$(ObjSuf): \ 753 modules/RunPUPPI.$(SrcSuf) \ 754 modules/RunPUPPI.h \ 755 external/PUPPI/puppiCleanContainer.hh \ 756 external/PUPPI/RecoObj.hh \ 757 external/PUPPI/puppiParticle.hh \ 758 external/PUPPI/puppiAlgoBin.hh \ 759 classes/DelphesClasses.h \ 760 classes/DelphesFactory.h \ 761 classes/DelphesFormula.h 737 762 tmp/modules/SimpleCalorimeter.$(ObjSuf): \ 738 763 modules/SimpleCalorimeter.$(SrcSuf) \ … … 869 894 tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \ 870 895 tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \ 896 tmp/external/PUPPI/puppiCleanContainer.$(ObjSuf) \ 871 897 tmp/modules/AngularSmearing.$(ObjSuf) \ 872 898 tmp/modules/BTagging.$(ObjSuf) \ … … 883 909 tmp/modules/ImpactParameterSmearing.$(ObjSuf) \ 884 910 tmp/modules/Isolation.$(ObjSuf) \ 911 tmp/modules/JetFlavorAssociation.$(ObjSuf) \ 885 912 tmp/modules/JetPileUpSubtractor.$(ObjSuf) \ 886 913 tmp/modules/LeptonDressing.$(ObjSuf) \ … … 891 918 tmp/modules/PileUpJetID.$(ObjSuf) \ 892 919 tmp/modules/PileUpMerger.$(ObjSuf) \ 920 tmp/modules/RunPUPPI.$(ObjSuf) \ 893 921 tmp/modules/SimpleCalorimeter.$(ObjSuf) \ 894 922 tmp/modules/StatusPidFilter.$(ObjSuf) \ … … 1080 1108 tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf): \ 1081 1109 external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(SrcSuf) 1110 tmp/external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(ObjSuf): \ 1111 external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(SrcSuf) \ 1112 external/fastjet/JetDefinition.hh \ 1113 external/fastjet/ClusterSequenceAreaBase.hh 1114 tmp/external/fastjet/contribs/RecursiveTools/Recluster.$(ObjSuf): \ 1115 external/fastjet/contribs/RecursiveTools/Recluster.$(SrcSuf) 1116 tmp/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(ObjSuf): \ 1117 external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(SrcSuf) \ 1118 external/fastjet/JetDefinition.hh \ 1119 external/fastjet/ClusterSequenceAreaBase.hh 1120 tmp/external/fastjet/contribs/RecursiveTools/SoftDrop.$(ObjSuf): \ 1121 external/fastjet/contribs/RecursiveTools/SoftDrop.$(SrcSuf) 1082 1122 tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf): \ 1083 1123 external/fastjet/contribs/SoftKiller/SoftKiller.$(SrcSuf) … … 1213 1253 external/fastjet/contribs/Nsubjettiness/Njettiness.hh \ 1214 1254 external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh \ 1215 external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh 1255 external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh \ 1256 external/fastjet/tools/Filter.hh \ 1257 external/fastjet/tools/Pruner.hh \ 1258 external/fastjet/contribs/RecursiveTools/SoftDrop.hh 1216 1259 tmp/modules/FastJetGridMedianEstimator.$(ObjSuf): \ 1217 1260 modules/FastJetGridMedianEstimator.$(SrcSuf) \ … … 1285 1328 tmp/external/fastjet/contribs/Nsubjettiness/Nsubjettiness.$(ObjSuf) \ 1286 1329 tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf) \ 1330 tmp/external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.$(ObjSuf) \ 1331 tmp/external/fastjet/contribs/RecursiveTools/Recluster.$(ObjSuf) \ 1332 tmp/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.$(ObjSuf) \ 1333 tmp/external/fastjet/contribs/RecursiveTools/SoftDrop.$(ObjSuf) \ 1287 1334 tmp/external/fastjet/contribs/SoftKiller/SoftKiller.$(ObjSuf) \ 1288 1335 tmp/external/fastjet/plugins/ATLASCone/ATLASConePlugin.$(ObjSuf) \ … … 1338 1385 display/Delphes3DGeometry.$(SrcSuf) \ 1339 1386 display/Delphes3DGeometry.h \ 1340 external/ExRootAnalysis/ExRootConfReader.h \1341 classes/DelphesClasses.h1387 classes/DelphesClasses.h \ 1388 external/ExRootAnalysis/ExRootConfReader.h 1342 1389 tmp/display/DelphesBranchElement.$(ObjSuf): \ 1343 1390 display/DelphesBranchElement.$(SrcSuf) \ … … 1352 1399 tmp/display/DelphesEventDisplay.$(ObjSuf): \ 1353 1400 display/DelphesEventDisplay.$(SrcSuf) \ 1354 external/ExRootAnalysis/ExRootConfReader.h \1355 external/ExRootAnalysis/ExRootTreeReader.h \1356 1401 display/DelphesCaloData.h \ 1357 1402 display/DelphesBranchElement.h \ 1358 1403 display/Delphes3DGeometry.h \ 1359 1404 display/DelphesEventDisplay.h \ 1360 classes/DelphesClasses.h 1405 display/DelphesDisplay.h \ 1406 display/Delphes3DGeometry.h \ 1407 display/DelphesHtmlSummary.h \ 1408 display/DelphesPlotSummary.h \ 1409 classes/DelphesClasses.h \ 1410 external/ExRootAnalysis/ExRootConfReader.h \ 1411 external/ExRootAnalysis/ExRootTreeReader.h 1361 1412 tmp/display/DelphesHtmlSummary.$(ObjSuf): \ 1362 1413 display/DelphesHtmlSummary.$(SrcSuf) \ … … 1541 1592 @touch $@ 1542 1593 1543 external/fastjet/Selector.hh: \1544 external/fastjet/PseudoJet.hh \1545 external/fastjet/RangeDefinition.hh1546 @touch $@1547 1548 1594 external/fastjet/internal/Dnn2piCylinder.hh: \ 1549 1595 external/fastjet/internal/DynamicNearestNeighbours.hh \ 1550 1596 external/fastjet/internal/DnnPlane.hh \ 1551 1597 external/fastjet/internal/numconsts.hh 1598 @touch $@ 1599 1600 external/fastjet/Selector.hh: \ 1601 external/fastjet/PseudoJet.hh \ 1602 external/fastjet/RangeDefinition.hh 1552 1603 @touch $@ 1553 1604 … … 1639 1690 @touch $@ 1640 1691 1692 modules/RunPUPPI.h: \ 1693 classes/DelphesModule.h 1694 @touch $@ 1695 1641 1696 modules/Cloner.h: \ 1642 1697 classes/DelphesModule.h … … 1675 1730 @touch $@ 1676 1731 1677 display/DelphesEventDisplay.h: \1678 external/ExRootAnalysis/ExRootTreeReader.h \1679 display/DelphesDisplay.h \1680 display/Delphes3DGeometry.h \1681 display/DelphesHtmlSummary.h \1682 display/DelphesPlotSummary.h1683 @touch $@1684 1685 1732 modules/TauTagging.h: \ 1686 1733 classes/DelphesModule.h \ … … 1725 1772 @touch $@ 1726 1773 1774 modules/JetFlavorAssociation.h: \ 1775 classes/DelphesModule.h \ 1776 classes/DelphesClasses.h 1777 @touch $@ 1778 1727 1779 modules/ParticlePropagator.h: \ 1728 1780 classes/DelphesModule.h … … 1757 1809 external/fastjet/AreaDefinition.hh \ 1758 1810 external/fastjet/ClusterSequenceAreaBase.hh 1811 @touch $@ 1812 1813 external/PUPPI/puppiCleanContainer.hh: \ 1814 external/PUPPI/RecoObj.hh \ 1815 external/PUPPI/puppiParticle.hh \ 1816 external/PUPPI/puppiAlgoBin.hh \ 1817 external/fastjet/internal/base.hh \ 1818 external/fastjet/PseudoJet.hh 1759 1819 @touch $@ 1760 1820 -
classes/ClassesLinkDef.h
rf53a4d2 r28027d5 46 46 #pragma link C++ class LHCOEvent+; 47 47 #pragma link C++ class LHEFEvent+; 48 #pragma link C++ class LHEFWeight+; 48 49 #pragma link C++ class HepMCEvent+; 49 50 #pragma link C++ class GenParticle+; -
classes/DelphesClasses.cc
rf53a4d2 r28027d5 120 120 PID(0), Status(0), M1(-1), M2(-1), D1(-1), D2(-1), 121 121 Charge(0), Mass(0.0), 122 IsPU(0), IsConstituent(0), 123 BTag(0), TauTag(0), Eem(0.0), Ehad(0.0), 122 IsPU(0), IsRecoPU(0), IsConstituent(0), 123 Flavor(0), FlavorAlgo(0), FlavorPhys(0), 124 BTag(0), BTagAlgo(0), BTagPhys(0), 125 TauTag(0), Eem(0.0), Ehad(0.0), 124 126 DeltaEta(0.0), DeltaPhi(0.0), 125 127 Momentum(0.0, 0.0, 0.0, 0.0), … … 129 131 NCharged(0), 130 132 NNeutrals(0), 133 NSubJetsTrimmed(0), 134 NSubJetsPruned(0), 135 NSubJetsSoftDropped(0), 131 136 Beta(0), 132 137 BetaStar(0), 133 138 MeanSqDeltaR(0), 134 139 PTD(0), 140 Ntimes(-1), 141 IsolationVar(-999), 142 IsolationVarRhoCorr(-999), 143 SumPtCharged(-999), 144 SumPtNeutral(-999), 145 SumPtChargedPU(-999), 146 SumPt(-999), 135 147 fFactory(0), 136 148 fArray(0) 137 149 { 150 int i; 138 151 Edges[0] = 0.0; 139 152 Edges[1] = 0.0; … … 150 163 Tau[3] = 0.0; 151 164 Tau[4] = 0.0; 165 for(i = 0; i < 5; ++i) 166 { 167 TrimmedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0); 168 PrunedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0); 169 SoftDroppedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0); 170 } 152 171 } 153 172 … … 224 243 object.IsPU = IsPU; 225 244 object.IsConstituent = IsConstituent; 245 object.Flavor = Flavor; 246 object.FlavorAlgo = FlavorAlgo; 247 object.FlavorPhys = FlavorPhys; 226 248 object.BTag = BTag; 249 object.BTagAlgo = BTagAlgo; 250 object.BTagPhys = BTagPhys; 227 251 object.TauTag = TauTag; 228 252 object.Eem = Eem; … … 249 273 object.MeanSqDeltaR = MeanSqDeltaR; 250 274 object.PTD = PTD; 275 object.Ntimes = Ntimes; 276 object.IsolationVar = IsolationVar; 277 object.IsolationVarRhoCorr = IsolationVarRhoCorr; 278 object.SumPtCharged = SumPtCharged; 279 object.SumPtNeutral = SumPtNeutral; 280 object.SumPtChargedPU = SumPtChargedPU; 281 object.SumPt = SumPt; 282 251 283 object.FracPt[0] = FracPt[0]; 252 284 object.FracPt[1] = FracPt[1]; … … 259 291 object.Tau[3] = Tau[3]; 260 292 object.Tau[4] = Tau[4]; 293 294 object.TrimmedP4[0] = TrimmedP4[0]; 295 object.TrimmedP4[1] = TrimmedP4[1]; 296 object.TrimmedP4[2] = TrimmedP4[2]; 297 object.TrimmedP4[3] = TrimmedP4[3]; 298 object.TrimmedP4[4] = TrimmedP4[4]; 299 object.PrunedP4[0] = PrunedP4[0]; 300 object.PrunedP4[1] = PrunedP4[1]; 301 object.PrunedP4[2] = PrunedP4[2]; 302 object.PrunedP4[3] = PrunedP4[3]; 303 object.PrunedP4[4] = PrunedP4[4]; 304 object.SoftDroppedP4[0] = SoftDroppedP4[0]; 305 object.SoftDroppedP4[1] = SoftDroppedP4[1]; 306 object.SoftDroppedP4[2] = SoftDroppedP4[2]; 307 object.SoftDroppedP4[3] = SoftDroppedP4[3]; 308 object.SoftDroppedP4[4] = SoftDroppedP4[4]; 309 310 object.NSubJetsTrimmed = NSubJetsTrimmed; 311 object.NSubJetsPruned = NSubJetsPruned; 312 object.NSubJetsSoftDropped = NSubJetsSoftDropped; 261 313 262 314 object.fFactory = fFactory; 263 315 object.fArray = 0; 316 317 // Copy cluster timing info 318 copy(Ecal_E_t.begin(),Ecal_E_t.end(),back_inserter(object.Ecal_E_t)); 264 319 265 320 if(fArray && fArray->GetEntriesFast() > 0) … … 278 333 void Candidate::Clear(Option_t* option) 279 334 { 335 int i; 280 336 SetUniqueID(0); 281 337 ResetBit(kIsReferenced); … … 287 343 IsPU = 0; 288 344 IsConstituent = 0; 345 Flavor = 0; 346 FlavorAlgo = 0; 347 FlavorPhys = 0; 289 348 BTag = 0; 349 BTagAlgo = 0; 350 BTagPhys = 0; 290 351 TauTag = 0; 291 352 Eem = 0.0; … … 311 372 MeanSqDeltaR = 0.0; 312 373 PTD = 0.0; 374 375 Ntimes = 0; 376 Ecal_E_t.clear(); 377 378 IsolationVar = -999; 379 IsolationVarRhoCorr = -999; 380 SumPtCharged = -999; 381 SumPtNeutral = -999; 382 SumPtChargedPU = -999; 383 SumPt = -999; 384 313 385 FracPt[0] = 0.0; 314 386 FracPt[1] = 0.0; … … 321 393 Tau[3] = 0.0; 322 394 Tau[4] = 0.0; 323 395 396 for(i = 0; i < 5; ++i) 397 { 398 TrimmedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0); 399 PrunedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0); 400 SoftDroppedP4[i].SetXYZT(0.0, 0.0, 0.0, 0.0); 401 } 402 403 NSubJetsTrimmed = 0; 404 NSubJetsPruned = 0; 405 NSubJetsSoftDropped = 0; 406 324 407 fArray = 0; 325 408 } -
classes/DelphesClasses.h
rf53a4d2 r28027d5 84 84 //--------------------------------------------------------------------------- 85 85 86 class LHEFWeight: public TObject 87 { 88 public: 89 Int_t ID; // weight ID 90 Float_t Weight; // weight value 91 92 ClassDef(LHEFWeight, 1) 93 }; 94 95 //--------------------------------------------------------------------------- 96 86 97 class HepMCEvent: public Event 87 98 { … … 231 242 TRefArray Particles; // references to generated particles 232 243 244 // Isolation variables 245 246 Float_t IsolationVar; 247 Float_t IsolationVarRhoCorr; 248 Float_t SumPtCharged; 249 Float_t SumPtNeutral; 250 Float_t SumPtChargedPU; 251 Float_t SumPt; 252 233 253 static CompBase *fgCompare; //! 234 254 const CompBase *GetCompare() const { return fgCompare; } … … 257 277 TRef Particle; // reference to generated particle 258 278 279 // Isolation variables 280 281 Float_t IsolationVar; 282 Float_t IsolationVarRhoCorr; 283 Float_t SumPtCharged; 284 Float_t SumPtNeutral; 285 Float_t SumPtChargedPU; 286 Float_t SumPt; 287 259 288 static CompBase *fgCompare; //! 260 289 const CompBase *GetCompare() const { return fgCompare; } … … 281 310 TRef Particle; // reference to generated particle 282 311 312 // Isolation variables 313 314 Float_t IsolationVar; 315 Float_t IsolationVarRhoCorr; 316 Float_t SumPtCharged; 317 Float_t SumPtNeutral; 318 Float_t SumPtChargedPU; 319 Float_t SumPt; 320 283 321 static CompBase *fgCompare; //! 284 322 const CompBase *GetCompare() const { return fgCompare; } … … 306 344 Float_t DeltaPhi; // jet radius in azimuthal angle 307 345 346 UInt_t Flavor; 347 UInt_t FlavorAlgo; 348 UInt_t FlavorPhys; 349 308 350 UInt_t BTag; // 0 or 1 for a jet that has been tagged as containing a heavy quark 351 UInt_t BTagAlgo; 352 UInt_t BTagPhys; 353 309 354 UInt_t TauTag; // 0 or 1 for a jet that has been tagged as a tau 310 355 … … 313 358 Float_t EhadOverEem; // ratio of the hadronic versus electromagnetic energy deposited in the calorimeter 314 359 315 Int_t NCharged; // number of charged constituents 316 Int_t NNeutrals; // number of neutral constituents 317 Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents) 318 Float_t BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents) 319 Float_t MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent 320 Float_t PTD; // average pt between constituent and jet weighted by pt of constituent 321 Float_t FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents) 322 323 Float_t Tau1; // 1-subjettiness 324 Float_t Tau2; // 2-subjettiness 325 Float_t Tau3; // 3-subjettiness 326 Float_t Tau4; // 4-subjettiness 327 Float_t Tau5; // 5-subjettiness 360 Int_t NCharged; // number of charged constituents 361 Int_t NNeutrals; // number of neutral constituents 362 Float_t Beta; // (sum pt of charged pile-up constituents)/(sum pt of charged constituents) 363 Float_t BetaStar; // (sum pt of charged constituents coming from hard interaction)/(sum pt of charged constituents) 364 Float_t MeanSqDeltaR; // average distance (squared) between constituent and jet weighted by pt (squared) of constituent 365 Float_t PTD; // average pt between constituent and jet weighted by pt of constituent 366 Float_t FracPt[5]; // (sum pt of constituents within a ring 0.1*i < DeltaR < 0.1*(i+1))/(sum pt of constituents) 367 368 Float_t Tau[5]; // N-subjettiness 369 370 TLorentzVector TrimmedP4[5]; // first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta 371 TLorentzVector PrunedP4[5]; // first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta 372 TLorentzVector SoftDroppedP4[5]; // first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta 373 374 Int_t NSubJetsTrimmed; // number of subjets trimmed 375 Int_t NSubJetsPruned; // number of subjets pruned 376 Int_t NSubJetsSoftDropped; // number of subjets soft-dropped 328 377 329 378 TRefArray Constituents; // references to constituents … … 334 383 335 384 TLorentzVector P4() const; 336 337 ClassDef(Jet, 2) 385 TLorentzVector Area; 386 387 ClassDef(Jet, 3) 338 388 }; 339 389 … … 392 442 Float_t E; // calorimeter tower energy 393 443 394 Float_t T; //particle arrival time of flight 444 Float_t T; // ecal deposit time, averaged by sqrt(EM energy) over all particles, not smeared 445 Int_t Ntimes; // number of hits contributing to time measurement 395 446 396 447 Float_t Eem; // calorimeter tower electromagnetic energy … … 452 503 453 504 Int_t IsPU; 505 Int_t IsRecoPU; 506 454 507 Int_t IsConstituent; 455 508 509 UInt_t Flavor; 510 UInt_t FlavorAlgo; 511 UInt_t FlavorPhys; 512 456 513 UInt_t BTag; 514 UInt_t BTagAlgo; 515 UInt_t BTagPhys; 516 457 517 UInt_t TauTag; 458 518 … … 482 542 Float_t FracPt[5]; 483 543 544 //Timing information 545 546 Int_t Ntimes; 547 std::vector<std::pair<Float_t,Float_t> > Ecal_E_t; 548 549 // Isolation variables 550 551 Float_t IsolationVar; 552 Float_t IsolationVarRhoCorr; 553 Float_t SumPtCharged; 554 Float_t SumPtNeutral; 555 Float_t SumPtChargedPU; 556 Float_t SumPt; 557 484 558 // N-subjettiness variables 485 559 486 560 Float_t Tau[5]; 561 562 // Other Substructure variables 563 564 TLorentzVector TrimmedP4[5]; // first entry (i = 0) is the total Trimmed Jet 4-momenta and from i = 1 to 4 are the trimmed subjets 4-momenta 565 TLorentzVector PrunedP4[5]; // first entry (i = 0) is the total Pruned Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta 566 TLorentzVector SoftDroppedP4[5]; // first entry (i = 0) is the total SoftDropped Jet 4-momenta and from i = 1 to 4 are the pruned subjets 4-momenta 567 568 Int_t NSubJetsTrimmed; // number of subjets trimmed 569 Int_t NSubJetsPruned; // number of subjets pruned 570 Int_t NSubJetsSoftDropped; // number of subjets soft-dropped 571 487 572 488 573 static CompBase *fgCompare; //! … … 504 589 void SetFactory(DelphesFactory *factory) { fFactory = factory; } 505 590 506 ClassDef(Candidate, 2)591 ClassDef(Candidate, 3) 507 592 }; 508 593 -
classes/DelphesLHEFReader.cc
rf53a4d2 r28027d5 82 82 fEventCounter = -1; 83 83 fParticleCounter = -1; 84 f RwgtList.clear();84 fWeightList.clear(); 85 85 } 86 86 … … 99 99 TObjArray *partonOutputArray) 100 100 { 101 int rc ;101 int rc, id; 102 102 char *pch; 103 103 double weight; … … 158 158 else if(strstr(fBuffer, "<wgt")) 159 159 { 160 pch = str str(fBuffer, ">");160 pch = strpbrk(fBuffer, "\"'"); 161 161 if(!pch) 162 162 { … … 165 165 } 166 166 167 DelphesStream bufferStream(pch + 1); 168 rc = bufferStream.ReadDbl(weight); 167 DelphesStream idStream(pch + 1); 168 rc = idStream.ReadInt(id); 169 170 pch = strchr(fBuffer, '>'); 171 if(!pch) 172 { 173 cerr << "** ERROR: " << "invalid weight format" << endl; 174 return kFALSE; 175 } 176 177 DelphesStream weightStream(pch + 1); 178 rc = weightStream.ReadDbl(weight); 169 179 170 180 if(!rc) … … 174 184 } 175 185 176 f RwgtList.push_back(weight);186 fWeightList.push_back(make_pair(id, weight)); 177 187 } 178 188 else if(strstr(fBuffer, "</event>")) … … 206 216 //--------------------------------------------------------------------------- 207 217 208 void DelphesLHEFReader::AnalyzeRwgt(ExRootTreeBranch *branch) 209 { 210 Weight *element; 211 vector<double>::const_iterator itRwgtList; 212 213 for(itRwgtList = fRwgtList.begin(); itRwgtList != fRwgtList.end(); ++itRwgtList) 214 { 215 element = static_cast<Weight *>(branch->NewEntry()); 216 217 element->Weight = *itRwgtList; 218 void DelphesLHEFReader::AnalyzeWeight(ExRootTreeBranch *branch) 219 { 220 LHEFWeight *element; 221 vector< pair< int, double > >::const_iterator itWeightList; 222 223 for(itWeightList = fWeightList.begin(); itWeightList != fWeightList.end(); ++itWeightList) 224 { 225 element = static_cast<LHEFWeight *>(branch->NewEntry()); 226 227 element->ID = itWeightList->first; 228 element->Weight = itWeightList->second; 218 229 } 219 230 } -
classes/DelphesLHEFReader.h
rf53a4d2 r28027d5 31 31 32 32 #include <vector> 33 #include <utility> 33 34 34 35 class TObjArray; … … 58 59 TStopwatch *readStopWatch, TStopwatch *procStopWatch); 59 60 60 void Analyze Rwgt(ExRootTreeBranch *branch);61 void AnalyzeWeight(ExRootTreeBranch *branch); 61 62 62 63 private: … … 83 84 double fPx, fPy, fPz, fE, fMass; 84 85 85 std::vector< double> fRwgtList;86 std::vector< std::pair< int, double > > fWeightList; 86 87 }; 87 88 -
doc/genMakefile.tcl
rf53a4d2 r28027d5 283 283 dictDeps {DISPLAY_DICT} {display/DisplayLinkDef.h} 284 284 285 sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc} 285 sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc} {external/PUPPI/*.cc} 286 286 287 287 sourceDeps {FASTJET} {modules/FastJet*.cc} {external/fastjet/*.cc} {external/fastjet/tools/*.cc} {external/fastjet/plugins/*/*.cc} {external/fastjet/contribs/*/*.cc} -
modules/BTagging.cc
rf53a4d2 r28027d5 22 22 * Determines origin of jet, 23 23 * applies b-tagging efficiency (miss identification rate) formulas 24 * and sets b-tagging flags 24 * and sets b-tagging flags 25 25 * 26 26 * \author P. Demin - UCL, Louvain-la-Neuve … … 46 46 #include "TLorentzVector.h" 47 47 48 #include <algorithm> 48 #include <algorithm> 49 49 #include <stdexcept> 50 50 #include <iostream> … … 62 62 63 63 Int_t GetCategory(TObject *object); 64 64 65 65 Double_t fEtaMax, fPTMin; 66 66 }; … … 75 75 76 76 if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1; 77 77 78 78 pdgCode = TMath::Abs(parton->PID); 79 79 if(pdgCode != 21 && pdgCode > 5) return -1; … … 117 117 param = GetParam("EfficiencyFormula"); 118 118 size = param.GetSize(); 119 119 120 120 fEfficiencyMap.clear(); 121 121 for(i = 0; i < size/2; ++i) … … 143 143 144 144 fFilter = new ExRootFilter(fPartonInputArray); 145 145 146 146 fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets")); 147 147 fItJetInputArray = fJetInputArray->MakeIterator(); … … 170 170 void BTagging::Process() 171 171 { 172 Candidate *jet , *parton;172 Candidate *jet; 173 173 Double_t pt, eta, phi, e; 174 174 TObjArray *partonArray; 175 175 map< Int_t, DelphesFormula * >::iterator itEfficiencyMap; 176 176 DelphesFormula *formula; 177 Int_t pdgCode, pdgCodeMax;178 177 179 178 // select quark and gluons 180 179 fFilter->Reset(); 181 180 partonArray = fFilter->GetSubArray(fClassifier, 0); 182 181 183 182 if(partonArray == 0) return; 184 183 185 184 TIter itPartonArray(partonArray); 186 185 187 186 // loop over all input jets 188 187 fItJetInputArray->Reset(); … … 190 189 { 191 190 const TLorentzVector &jetMomentum = jet->Momentum; 192 pdgCodeMax = -1;193 191 eta = jetMomentum.Eta(); 194 192 phi = jetMomentum.Phi(); 195 193 pt = jetMomentum.Pt(); 196 194 e = jetMomentum.E(); 197 198 // loop over all input partons199 itPartonArray.Reset();200 while((parton = static_cast<Candidate*>(itPartonArray.Next())))201 {202 pdgCode = TMath::Abs(parton->PID);203 if(pdgCode == 21) pdgCode = 0;204 if(jetMomentum.DeltaR(parton->Momentum) <= fDeltaR)205 {206 if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;207 }208 }209 if(pdgCodeMax == 0) pdgCodeMax = 21;210 if(pdgCodeMax == -1) pdgCodeMax = 0;211 195 212 196 // find an efficency formula 213 itEfficiencyMap = fEfficiencyMap.find( pdgCodeMax);197 itEfficiencyMap = fEfficiencyMap.find(jet->Flavor); 214 198 if(itEfficiencyMap == fEfficiencyMap.end()) 215 199 { … … 220 204 // apply an efficency formula 221 205 jet->BTag |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber; 222 } 223 } 224 225 //------------------------------------------------------------------------------ 206 207 // find an efficency formula for algo flavor definition 208 itEfficiencyMap = fEfficiencyMap.find(jet->FlavorAlgo); 209 if(itEfficiencyMap == fEfficiencyMap.end()) 210 { 211 itEfficiencyMap = fEfficiencyMap.find(0); 212 } 213 formula = itEfficiencyMap->second; 214 215 // apply an efficency formula 216 jet->BTagAlgo |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber; 217 218 // find an efficency formula for phys flavor definition 219 itEfficiencyMap = fEfficiencyMap.find(jet->FlavorPhys); 220 if(itEfficiencyMap == fEfficiencyMap.end()) 221 { 222 itEfficiencyMap = fEfficiencyMap.find(0); 223 } 224 formula = itEfficiencyMap->second; 225 226 // apply an efficency formula 227 jet->BTagPhys |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber; 228 } 229 } 230 231 //------------------------------------------------------------------------------ -
modules/Calorimeter.cc
rf53a4d2 r28027d5 150 150 */ 151 151 152 // read min E value for timing measurement in ECAL 153 fTimingEMin = GetDouble("TimingEMin",4.); 154 // For timing 155 // So far this flag needs to be false 156 // Curved extrapolation not supported 157 fElectronsFromTrack = false; 158 159 152 160 // read min E value for towers to be saved 153 161 fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0); … … 356 364 fTrackHCalEnergy = 0.0; 357 365 358 fTowerECalTime = 0.0;359 fTowerHCalTime = 0.0;360 361 fTrackECalTime = 0.0;362 fTrackHCalTime = 0.0;363 364 fTowerECalTimeWeight = 0.0;365 fTowerHCalTimeWeight = 0.0;366 367 366 fTowerTrackHits = 0; 368 367 fTowerPhotonHits = 0; … … 380 379 position = track->Position; 381 380 382 383 381 ecalEnergy = momentum.E() * fTrackECalFractions[number]; 384 382 hcalEnergy = momentum.E() * fTrackHCalFractions[number]; … … 387 385 fTrackHCalEnergy += hcalEnergy; 388 386 389 fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T(); 390 fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T(); 391 392 fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy); 393 fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy); 387 bool dbg_scz = false; 388 if (dbg_scz) { 389 cout << " Calorimeter input track has x y z t " << track->Position.X() << " " << track->Position.Y() << " " << track->Position.Z() << " " << track->Position.T() 390 << endl; 391 Candidate *prt = static_cast<Candidate*>(track->GetCandidates()->Last()); 392 const TLorentzVector &ini = prt->Position; 393 394 cout << " and parent has x y z t " << ini.X() << " " << ini.Y() << " " << ini.Z() << " " << ini.T(); 395 396 } 397 398 if (ecalEnergy > fTimingEMin && fTower) { 399 if (fElectronsFromTrack) { 400 // cout << " SCZ Debug pushing back track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU 401 // << " PID=" << track->PID << endl; 402 fTower->Ecal_E_t.push_back(std::make_pair<float,float>(ecalEnergy,track->Position.T())); 403 } else { 404 // cout << " Skipping track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU 405 // << " PID=" << track->PID << endl; 406 } 407 } 408 394 409 395 410 fTowerTrackArray->Add(track); … … 397 412 continue; 398 413 } 399 414 400 415 // check for photon and electron hits in current tower 401 416 if(flags & 2) ++fTowerPhotonHits; … … 412 427 fTowerHCalEnergy += hcalEnergy; 413 428 414 fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T(); 415 fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T(); 416 417 fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy); 418 fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy); 419 429 if (ecalEnergy > fTimingEMin && fTower) { 430 if (abs(particle->PID) != 11 || !fElectronsFromTrack) { 431 // cout << " SCZ Debug About to push back particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU 432 // << " PID=" << particle->PID << endl; 433 fTower->Ecal_E_t.push_back(std::make_pair<Float_t,Float_t>(ecalEnergy,particle->Position.T())); 434 } else { 435 436 // N.B. Only charged particle set to leave ecal energy is the electrons 437 // cout << " SCZ Debug To avoid double-counting, skipping particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU 438 // << " PID=" << particle->PID << endl; 439 440 } 441 } 420 442 421 443 fTower->AddCandidate(particle); … … 434 456 Double_t ecalEnergy, hcalEnergy; 435 457 Double_t ecalSigma, hcalSigma; 436 Double_t ecalTime, hcalTime, time; 437 458 438 459 if(!fTower) return; 439 460 … … 444 465 hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma); 445 466 446 ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;447 hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;448 449 467 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy); 450 468 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy); … … 454 472 455 473 energy = ecalEnergy + hcalEnergy; 456 time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 457 474 458 475 if(fSmearTowerCenter) 459 476 { … … 469 486 pt = energy / TMath::CosH(eta); 470 487 471 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time); 488 // Time calculation for tower 489 fTower->Ntimes = 0; 490 Float_t tow_sumT = 0; 491 Float_t tow_sumW = 0; 492 493 for (Int_t i = 0 ; i < fTower->Ecal_E_t.size() ; i++) 494 { 495 Float_t w = TMath::Sqrt(fTower->Ecal_E_t[i].first); 496 tow_sumT += w*fTower->Ecal_E_t[i].second; 497 tow_sumW += w; 498 fTower->Ntimes++; 499 } 500 501 if (tow_sumW > 0.) { 502 fTower->Position.SetPtEtaPhiE(1.0, eta, phi,tow_sumT/tow_sumW); 503 } else { 504 fTower->Position.SetPtEtaPhiE(1.0,eta,phi,999999.); 505 } 506 507 472 508 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); 473 509 fTower->Eem = ecalEnergy; -
modules/Calorimeter.h
rf53a4d2 r28027d5 60 60 Double_t fTrackECalEnergy, fTrackHCalEnergy; 61 61 62 Double_t fTowerECalTime, fTowerHCalTime; 63 Double_t fTrackECalTime, fTrackHCalTime; 64 65 Double_t fTowerECalTimeWeight, fTowerHCalTimeWeight; 66 Double_t fTrackECalTimeWeight, fTrackHCalTimeWeight; 67 62 Double_t fTimingEMin; 63 Bool_t fElectronsFromTrack; // for timing 64 68 65 Int_t fTowerTrackHits, fTowerPhotonHits; 69 66 -
modules/FastJetFinder.cc
rf53a4d2 r28027d5 66 66 #include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh" 67 67 68 #include "fastjet/tools/Filter.hh" 69 #include "fastjet/tools/Pruner.hh" 70 #include "fastjet/contribs/RecursiveTools/SoftDrop.hh" 71 68 72 using namespace std; 69 73 using namespace fastjet; … … 121 125 fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8) 122 126 127 //-- Trimming parameters -- 128 129 fComputeTrimming = GetBool("ComputeTrimming", false); 130 fRTrim = GetDouble("RTrim", 0.2); 131 fPtFracTrim = GetDouble("PtFracTrim", 0.05); 132 133 134 //-- Pruning parameters -- 135 136 fComputePruning = GetBool("ComputePruning", false); 137 fZcutPrun = GetDouble("ZcutPrun", 0.1); 138 fRcutPrun = GetDouble("RcutPrun", 0.5); 139 fRPrun = GetDouble("RPrun", 0.8); 140 141 //-- SoftDrop parameters -- 142 143 fComputeSoftDrop = GetBool("ComputeSoftDrop", false); 144 fBetaSoftDrop = GetDouble("BetaSoftDrop", 0.0); 145 fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1); 146 fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8); 147 148 123 149 // --- Jet Area Parameters --- 124 150 fAreaAlgorithm = GetInt("AreaAlgorithm", 0); … … 260 286 PseudoJet jet, area; 261 287 ClusterSequence *sequence; 262 vector< PseudoJet > inputList, outputList ;288 vector< PseudoJet > inputList, outputList, subjets; 263 289 vector< PseudoJet >::iterator itInputList, itOutputList; 264 290 vector< TEstimatorStruct >::iterator itEstimators; … … 352 378 candidate->DeltaEta = detaMax; 353 379 candidate->DeltaPhi = dphiMax; 354 380 381 //------------------------------------ 382 // Trimming 383 //------------------------------------ 384 385 if(fComputeTrimming) 386 { 387 388 fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim)); 389 fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList); 390 391 trimmed_jet = join(trimmed_jet.constituents()); 392 393 candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m()); 394 395 // four hardest subjets 396 subjets.clear(); 397 subjets = trimmed_jet.pieces(); 398 subjets = sorted_by_pt(subjets); 399 400 candidate->NSubJetsTrimmed = subjets.size(); 401 402 for (size_t i = 0; i < subjets.size() and i < 4; i++){ 403 if(subjets.at(i).pt() < 0) continue ; 404 candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m()); 405 } 406 } 407 408 409 //------------------------------------ 410 // Pruning 411 //------------------------------------ 412 413 414 if(fComputePruning) 415 { 416 417 fastjet::Pruner pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm,fRPrun),fZcutPrun,fRcutPrun); 418 fastjet::PseudoJet pruned_jet = pruner(*itOutputList); 419 420 candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m()); 421 422 // four hardest subjet 423 subjets.clear(); 424 subjets = pruned_jet.pieces(); 425 subjets = sorted_by_pt(subjets); 426 427 candidate->NSubJetsPruned = subjets.size(); 428 429 for (size_t i = 0; i < subjets.size() and i < 4; i++){ 430 if(subjets.at(i).pt() < 0) continue ; 431 candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m()); 432 } 433 434 } 435 436 //------------------------------------ 437 // SoftDrop 438 //------------------------------------ 439 440 if(fComputeSoftDrop) 441 { 442 443 contrib::SoftDrop softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop); 444 fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList); 445 446 candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m()); 447 448 // four hardest subjet 449 450 subjets.clear(); 451 subjets = softdrop_jet.pieces(); 452 subjets = sorted_by_pt(subjets); 453 candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size(); 454 455 for (size_t i = 0; i < subjets.size() and i < 4; i++){ 456 if(subjets.at(i).pt() < 0) continue ; 457 candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m()); 458 } 459 } 460 355 461 // --- compute N-subjettiness with N = 1,2,3,4,5 ---- 356 462 -
modules/FastJetFinder.h
rf53a4d2 r28027d5 83 83 Int_t fN ; 84 84 85 //-- Trimming parameters -- 86 87 Bool_t fComputeTrimming; 88 Double_t fRTrim; 89 Double_t fPtFracTrim; 90 91 //-- Pruning parameters -- 92 93 Bool_t fComputePruning; 94 Double_t fZcutPrun; 95 Double_t fRcutPrun; 96 Double_t fRPrun; 97 98 //-- SoftDrop parameters -- 99 100 Bool_t fComputeSoftDrop; 101 Double_t fBetaSoftDrop; 102 Double_t fSymmetryCutSoftDrop; 103 Double_t fR0SoftDrop; 104 85 105 // --- FastJet Area method -------- 86 106 -
modules/Isolation.cc
rf53a4d2 r28027d5 25 25 * the transverse momenta fraction within (PTRatioMin, PTRatioMax]. 26 26 * 27 * \author P. Demin - UCL, Louvain-la-Neuve27 * \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve 28 28 * 29 29 */ … … 109 109 fUsePTSum = GetBool("UsePTSum", false); 110 110 111 fVetoLeptons = GetBool("VetoLeptons", true); 112 111 113 fClassifier->fPTMin = GetDouble("PTMin", 0.5); 114 112 115 113 116 // import input array(s) … … 153 156 Candidate *candidate, *isolation, *object; 154 157 TObjArray *isolationArray; 155 Double_t sum , ratio;158 Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU, sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr; 156 159 Int_t counter; 157 160 Double_t eta = 0.0; … … 180 183 181 184 // loop over all input tracks 182 sum = 0.0; 185 186 sumNeutral = 0.0; 187 sumCharged = 0.0; 188 sumChargedPU = 0.0; 189 sumAllParticles = 0.0; 190 183 191 counter = 0; 184 192 itIsolationArray.Reset(); 193 185 194 while((isolation = static_cast<Candidate*>(itIsolationArray.Next()))) 186 195 { … … 188 197 189 198 if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax && 190 candidate->GetUniqueID() != isolation->GetUniqueID()) 199 candidate->GetUniqueID() != isolation->GetUniqueID() && 200 ( !fVetoLeptons || (TMath::Abs(candidate->PID) != 11 && (TMath::Abs(candidate->PID) != 13)) ) ) 191 201 { 192 sum += isolationMomentum.Pt(); 202 203 sumAllParticles += isolationMomentum.Pt(); 204 if(isolation->Charge !=0) 205 { 206 sumCharged += isolationMomentum.Pt(); 207 if(isolation->IsRecoPU != 0) sumChargedPU += isolationMomentum.Pt(); 208 } 209 210 else sumNeutral += isolationMomentum.Pt(); 211 193 212 ++counter; 194 213 } … … 209 228 } 210 229 211 // correct sum for pile-up contamination 212 sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi(); 213 214 ratio = sum/candidateMomentum.Pt(); 215 if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue; 216 230 231 // correct sum for pile-up contamination 232 sumDBeta = sumCharged + TMath::Max(sumNeutral-0.5*sumChargedPU,0.0); 233 sumRhoCorr = sumCharged + TMath::Max(sumNeutral-TMath::Max(rho,0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(),0.0); 234 ratioDBeta = sumDBeta/candidateMomentum.Pt(); 235 ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt(); 236 237 candidate->IsolationVar = ratioDBeta; 238 candidate->IsolationVarRhoCorr = ratioRhoCorr; 239 candidate->SumPtCharged = sumCharged; 240 candidate->SumPtNeutral = sumNeutral; 241 candidate->SumPtChargedPU = sumChargedPU; 242 candidate->SumPt = sumAllParticles; 243 244 if((fUsePTSum && sumDBeta > fPTSumMax) || ratioDBeta > fPTRatioMax) continue; 217 245 fOutputArray->Add(candidate); 218 246 } -
modules/Isolation.h
rf53a4d2 r28027d5 59 59 Bool_t fUsePTSum; 60 60 61 Bool_t fVetoLeptons; 62 61 63 IsolationClassifier *fClassifier; //! 62 64 -
modules/ModulesLinkDef.h
rf53a4d2 r28027d5 58 58 #include "modules/Weighter.h" 59 59 #include "modules/Hector.h" 60 #include "modules/RunPUPPI.h" 61 #include "modules/JetFlavorAssociation.h" 60 62 #include "modules/ExampleModule.h" 61 63 … … 98 100 #pragma link C++ class Weighter+; 99 101 #pragma link C++ class Hector+; 102 #pragma link C++ class RunPUPPI+; 103 #pragma link C++ class JetFlavorAssociation+; 100 104 #pragma link C++ class ExampleModule+; 101 105 -
modules/PileUpJetID.cc
rf53a4d2 r28027d5 1 /*2 * Delphes: a framework for fast simulation of a generic collider experiment3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium4 *5 * This program is free software: you can redistribute it and/or modify6 * it under the terms of the GNU General Public License as published by7 * the Free Software Foundation, either version 3 of the License, or8 * (at your option) any later version.9 *10 * This program is distributed in the hope that it will be useful,11 * but WITHOUT ANY WARRANTY; without even the implied warranty of12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the13 * GNU General Public License for more details.14 *15 * You should have received a copy of the GNU General Public License16 * along with this program. If not, see <http://www.gnu.org/licenses/>.17 */18 19 1 20 2 /** \class PileUpJetID 21 3 * 22 * CMS PileUp Jet ID Variables , based on http://cds.cern.ch/record/15815834 * CMS PileUp Jet ID Variables 23 5 * 24 * \author S. Zenz , December 20136 * \author S. Zenz 25 7 * 26 8 */ … … 41 23 #include "TRandom3.h" 42 24 #include "TObjArray.h" 43 #include "TDatabasePDG.h"25 //#include "TDatabasePDG.h" 44 26 #include "TLorentzVector.h" 45 27 … … 54 36 55 37 PileUpJetID::PileUpJetID() : 56 fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0) ,fItVertexInputArray(0)38 fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0) 57 39 { 58 40 … … 70 52 void PileUpJetID::Init() 71 53 { 54 72 55 fJetPTMin = GetDouble("JetPTMin", 20.0); 73 56 fParameterR = GetDouble("ParameterR", 0.5); 74 57 fUseConstituents = GetInt("UseConstituents", 0); 75 58 59 60 /* 61 Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5 62 Double_t fBetaMinBarrel; // |eta| < 2.5 63 Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0 64 Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0 65 Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0 66 */ 67 68 fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1); 69 fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1); 70 fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1); 71 fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1); 72 fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1); 73 fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0); 74 fNeutralPTMin = GetDouble("NeutralPTMin", 2.0); 75 76 77 78 cout << " set MeanSqDeltaRMaxBarrel " << fMeanSqDeltaRMaxBarrel << endl; 79 cout << " set BetaMinBarrel " << fBetaMinBarrel << endl; 80 cout << " set MeanSqDeltaRMaxEndcap " << fMeanSqDeltaRMaxEndcap << endl; 81 cout << " set BetaMinEndcap " << fBetaMinEndcap << endl; 82 cout << " set MeanSqDeltaRMaxForward " << fMeanSqDeltaRMaxForward << endl; 83 84 85 76 86 fAverageEachTower = false; // for timing 77 87 … … 81 91 fItJetInputArray = fJetInputArray->MakeIterator(); 82 92 83 fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks")); 93 94 // cout << "BeforE SCZ additions in init" << endl; 95 // cout << GetString("TrackInputArray", "ParticlePropagator/tracks") << endl; 96 // cout << GetString("EFlowTrackInputArray", "ParticlePropagator/tracks") << endl; 97 98 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks")); 84 99 fItTrackInputArray = fTrackInputArray->MakeIterator(); 85 100 86 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", " Calorimeter/eflowTowers"));101 fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks")); 87 102 fItNeutralInputArray = fNeutralInputArray->MakeIterator(); 88 103 89 fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));90 fItVertexInputArray = fVertexInputArray->MakeIterator();91 92 fZVertexResolution = GetDouble("ZVertexResolution", 0.005)*1.0E3;93 104 94 105 // create output array(s) 95 106 96 107 fOutputArray = ExportArray(GetString("OutputArray", "jets")); 108 109 fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers")); 110 111 112 // cout << " end of INIT " << endl; 113 97 114 } 98 115 … … 101 118 void PileUpJetID::Finish() 102 119 { 120 // cout << "In finish" << endl; 121 103 122 if(fItJetInputArray) delete fItJetInputArray; 104 123 if(fItTrackInputArray) delete fItTrackInputArray; 105 124 if(fItNeutralInputArray) delete fItNeutralInputArray; 106 if(fItVertexInputArray) delete fItVertexInputArray; 125 107 126 } 108 127 … … 111 130 void PileUpJetID::Process() 112 131 { 132 // cout << "start of process" << endl; 133 113 134 Candidate *candidate, *constituent; 114 135 TLorentzVector momentum, area; 115 Int_t i, nc, nn; 116 Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq; 117 Double_t dr, pt, pt_ann[5]; 118 Double_t zvtx = 0.0; 119 120 Candidate *track; 121 122 // find z position of primary vertex 123 124 fItVertexInputArray->Reset(); 125 while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next()))) 126 { 127 if(!candidate->IsPU) 128 { 129 zvtx = candidate->Position.Z(); 130 break; 131 } 132 } 136 137 // cout << "BeforE SCZ additions in process" << endl; 138 139 // SCZ 140 Candidate *trk; 133 141 134 142 // loop over all input candidates … … 139 147 area = candidate->Area; 140 148 141 sumpt = 0.0; 142 sumptch = 0.0; 143 sumptchpv = 0.0; 144 sumptchpu = 0.0; 145 sumdrsqptsq = 0.0; 146 sumptsq = 0.0; 147 nc = 0; 148 nn = 0; 149 150 for(i = 0; i < 5; ++i) 151 { 152 pt_ann[i] = 0.0; 153 } 154 155 if(fUseConstituents) 156 { 149 float sumT0 = 0.; 150 float sumT1 = 0.; 151 float sumT10 = 0.; 152 float sumT20 = 0.; 153 float sumT30 = 0.; 154 float sumT40 = 0.; 155 float sumWeightsForT = 0.; 156 candidate->Ntimes = 0; 157 158 float sumpt = 0.; 159 float sumptch = 0.; 160 float sumptchpv = 0.; 161 float sumptchpu = 0.; 162 float sumdrsqptsq = 0.; 163 float sumptsq = 0.; 164 int nc = 0; 165 int nn = 0; 166 float pt_ann[5]; 167 168 for (int i = 0 ; i < 5 ; i++) { 169 pt_ann[i] = 0.; 170 } 171 172 if (fUseConstituents) { 157 173 TIter itConstituents(candidate->GetCandidates()); 158 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) 159 { 160 pt = constituent->Momentum.Pt(); 161 dr = candidate->Momentum.DeltaR(constituent->Momentum); 162 sumpt += pt; 163 sumdrsqptsq += dr*dr*pt*pt; 164 sumptsq += pt*pt; 165 if(constituent->Charge == 0) 166 { 167 // neutrals 168 ++nn; 169 } 170 else 171 { 172 // charged 173 if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution) 174 { 175 sumptchpu += pt; 176 } 177 else 178 { 179 sumptchpv += pt; 180 } 181 sumptch += pt; 182 ++nc; 183 } 184 for(i = 0; i < 5; ++i) 185 { 186 if(dr > 0.1*i && dr < 0.1*(i + 1)) 187 { 188 pt_ann[i] += pt; 189 } 190 } 191 } 192 } 193 else 194 { 174 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) { 175 float pt = constituent->Momentum.Pt(); 176 float dr = candidate->Momentum.DeltaR(constituent->Momentum); 177 // cout << " There exists a constituent with dr=" << dr << endl; 178 sumpt += pt; 179 sumdrsqptsq += dr*dr*pt*pt; 180 sumptsq += pt*pt; 181 if (constituent->Charge == 0) { 182 nn++; 183 } else { 184 if (constituent->IsRecoPU) { 185 sumptchpu += pt; 186 } else { 187 sumptchpv += pt; 188 } 189 sumptch += pt; 190 nc++; 191 } 192 for (int i = 0 ; i < 5 ; i++) { 193 if (dr > 0.1*i && dr < 0.1*(i+1)) { 194 pt_ann[i] += pt; 195 } 196 } 197 float tow_sumT = 0; 198 float tow_sumW = 0; 199 for (int i = 0 ; i < constituent->Ecal_E_t.size() ; i++) { 200 float w = TMath::Sqrt(constituent->Ecal_E_t[i].first); 201 if (fAverageEachTower) { 202 tow_sumT += w*constituent->Ecal_E_t[i].second; 203 tow_sumW += w; 204 } else { 205 sumT0 += w*constituent->Ecal_E_t[i].second; 206 sumT1 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.001); 207 sumT10 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.010); 208 sumT20 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.020); 209 sumT30 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.030); 210 sumT40 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.040); 211 sumWeightsForT += w; 212 candidate->Ntimes++; 213 } 214 } 215 if (fAverageEachTower && tow_sumW > 0.) { 216 sumT0 += tow_sumT; 217 sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001); 218 sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010); 219 sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020); 220 sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030); 221 sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040); 222 sumWeightsForT += tow_sumW; 223 candidate->Ntimes++; 224 } 225 } 226 } else { 195 227 // Not using constituents, using dr 196 228 fItTrackInputArray->Reset(); 197 while((track = static_cast<Candidate*>(fItTrackInputArray->Next()))) 198 { 199 if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR) 200 { 201 pt = track->Momentum.Pt(); 202 sumpt += pt; 203 sumptch += pt; 204 if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution) 205 { 206 sumptchpu += pt; 207 } 208 else 209 { 210 sumptchpv += pt; 211 } 212 dr = candidate->Momentum.DeltaR(track->Momentum); 213 sumdrsqptsq += dr*dr*pt*pt; 214 sumptsq += pt*pt; 215 nc++; 216 for(i = 0; i < 5; ++i) 217 { 218 if(dr > 0.1*i && dr < 0.1*(i + 1)) 219 { 220 pt_ann[i] += pt; 221 } 222 } 223 } 224 } 225 229 while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) { 230 if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) { 231 float pt = trk->Momentum.Pt(); 232 sumpt += pt; 233 sumptch += pt; 234 if (trk->IsRecoPU) { 235 sumptchpu += pt; 236 } else { 237 sumptchpv += pt; 238 } 239 float dr = candidate->Momentum.DeltaR(trk->Momentum); 240 sumdrsqptsq += dr*dr*pt*pt; 241 sumptsq += pt*pt; 242 nc++; 243 for (int i = 0 ; i < 5 ; i++) { 244 if (dr > 0.1*i && dr < 0.1*(i+1)) { 245 pt_ann[i] += pt; 246 } 247 } 248 } 249 } 226 250 fItNeutralInputArray->Reset(); 227 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) 228 { 229 if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) 230 { 231 pt = constituent->Momentum.Pt(); 232 sumpt += pt; 233 dr = candidate->Momentum.DeltaR(constituent->Momentum); 234 sumdrsqptsq += dr*dr*pt*pt; 235 sumptsq += pt*pt; 236 nn++; 237 for(i = 0; i < 5; ++i) 238 { 239 if(dr > 0.1*i && dr < 0.1*(i + 1)) 240 { 241 pt_ann[i] += pt; 242 } 243 } 244 } 245 } 246 } 247 248 if(sumptch > 0.0) 249 { 250 candidate->Beta = sumptchpu/sumptch; 251 candidate->BetaStar = sumptchpv/sumptch; 252 } 253 else 254 { 255 candidate->Beta = -999.0; 256 candidate->BetaStar = -999.0; 257 } 258 if(sumptsq > 0.0) 259 { 251 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) { 252 if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) { 253 float pt = constituent->Momentum.Pt(); 254 sumpt += pt; 255 float dr = candidate->Momentum.DeltaR(constituent->Momentum); 256 sumdrsqptsq += dr*dr*pt*pt; 257 sumptsq += pt*pt; 258 nn++; 259 for (int i = 0 ; i < 5 ; i++) { 260 if (dr > 0.1*i && dr < 0.1*(i+1)) { 261 pt_ann[i] += pt; 262 } 263 } 264 } 265 } 266 } 267 268 if (sumptch > 0.) { 269 candidate->Beta = sumptchpv/sumptch; 270 candidate->BetaStar = sumptchpu/sumptch; 271 } else { 272 candidate->Beta = -999.; 273 candidate->BetaStar = -999.; 274 } 275 if (sumptsq > 0.) { 260 276 candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq; 261 } 262 else 263 { 264 candidate->MeanSqDeltaR = -999.0; 277 } else { 278 candidate->MeanSqDeltaR = -999.; 265 279 } 266 280 candidate->NCharged = nc; 267 281 candidate->NNeutrals = nn; 268 if(sumpt > 0.0) 269 { 282 if (sumpt > 0.) { 270 283 candidate->PTD = TMath::Sqrt(sumptsq) / sumpt; 271 for(i = 0; i < 5; ++i) 272 { 284 for (int i = 0 ; i < 5 ; i++) { 273 285 candidate->FracPt[i] = pt_ann[i]/sumpt; 274 286 } 275 } 276 else 277 { 278 candidate->PTD = -999.0; 279 for(i = 0; i < 5; ++i) 280 { 281 candidate->FracPt[i] = -999.0; 287 } else { 288 candidate->PTD = -999.; 289 for (int i = 0 ; i < 5 ; i++) { 290 candidate->FracPt[i] = -999.; 282 291 } 283 292 } 284 293 285 294 fOutputArray->Add(candidate); 295 296 // New stuff 297 /* 298 fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1); 299 fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1); 300 fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1); 301 fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1); 302 fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1); 303 */ 304 305 bool passId = false; 306 if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) { 307 if (fabs(candidate->Momentum.Eta())<1.5) { 308 passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel)); 309 } else if (fabs(candidate->Momentum.Eta())<4.0) { 310 passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap)); 311 } else { 312 passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward); 313 } 314 } 315 316 // cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt() 317 // << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl; 318 319 if (passId) { 320 if (fUseConstituents) { 321 TIter itConstituents(candidate->GetCandidates()); 322 while((constituent = static_cast<Candidate*>(itConstituents.Next()))) { 323 if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) { 324 fNeutralsInPassingJets->Add(constituent); 325 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl; 326 } 327 } 328 } else { // use DeltaR 329 fItNeutralInputArray->Reset(); 330 while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) { 331 if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) { 332 fNeutralsInPassingJets->Add(constituent); 333 // cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl; 334 } 335 } 336 } 337 } 338 339 286 340 } 287 341 } 288 342 289 343 //------------------------------------------------------------------------------ 290 -
modules/PileUpJetID.h
rf53a4d2 r28027d5 1 /*2 * Delphes: a framework for fast simulation of a generic collider experiment3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium4 *5 * This program is free software: you can redistribute it and/or modify6 * it under the terms of the GNU General Public License as published by7 * the Free Software Foundation, either version 3 of the License, or8 * (at your option) any later version.9 *10 * This program is distributed in the hope that it will be useful,11 * but WITHOUT ANY WARRANTY; without even the implied warranty of12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the13 * GNU General Public License for more details.14 *15 * You should have received a copy of the GNU General Public License16 * along with this program. If not, see <http://www.gnu.org/licenses/>.17 */18 19 1 #ifndef PileUpJetID_h 20 2 #define PileUpJetID_h … … 22 4 /** \class PileUpJetID 23 5 * 24 * CMS PileUp Jet ID Variables , based on http://cds.cern.ch/record/15815836 * CMS PileUp Jet ID Variables 25 7 * 26 * \author S. Zenz , December 20138 * \author S. Zenz 27 9 * 28 10 */ … … 34 16 35 17 class TObjArray; 18 class DelphesFormula; 36 19 37 20 class PileUpJetID: public DelphesModule … … 51 34 Double_t fParameterR; 52 35 36 Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5 37 Double_t fBetaMinBarrel; // |eta| < 2.5 38 Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0 39 Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0 40 Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0 41 42 Double_t fNeutralPTMin; 43 Double_t fJetPTMinForNeutrals; 44 45 /* 46 JAY 47 --- 48 49 |Eta|<1.5 50 51 meanSqDeltaR betaStar SigEff BgdEff 52 0.13 0.92 96% 8% 53 0.13 0.95 97% 16% 54 0.13 0.97 98% 27% 55 56 |Eta|>1.5 57 58 meanSqDeltaR betaStar SigEff BgdEff 59 0.14 0.91 95% 15% 60 0.14 0.94 97% 19% 61 0.14 0.97 98% 29% 62 63 BRYAN 64 ----- 65 66 Barrel (MeanSqDR, Beta, sig eff, bg eff): 67 0.10, 0.08, 90%, 8% 68 0.11, 0.12, 90%, 6% 69 0.13, 0.16, 89%, 5% 70 71 Endcap (MeanSqDR, Beta, sig eff, bg eff): 72 0.07, 0.06, 89%, 4% 73 0.08, 0.08, 92%, 6% 74 0.09, 0.08, 95%, 10% 75 0.10, 0.08, 97%, 13% 76 77 SETH GUESSES FOR |eta| > 4.0 78 ---------------------------- 79 80 MeanSqDeltaR 81 0.07 82 0.10 83 0.14 84 0.2 85 */ 86 53 87 // If set to true, may have weird results for PFCHS 54 88 // If set to false, uses everything within dR < fParameterR even if in other jets &c. 55 89 // Results should be very similar for PF 56 Int_t fUseConstituents; 90 Int_t fUseConstituents; 57 91 58 92 Bool_t fAverageEachTower; … … 62 96 const TObjArray *fJetInputArray; //! 63 97 64 const TObjArray *fTrackInputArray; // !65 const TObjArray *fNeutralInputArray; //!98 const TObjArray *fTrackInputArray; // SCZ 99 const TObjArray *fNeutralInputArray; 66 100 67 TIterator *fItTrackInputArray; // !68 TIterator *fItNeutralInputArray; // !101 TIterator *fItTrackInputArray; // SCZ 102 TIterator *fItNeutralInputArray; // SCZ 69 103 70 104 TObjArray *fOutputArray; //! 105 TObjArray *fNeutralsInPassingJets; // SCZ 71 106 72 TIterator *fItVertexInputArray; //!73 const TObjArray *fVertexInputArray; //!74 107 75 Double_t fZVertexResolution; 76 77 ClassDef(PileUpJetID, 1) 108 ClassDef(PileUpJetID, 2) 78 109 }; 79 110 80 111 #endif 81 -
modules/PileUpMerger.cc
rf53a4d2 r28027d5 80 80 fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09); 81 81 82 fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0); 83 fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0); 84 fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0); 85 fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0); 86 82 87 // read vertex smearing formula 83 88 … … 111 116 TParticlePDG *pdgParticle; 112 117 Int_t pid; 113 Float_t x, y, z, t ;118 Float_t x, y, z, t, vx, vy; 114 119 Float_t px, py, pz, e; 115 120 Double_t dz, dphi, dt; 116 Int_t numberOfEvents, event ;121 Int_t numberOfEvents, event, numberOfParticles; 117 122 Long64_t allEntries, entry; 118 Candidate *candidate, *vertex candidate;123 Candidate *candidate, *vertex; 119 124 DelphesFactory *factory; 120 125 … … 123 128 fItInputArray->Reset(); 124 129 125 // --- Deal with Primary vertex first ------130 // --- Deal with primary vertex first ------ 126 131 127 132 fFunction->GetRandom2(dz, dt); … … 129 134 dt *= c_light*1.0E3; // necessary in order to make t in mm/c 130 135 dz *= 1.0E3; // necessary in order to make z in mm 131 136 vx = 0.0; 137 vy = 0.0; 138 numberOfParticles = fInputArray->GetEntriesFast(); 132 139 while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) 133 140 { 141 vx += candidate->Position.X(); 142 vy += candidate->Position.Y(); 134 143 z = candidate->Position.Z(); 135 144 t = candidate->Position.T(); … … 139 148 } 140 149 150 if(numberOfParticles > 0) 151 { 152 vx /= numberOfParticles; 153 vy /= numberOfParticles; 154 } 155 141 156 factory = GetFactory(); 142 157 143 vertex candidate= factory->NewCandidate();144 vertex candidate->Position.SetXYZT(0.0, 0.0, dz, dt);145 fVertexOutputArray->Add(vertex candidate);158 vertex = factory->NewCandidate(); 159 vertex->Position.SetXYZT(vx, vy, dz, dt); 160 fVertexOutputArray->Add(vertex); 146 161 147 162 // --- Then with pile-up vertices ------ … … 181 196 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi()); 182 197 183 vertexcandidate = factory->NewCandidate(); 184 vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt); 185 vertexcandidate->IsPU = 1; 186 187 fVertexOutputArray->Add(vertexcandidate); 188 198 vx = 0.0; 199 vy = 0.0; 200 numberOfParticles = 0; 189 201 while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e)) 190 202 { … … 204 216 candidate->Momentum.RotateZ(dphi); 205 217 218 x -= fInputBeamSpotX; 219 y -= fInputBeamSpotY; 206 220 candidate->Position.SetXYZT(x, y, z + dz, t + dt); 207 221 candidate->Position.RotateZ(dphi); 222 candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0); 223 224 vx += candidate->Position.X(); 225 vy += candidate->Position.Y(); 226 ++numberOfParticles; 208 227 209 228 fParticleOutputArray->Add(candidate); 210 229 } 211 } 212 } 213 214 //------------------------------------------------------------------------------ 215 230 231 if(numberOfParticles > 0) 232 { 233 vx /= numberOfParticles; 234 vy /= numberOfParticles; 235 } 236 237 vertex = factory->NewCandidate(); 238 vertex->Position.SetXYZT(vx, vy, dz, dt); 239 vertex->IsPU = 1; 240 241 fVertexOutputArray->Add(vertex); 242 } 243 } 244 245 //------------------------------------------------------------------------------ -
modules/PileUpMerger.h
rf53a4d2 r28027d5 53 53 Double_t fTVertexSpread; 54 54 55 Double_t fInputBeamSpotX; 56 Double_t fInputBeamSpotY; 57 Double_t fOutputBeamSpotX; 58 Double_t fOutputBeamSpotY; 59 55 60 DelphesTF2 *fFunction; //! 56 61 -
modules/PileUpMergerPythia8.cc
rf53a4d2 r28027d5 29 29 #include "classes/DelphesClasses.h" 30 30 #include "classes/DelphesFactory.h" 31 #include "classes/Delphes Formula.h"31 #include "classes/DelphesTF2.h" 32 32 #include "classes/DelphesPileUpReader.h" 33 33 … … 56 56 57 57 PileUpMergerPythia8::PileUpMergerPythia8() : 58 fPythia(0), fItInputArray(0) 59 { 58 fFunction(0), fPythia(0), fItInputArray(0) 59 { 60 fFunction = new DelphesTF2; 60 61 } 61 62 … … 64 65 PileUpMergerPythia8::~PileUpMergerPythia8() 65 66 { 67 delete fFunction; 66 68 } 67 69 … … 72 74 const char *fileName; 73 75 76 fPileUpDistribution = GetInt("PileUpDistribution", 0); 77 74 78 fMeanPileUp = GetDouble("MeanPileUp", 10); 75 fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3; 79 80 fZVertexSpread = GetDouble("ZVertexSpread", 0.15); 81 fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09); 82 83 fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0); 84 fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0); 85 fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0); 86 fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0); 76 87 77 88 fPTMin = GetDouble("PTMin", 0.0); 89 90 fFunction->Compile(GetString("VertexDistributionFormula", "0.0")); 91 fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread); 78 92 79 93 fileName = GetString("ConfigFile", "MinBias.cmnd"); … … 86 100 87 101 // create output arrays 88 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles")); 102 fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles")); 103 fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices")); 89 104 } 90 105 … … 103 118 TParticlePDG *pdgParticle; 104 119 Int_t pid, status; 105 Float_t x, y, z, t ;120 Float_t x, y, z, t, vx, vy; 106 121 Float_t px, py, pz, e; 107 Double_t dz, dphi ;108 Int_t poisson, event, i;109 Candidate *candidate ;122 Double_t dz, dphi, dt; 123 Int_t numberOfEvents, event, numberOfParticles, i; 124 Candidate *candidate, *vertex; 110 125 DelphesFactory *factory; 111 126 127 const Double_t c_light = 2.99792458E8; 128 112 129 fItInputArray->Reset(); 130 131 // --- Deal with primary vertex first ------ 132 133 fFunction->GetRandom2(dz, dt); 134 135 dt *= c_light*1.0E3; // necessary in order to make t in mm/c 136 dz *= 1.0E3; // necessary in order to make z in mm 137 vx = 0.0; 138 vy = 0.0; 139 numberOfParticles = fInputArray->GetEntriesFast(); 113 140 while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) 114 141 { 115 fOutputArray->Add(candidate); 142 vx += candidate->Position.X(); 143 vy += candidate->Position.Y(); 144 z = candidate->Position.Z(); 145 t = candidate->Position.T(); 146 candidate->Position.SetZ(z + dz); 147 candidate->Position.SetT(t + dt); 148 fParticleOutputArray->Add(candidate); 149 } 150 151 if(numberOfParticles > 0) 152 { 153 vx /= numberOfParticles; 154 vy /= numberOfParticles; 116 155 } 117 156 118 157 factory = GetFactory(); 119 158 120 poisson = gRandom->Poisson(fMeanPileUp); 121 122 for(event = 0; event < poisson; ++event) 159 vertex = factory->NewCandidate(); 160 vertex->Position.SetXYZT(vx, vy, dz, dt); 161 fVertexOutputArray->Add(vertex); 162 163 // --- Then with pile-up vertices ------ 164 165 switch(fPileUpDistribution) 166 { 167 case 0: 168 numberOfEvents = gRandom->Poisson(fMeanPileUp); 169 break; 170 case 1: 171 numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1); 172 break; 173 default: 174 numberOfEvents = gRandom->Poisson(fMeanPileUp); 175 break; 176 } 177 178 for(event = 0; event < numberOfEvents; ++event) 123 179 { 124 180 while(!fPythia->next()); 125 181 126 dz = gRandom->Gaus(0.0, fZVertexSpread); 182 // --- Pile-up vertex smearing 183 184 fFunction->GetRandom2(dz, dt); 185 186 dt *= c_light*1.0E3; // necessary in order to make t in mm/c 187 dz *= 1.0E3; // necessary in order to make z in mm 188 127 189 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi()); 128 190 129 for(i = 1; i < fPythia->event.size(); ++i) 191 vx = 0.0; 192 vy = 0.0; 193 numberOfParticles = fPythia->event.size(); 194 for(i = 1; i < numberOfParticles; ++i) 130 195 { 131 196 Pythia8::Particle &particle = fPythia->event[i]; … … 143 208 candidate->PID = pid; 144 209 145 candidate->Status = status;210 candidate->Status = 1; 146 211 147 212 pdgParticle = pdg->GetParticle(pid); … … 154 219 candidate->Momentum.RotateZ(dphi); 155 220 156 candidate->Position.SetXYZT(x, y, z + dz, t); 221 x -= fInputBeamSpotX; 222 y -= fInputBeamSpotY; 223 candidate->Position.SetXYZT(x, y, z + dz, t + dt); 157 224 candidate->Position.RotateZ(dphi); 158 159 fOutputArray->Add(candidate); 225 candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0); 226 227 vx += candidate->Position.X(); 228 vy += candidate->Position.Y(); 229 230 fParticleOutputArray->Add(candidate); 160 231 } 161 } 162 } 163 164 //------------------------------------------------------------------------------ 165 232 233 if(numberOfParticles > 0) 234 { 235 vx /= numberOfParticles; 236 vy /= numberOfParticles; 237 } 238 239 vertex = factory->NewCandidate(); 240 vertex->Position.SetXYZT(vx, vy, dz, dt); 241 vertex->IsPU = 1; 242 243 fVertexOutputArray->Add(vertex); 244 } 245 } 246 247 //------------------------------------------------------------------------------ -
modules/PileUpMergerPythia8.h
rf53a4d2 r28027d5 31 31 32 32 class TObjArray; 33 class DelphesTF2; 33 34 34 35 namespace Pythia8 … … 50 51 private: 51 52 53 Int_t fPileUpDistribution; 52 54 Double_t fMeanPileUp; 55 53 56 Double_t fZVertexSpread; 57 Double_t fTVertexSpread; 58 59 Double_t fInputBeamSpotX; 60 Double_t fInputBeamSpotY; 61 Double_t fOutputBeamSpotX; 62 Double_t fOutputBeamSpotY; 63 54 64 Double_t fPTMin; 65 66 DelphesTF2 *fFunction; //! 55 67 56 68 Pythia8::Pythia *fPythia; //! … … 60 72 const TObjArray *fInputArray; //! 61 73 62 TObjArray *fOutputArray; //! 74 TObjArray *fParticleOutputArray; //! 75 TObjArray *fVertexOutputArray; //! 63 76 64 77 ClassDef(PileUpMergerPythia8, 1) -
modules/TrackPileUpSubtractor.cc
rf53a4d2 r28027d5 74 74 fZVertexResolution = GetDouble("ZVertexResolution", 0.005)*1.0E3; 75 75 76 fPTMin = GetDouble("PTMin", 0.); 76 77 // import arrays with output from other modules 77 78 78 79 ExRootConfParam param = GetParam("InputArray"); 79 80 Long_t i, size; … … 147 148 // assume perfect pile-up subtraction for tracks outside fZVertexResolution 148 149 149 if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue; 150 151 array->Add(candidate); 150 if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) candidate->IsRecoPU = 1; 151 else 152 { 153 candidate->IsRecoPU = 0; 154 if( candidate->Momentum.Pt() > fPTMin) array->Add(candidate); 155 } 152 156 } 153 157 } -
modules/TrackPileUpSubtractor.h
rf53a4d2 r28027d5 50 50 Double_t fZVertexResolution; 51 51 52 Double_t fPTMin; 53 52 54 std::map< TIterator *, TObjArray * > fInputMap; //! 53 55 -
modules/TreeWriter.cc
rf53a4d2 r28027d5 291 291 entry->ZOuter = position.Z(); 292 292 entry->TOuter = position.T()*1.0E-3/c_light; 293 293 294 294 entry->Dxy = candidate->Dxy; 295 295 entry->SDxy = candidate->SDxy ; … … 297 297 entry->Yd = candidate->Yd; 298 298 entry->Zd = candidate->Zd; 299 299 300 300 const TLorentzVector &momentum = candidate->Momentum; 301 301 … … 362 362 363 363 entry->T = position.T()*1.0E-3/c_light; 364 entry->Ntimes = candidate->Ntimes; 364 365 365 366 FillParticles(candidate, &entry->Particles); … … 403 404 entry->T = position.T()*1.0E-3/c_light; 404 405 406 // Isolation variables 407 408 entry->IsolationVar = candidate->IsolationVar; 409 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ; 410 entry->SumPtCharged = candidate->SumPtCharged ; 411 entry->SumPtNeutral = candidate->SumPtNeutral ; 412 entry->SumPtChargedPU = candidate->SumPtChargedPU ; 413 entry->SumPt = candidate->SumPt ; 414 405 415 entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9; 406 416 … … 442 452 entry->T = position.T()*1.0E-3/c_light; 443 453 454 // Isolation variables 455 456 entry->IsolationVar = candidate->IsolationVar; 457 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ; 458 entry->SumPtCharged = candidate->SumPtCharged ; 459 entry->SumPtNeutral = candidate->SumPtNeutral ; 460 entry->SumPtChargedPU = candidate->SumPtChargedPU ; 461 entry->SumPt = candidate->SumPt ; 462 463 444 464 entry->Charge = candidate->Charge; 445 465 … … 469 489 const TLorentzVector &momentum = candidate->Momentum; 470 490 const TLorentzVector &position = candidate->Position; 471 472 491 473 492 pt = momentum.Pt(); … … 488 507 entry->T = position.T()*1.0E-3/c_light; 489 508 509 // Isolation variables 510 511 entry->IsolationVar = candidate->IsolationVar; 512 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ; 513 entry->SumPtCharged = candidate->SumPtCharged ; 514 entry->SumPtNeutral = candidate->SumPtNeutral ; 515 entry->SumPtChargedPU = candidate->SumPtChargedPU ; 516 entry->SumPt = candidate->SumPt ; 517 490 518 entry->Charge = candidate->Charge; 491 519 … … 504 532 Double_t ecalEnergy, hcalEnergy; 505 533 const Double_t c_light = 2.99792458E8; 534 Int_t i; 506 535 507 536 array->Sort(); … … 532 561 entry->Mass = momentum.M(); 533 562 563 entry->Area = candidate->Area; 564 534 565 entry->DeltaEta = candidate->DeltaEta; 535 566 entry->DeltaPhi = candidate->DeltaPhi; 536 567 568 entry->Flavor = candidate->Flavor; 569 entry->FlavorAlgo = candidate->FlavorAlgo; 570 entry->FlavorPhys = candidate->FlavorPhys; 571 537 572 entry->BTag = candidate->BTag; 573 574 entry->BTagAlgo = candidate->BTagAlgo; 575 entry->BTagPhys = candidate->BTagPhys; 576 538 577 entry->TauTag = candidate->TauTag; 539 578 … … 561 600 entry->MeanSqDeltaR = candidate->MeanSqDeltaR; 562 601 entry->PTD = candidate->PTD; 563 entry->FracPt[0] = candidate->FracPt[0]; 564 entry->FracPt[1] = candidate->FracPt[1]; 565 entry->FracPt[2] = candidate->FracPt[2]; 566 entry->FracPt[3] = candidate->FracPt[3]; 567 entry->FracPt[4] = candidate->FracPt[4]; 568 569 //--- N-subjettiness variables ---- 570 571 entry->Tau1 = candidate->Tau[0]; 572 entry->Tau2 = candidate->Tau[1]; 573 entry->Tau3 = candidate->Tau[2]; 574 entry->Tau4 = candidate->Tau[3]; 575 entry->Tau5 = candidate->Tau[4]; 576 602 603 //--- Sub-structure variables ---- 604 605 entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed; 606 entry->NSubJetsPruned = candidate->NSubJetsPruned; 607 entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped; 608 609 for(i = 0; i < 5; i++) 610 { 611 entry->FracPt[i] = candidate -> FracPt[i]; 612 entry->Tau[i] = candidate -> Tau[i]; 613 entry->TrimmedP4[i] = candidate -> TrimmedP4[i]; 614 entry->PrunedP4[i] = candidate -> PrunedP4[i]; 615 entry->SoftDroppedP4[i] = candidate -> SoftDroppedP4[i]; 616 } 617 577 618 FillParticles(candidate, &entry->Particles); 578 619 } -
readers/DelphesLHEF.cpp
rf53a4d2 r28027d5 63 63 TStopwatch readStopWatch, procStopWatch; 64 64 ExRootTreeWriter *treeWriter = 0; 65 ExRootTreeBranch *branchEvent = 0, *branch Rwgt = 0;65 ExRootTreeBranch *branchEvent = 0, *branchWeight = 0; 66 66 ExRootConfReader *confReader = 0; 67 67 Delphes *modularDelphes = 0; … … 103 103 104 104 branchEvent = treeWriter->NewBranch("Event", LHEFEvent::Class()); 105 branch Rwgt = treeWriter->NewBranch("Rwgt", Weight::Class());105 branchWeight = treeWriter->NewBranch("Weight", Weight::Class()); 106 106 107 107 confReader = new ExRootConfReader; … … 196 196 197 197 reader->AnalyzeEvent(branchEvent, eventCounter, &readStopWatch, &procStopWatch); 198 reader->Analyze Rwgt(branchRwgt);198 reader->AnalyzeWeight(branchWeight); 199 199 200 200 treeWriter->Fill(); -
readers/DelphesPythia8.cpp
rf53a4d2 r28027d5 20 20 #include <iostream> 21 21 #include <sstream> 22 #include <string> 22 23 23 24 #include <signal.h> … … 38 39 #include "classes/DelphesClasses.h" 39 40 #include "classes/DelphesFactory.h" 41 #include "classes/DelphesLHEFReader.h" 40 42 41 43 #include "ExRootAnalysis/ExRootTreeWriter.h" … … 149 151 char appName[] = "DelphesPythia8"; 150 152 stringstream message; 153 FILE *inputFile = 0; 151 154 TFile *outputFile = 0; 152 155 TStopwatch readStopWatch, procStopWatch; 153 156 ExRootTreeWriter *treeWriter = 0; 154 157 ExRootTreeBranch *branchEvent = 0; 158 ExRootTreeBranch *branchEventLHEF = 0, *branchWeightLHEF = 0; 155 159 ExRootConfReader *confReader = 0; 156 160 Delphes *modularDelphes = 0; 157 161 DelphesFactory *factory = 0; 158 162 TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0; 163 TObjArray *stableParticleOutputArrayLHEF = 0, *allParticleOutputArrayLHEF = 0, *partonOutputArrayLHEF = 0; 164 DelphesLHEFReader *reader = 0; 159 165 Long64_t eventCounter, errorCounter; 160 166 Long64_t numberOfEvents, timesAllowErrors; … … 205 211 partonOutputArray = modularDelphes->ExportArray("partons"); 206 212 207 modularDelphes->InitTask(); 208 209 // Initialize pythia 213 // Initialize Pythia 210 214 pythia = new Pythia8::Pythia; 211 215 … … 221 225 numberOfEvents = pythia->mode("Main:numberOfEvents"); 222 226 timesAllowErrors = pythia->mode("Main:timesAllowErrors"); 227 228 inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r"); 229 if(inputFile) 230 { 231 reader = new DelphesLHEFReader; 232 reader->SetInputFile(inputFile); 233 234 branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class()); 235 branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class()); 236 237 allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF"); 238 stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF"); 239 partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF"); 240 } 241 242 modularDelphes->InitTask(); 223 243 224 244 pythia->init(); … … 234 254 for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter) 235 255 { 256 while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF, 257 stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady()); 258 236 259 if(!pythia->next()) 237 260 { 238 261 // If failure because reached end of file then exit event loop 239 if 262 if(pythia->info.atEndOfFile()) 240 263 { 241 264 cerr << "Aborted since reached end of Les Houches Event File" << endl; … … 244 267 245 268 // First few failures write off as "acceptable" errors, then quit 246 if (++errorCounter < timesAllowErrors) continue; 247 cerr << "Event generation aborted prematurely, owing to error!" << endl; 248 break; 269 if(++errorCounter > timesAllowErrors) 270 { 271 cerr << "Event generation aborted prematurely, owing to error!" << endl; 272 break; 273 } 274 275 modularDelphes->Clear(); 276 reader->Clear(); 277 continue; 249 278 } 250 279 … … 258 287 procStopWatch.Stop(); 259 288 289 if(reader) 290 { 291 reader->AnalyzeEvent(branchEventLHEF, eventCounter, &readStopWatch, &procStopWatch); 292 reader->AnalyzeWeight(branchWeightLHEF); 293 } 294 260 295 treeWriter->Fill(); 261 296 262 297 treeWriter->Clear(); 263 298 modularDelphes->Clear(); 299 if(reader) reader->Clear(); 264 300 265 301 readStopWatch.Start(); … … 277 313 cout << "** Exiting..." << endl; 278 314 315 delete reader; 279 316 delete pythia; 280 317 delete modularDelphes;
Note:
See TracChangeset
for help on using the changeset viewer.