Index: trunk/code/CMS_2013_I1256590_CONSTSUB.yoda =================================================================== --- trunk/code/CMS_2013_I1256590_CONSTSUB.yoda (revision 0) +++ trunk/code/CMS_2013_I1256590_CONSTSUB.yoda (revision 467) @@ -0,0 +1,206 @@ +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d01-x01-y03 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d01-x01-y03 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-02 2.500000e-02 2.500000e-02 1.02200e+00 4.60000e-02 4.60000e-02 +7.500000e-02 2.500000e-02 2.500000e-02 9.59000e-01 5.30000e-02 5.30000e-02 +1.250000e-01 2.500000e-02 2.500000e-02 8.62000e-01 6.40000e-02 6.50000e-02 +1.750000e-01 2.500000e-02 2.500000e-02 9.60000e-01 9.20000e-02 9.10000e-02 +2.250000e-01 2.500000e-02 2.500000e-02 1.20900e+00 1.37000e-01 1.38000e-01 +2.750000e-01 2.500000e-02 2.500000e-02 1.35700e+00 2.01000e-01 2.02000e-01 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d01-x01-y02 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d01-x01-y02 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-02 2.500000e-02 2.500000e-02 1.23630e+01 2.33000e-02 2.33000e-02 +7.500000e-02 2.500000e-02 2.500000e-02 4.75400e+00 1.34000e-2 1.34000e-02 +1.250000e-01 2.500000e-02 2.500000e-02 1.52800e+00 6.00000e-03 6.00000e-03 +1.750000e-01 2.500000e-02 2.500000e-02 7.36000e-01 4.00000e-03 4.00000e-03 +2.250000e-01 2.500000e-02 2.500000e-02 4.02000e-01 2.00000e-03 2.00000e-03 +2.750000e-01 2.500000e-02 2.500000e-02 2.17000e-01 2.00000e-03 2.00000e-03 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d01-x01-y01 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d01-x01-y01 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-02 2.500000e-02 2.500000e-02 1.26360e+01 9.84000e-01 9.83000e-01 +7.500000e-02 2.500000e-02 2.500000e-02 4.56000e+00 3.99000e-01 3.97000e-01 +1.250000e-01 2.500000e-02 2.500000e-02 1.31800e+00 1.41000e-01 1.40000e-01 +1.750000e-01 2.500000e-02 2.500000e-02 7.07000e-01 8.70000e-02 8.60000e-02 +2.250000e-01 2.500000e-02 2.500000e-02 4.86000e-01 6.50000e-02 6.60000e-02 +2.750000e-01 2.500000e-02 2.500000e-02 2.94000e-01 4.90000e-02 5.10000e-02 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d03-x01-y01 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d03-x01-y01 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-01 2.500000e-01 2.500000e-01 4.589493e-02 2.640026e-02 2.607364e-02 +7.500000e-01 2.500000e-01 2.500000e-01 3.302992e-01 1.372295e-01 1.308755e-01 +1.250000e+00 2.500000e-01 2.500000e-01 7.553553e-01 2.064520e-01 2.113501e-01 +1.750000e+00 2.500000e-01 2.500000e-01 1.219367e+00 1.948854e-01 1.904400e-01 +2.250000e+00 2.500000e-01 2.500000e-01 1.560549e+00 1.710532e-01 1.668595e-01 +2.750000e+00 2.500000e-01 2.500000e-01 1.804274e+00 2.437255e-01 2.817904e-01 +3.250000e+00 2.500000e-01 2.500000e-01 2.086064e+00 4.080724e-01 3.258002e-01 +3.750000e+00 2.500000e-01 2.500000e-01 2.447121e+00 4.787017e-01 4.655116e-01 +4.250000e+00 2.500000e-01 2.500000e-01 2.556013e+00 3.771229e-01 3.566198e-01 +4.750000e+00 2.500000e-01 2.500000e-01 1.273626e+00 4.731263e-01 4.790329e-01 +5.250000e+00 2.500000e-01 2.500000e-01 1.213483e-01 5.143852e-02 1.293548e-01 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d03-x01-y03 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d03-x01-y03 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-01 2.500000e-01 2.500000e-01 9.511278e-01 2.506266e-01 2.631579e-01 +7.500000e-01 2.500000e-01 2.500000e-01 9.949875e-01 1.817043e-01 1.879699e-01 +1.250000e+00 2.500000e-01 2.500000e-01 8.508772e-01 1.190476e-01 1.190476e-01 +1.750000e+00 2.500000e-01 2.500000e-01 8.132832e-01 1.002506e-01 9.398496e-02 +2.250000e+00 2.500000e-01 2.500000e-01 7.380952e-01 8.145363e-02 8.145363e-02 +2.750000e+00 2.500000e-01 2.500000e-01 7.506266e-01 8.145363e-02 8.145363e-02 +3.250000e+00 2.500000e-01 2.500000e-01 8.508772e-01 8.771930e-02 8.771930e-02 +3.750000e+00 2.500000e-01 2.500000e-01 1.095238e+00 1.190476e-01 1.127820e-01 +4.250000e+00 2.500000e-01 2.500000e-01 1.508772e+00 1.879699e-01 1.879699e-01 +4.750000e+00 2.500000e-01 2.500000e-01 1.765664e+00 3.320802e-01 3.320802e-01 +5.250000e+00 2.500000e-01 2.500000e-01 2.223058e+00 6.892231e-01 6.892231e-01 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d03-x01-y02 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d03-x01-y02 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-01 2.500000e-01 2.500000e-01 4.793716e-02 0.000000e+00 0.000000e+00 +7.500000e-01 2.500000e-01 2.500000e-01 3.302992e-01 0.000000e+00 0.000000e+00 +1.250000e+00 2.500000e-01 2.500000e-01 8.990452e-01 0.000000e+00 0.000000e+00 +1.750000e+00 2.500000e-01 2.500000e-01 1.494066e+00 0.000000e+00 0.000000e+00 +2.250000e+00 2.500000e-01 2.500000e-01 2.147498e+00 0.000000e+00 0.000000e+00 +2.750000e+00 2.500000e-01 2.500000e-01 2.447121e+00 0.000000e+00 0.000000e+00 +3.250000e+00 2.500000e-01 2.500000e-01 2.447121e+00 0.000000e+00 0.000000e+00 +3.750000e+00 2.500000e-01 2.500000e-01 2.243057e+00 0.000000e+00 0.000000e+00 +4.250000e+00 2.500000e-01 2.500000e-01 1.702521e+00 0.000000e+00 0.000000e+00 +4.750000e+00 2.500000e-01 2.500000e-01 7.231755e-01 0.000000e+00 0.000000e+00 +5.250000e+00 2.500000e-01 2.500000e-01 5.462544e-02 0.000000e+00 0.000000e+00 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d04-x01-y02 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d04-x01-y02 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +1.097304e+00 9.730405e-02 9.730405e-02 9.843526e-01 0.000000e+00 0.000000e+00 +1.310850e+00 1.162410e-01 1.162410e-01 1.032045e+00 0.000000e+00 0.000000e+00 +1.565954e+00 1.388630e-01 1.388630e-01 1.032045e+00 0.000000e+00 0.000000e+00 +1.870703e+00 1.658865e-01 1.658865e-01 9.689500e-01 0.000000e+00 0.000000e+00 +2.234760e+00 1.981700e-01 1.981700e-01 9.537884e-01 0.000000e+00 0.000000e+00 +2.669665e+00 2.367350e-01 2.367350e-01 8.407319e-01 0.000000e+00 0.000000e+00 +3.189208e+00 2.828065e-01 2.828065e-01 7.648242e-01 0.000000e+00 0.000000e+00 +3.809856e+00 3.378430e-01 3.378430e-01 6.132976e-01 0.000000e+00 0.000000e+00 +4.551290e+00 4.035905e-01 4.035905e-01 5.238152e-01 0.000000e+00 0.000000e+00 +5.437013e+00 4.821335e-01 4.821335e-01 4.473886e-01 0.000000e+00 0.000000e+00 +6.495107e+00 5.759600e-01 5.759600e-01 3.644549e-01 0.000000e+00 0.000000e+00 +7.759116e+00 6.880480e-01 6.880480e-01 2.922493e-01 0.000000e+00 0.000000e+00 +9.269112e+00 8.219480e-01 8.219480e-01 2.418587e-01 0.000000e+00 0.000000e+00 +1.107298e+01 9.819200e-01 9.819200e-01 1.764312e-01 0.000000e+00 0.000000e+00 +1.322790e+01 1.173000e+00 1.173000e+00 1.370839e-01 0.000000e+00 0.000000e+00 +1.580216e+01 1.401255e+00 1.401255e+00 1.015896e-01 0.000000e+00 0.000000e+00 +1.887738e+01 1.673970e+00 1.673970e+00 6.957701e-02 0.000000e+00 0.000000e+00 +2.255110e+01 1.999745e+00 1.999745e+00 4.917913e-02 0.000000e+00 0.000000e+00 +2.693975e+01 2.388910e+00 2.388910e+00 3.368195e-02 0.000000e+00 0.000000e+00 +3.218248e+01 2.853820e+00 2.853820e+00 2.306820e-02 0.000000e+00 0.000000e+00 +3.844549e+01 3.409190e+00 3.409190e+00 1.287031e-02 0.000000e+00 0.000000e+00 +4.592735e+01 4.072655e+00 4.072655e+00 6.532336e-03 0.000000e+00 0.000000e+00 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d04-x01-y03 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d04-x01-y03 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +1.097304e+00 9.730405e-02 9.730405e-02 1.050239e+00 5.598086e-01 5.645933e-01 +1.310850e+00 1.162410e-01 1.162410e-01 6.339713e-01 4.354067e-01 4.354067e-01 +1.565954e+00 1.388630e-01 1.388630e-01 6.818182e-01 3.636364e-01 3.684211e-01 +1.870703e+00 1.658865e-01 1.658865e-01 4.043062e-01 2.918660e-01 2.918660e-01 +2.234760e+00 1.981700e-01 1.981700e-01 2.224880e-01 2.344498e-01 2.392344e-01 +2.669665e+00 2.367350e-01 2.367350e-01 5.502392e-02 1.818182e-01 1.866029e-01 +3.189208e+00 2.828065e-01 2.828065e-01 -5.502392e-02 1.387560e-01 1.435407e-01 +3.809856e+00 3.378430e-01 3.378430e-01 -6.937799e-02 1.052632e-01 1.004785e-01 +4.551290e+00 4.035905e-01 4.035905e-01 -9.808612e-02 7.177033e-02 7.177033e-02 +5.437013e+00 4.821335e-01 4.821335e-01 -1.028708e-01 5.263158e-02 5.263158e-02 +6.495107e+00 5.759600e-01 5.759600e-01 -8.851675e-02 3.827751e-02 3.349282e-02 +7.759116e+00 6.880480e-01 6.880480e-01 -8.373206e-02 2.870813e-02 2.392344e-02 +9.269112e+00 8.219480e-01 8.219480e-01 -6.937799e-02 1.435407e-02 2.392344e-02 +1.107298e+01 9.819200e-01 9.819200e-01 -4.066986e-02 1.913876e-02 1.435407e-02 +1.322790e+01 1.173000e+00 1.173000e+00 -3.588517e-02 1.435407e-02 1.435407e-02 +1.580216e+01 1.401255e+00 1.401255e+00 -2.631579e-02 9.569378e-03 1.435407e-02 +1.887738e+01 1.673970e+00 1.673970e+00 -7.177033e-03 1.435407e-02 4.784689e-03 +2.255110e+01 1.999745e+00 1.999745e+00 -7.177033e-03 9.569378e-03 4.784689e-03 +2.693975e+01 2.388910e+00 2.388910e+00 -2.392344e-03 0.000000e+00 0.000000e+00 +3.218248e+01 2.853820e+00 2.853820e+00 2.392344e-03 0.000000e+00 0.000000e+00 +3.844549e+01 3.409190e+00 3.409190e+00 2.392344e-03 0.000000e+00 0.000000e+00 +4.592735e+01 4.072655e+00 4.072655e+00 -2.392344e-03 0.000000e+00 0.000000e+00 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d04-x01-y01 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d04-x01-y01 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +1.097304e+00 9.730405e-02 9.730405e-02 2.033383e+00 4.782024e-01 4.627072e-01 +1.310850e+00 1.162410e-01 1.162410e-01 1.656449e+00 3.489591e-01 3.451167e-01 +1.565954e+00 1.388630e-01 1.388630e-01 1.709530e+00 2.947629e-01 3.238530e-01 +1.870703e+00 1.658865e-01 1.658865e-01 1.370839e+00 2.541162e-01 2.341778e-01 +2.234760e+00 1.981700e-01 1.981700e-01 1.189440e+00 2.204898e-01 1.813987e-01 +2.669665e+00 2.367350e-01 2.367350e-01 8.954777e-01 1.774117e-01 1.696391e-01 +3.189208e+00 2.828065e-01 2.828065e-01 7.068302e-01 1.400369e-01 1.339018e-01 +3.809856e+00 3.378430e-01 3.378430e-01 5.491943e-01 1.088062e-01 1.040393e-01 +4.551290e+00 4.035905e-01 4.035905e-01 4.267141e-01 8.454043e-02 8.890473e-02 +5.437013e+00 4.821335e-01 4.821335e-01 3.421737e-01 6.779131e-02 5.845289e-02 +6.495107e+00 5.759600e-01 5.759600e-01 2.700890e-01 4.656965e-02 4.613880e-02 +7.759116e+00 6.880480e-01 6.880480e-01 2.065706e-01 3.561760e-02 3.528807e-02 +9.269112e+00 8.219480e-01 8.219480e-01 1.709530e-01 2.947629e-02 2.607165e-02 +1.107298e+01 9.819200e-01 9.819200e-01 1.328274e-01 2.115517e-02 2.269066e-02 +1.322790e+01 1.173000e+00 1.173000e+00 1.015896e-01 1.883195e-02 1.924512e-02 +1.580216e+01 1.401255e+00 1.401255e+00 7.648242e-02 1.611232e-02 1.593491e-02 +1.887738e+01 1.673970e+00 1.673970e+00 5.667932e-02 1.332960e-02 1.073732e-02 +2.255110e+01 1.999745e+00 1.999745e+00 4.006265e-02 1.083773e-02 1.149923e-02 +2.693975e+01 2.388910e+00 2.388910e+00 2.831749e-02 8.615028e-03 9.295896e-03 +3.218248e+01 2.853820e+00 2.853820e+00 1.909070e-02 6.620015e-03 6.266981e-03 +3.844549e+01 3.409190e+00 3.409190e+00 1.170828e-02 4.750580e-03 4.856209e-03 +4.592735e+01 4.072655e+00 4.072655e+00 6.636175e-03 2.933691e-03 3.207351e-03 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d02-x01-y01 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d02-x01-y01 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-02 2.500000e-02 2.500000e-02 6.375449e-01 0.000000e+00 0.000000e+00 +7.500000e-02 2.500000e-02 2.500000e-02 8.643001e-01 0.000000e+00 0.000000e+00 +1.250000e-01 2.500000e-02 2.500000e-02 9.278651e-01 0.000000e+00 0.000000e+00 +1.750000e-01 2.500000e-02 2.500000e-02 9.612714e-01 0.000000e+00 0.000000e+00 +2.250000e-01 2.500000e-02 2.500000e-02 9.858084e-01 0.000000e+00 0.000000e+00 +2.750000e-01 2.500000e-02 2.500000e-02 9.979320e-01 0.000000e+00 0.000000e+00 +# END YODA_SCATTER2D + +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_CONSTSUB/d02-x01-y02 +Path=/REF/CMS_2013_I1256590_CONSTSUB/d02-x01-y02 +Title= +Type=Scatter2D +# xval xerr- xerr+ yval yerr- yerr+ +2.500000e-02 2.500000e-02 2.500000e-02 6.357711e-01 0.000000e+00 0.000000e+00 +7.500000e-02 2.500000e-02 2.500000e-02 8.678478e-01 0.000000e+00 0.000000e+00 +1.250000e-01 2.500000e-02 2.500000e-02 9.366820e-01 0.000000e+00 0.000000e+00 +1.750000e-01 2.500000e-02 2.500000e-02 9.700882e-01 0.000000e+00 0.000000e+00 +2.250000e-01 2.500000e-02 2.500000e-02 9.892897e-01 0.000000e+00 0.000000e+00 +2.750000e-01 2.500000e-02 2.500000e-02 9.978656e-01 0.000000e+00 0.000000e+00 +# END YODA_SCATTER2D + Index: trunk/code/CMS_GROOMEDMASS_CONSTSUB.cc =================================================================== --- trunk/code/CMS_GROOMEDMASS_CONSTSUB.cc (revision 0) +++ trunk/code/CMS_GROOMEDMASS_CONSTSUB.cc (revision 467) @@ -0,0 +1,371 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/FastJets.hh" +//#include "Rivet/Tools/Logging.hh" +#include "Rivet/Tools/ParticleIdUtils.hh" +//#include "fastjet/AreaDefinition.hh" +#include "fastjet/ClusterSequence.hh" +//#include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/contrib/SoftDrop.hh" +#include "fastjet/contrib/Recluster.hh" +#include "HepMC/GenParticle.h" +#include "HepMC/GenEvent.h" + +namespace Rivet { + + using namespace fastjet; + + struct PJetComp { + bool operator() (const PseudoJet& pj1, const PseudoJet& pj2) const + {return pj1.pt() & scmap, + double & zcut, double & beta, double & zg, double & dR, int & Nstep, double & Mg, double & jetpt) { + PseudoJets subtrjet,pjets; + subtrjet = ConstSubJet(pj, scmap); + int Nsteptot(0); + PseudoJet pjet, parent1, parent2; + //pjet = contrib::Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(pj); + JetDefinition jet_def(cambridge_algorithm, JetDefinition::max_allowable_R); + ClusterSequence cs(subtrjet, jet_def); + pjets = sorted_by_pt(cs.exclusive_jets(1)); + pjet=pjets[0]; + jetpt = pjet.pt(); + while (pjet.has_parents(parent1,parent2)) { + //cout<<"unsubtracted pts: "< 1.) pt1 = 0.; + else pt1 = mom1.pT(); + if (mom2.E() < 0. || deltaR(mom2.eta(), mom2.phi(), Get4Mom(parent2).eta(), Get4Mom(parent2).phi()) > 1.) pt2 = 0.; + else pt2 = mom2.pT(); + //cout<<"subtracted pts: "< zcut*pow(dR/_jetR,beta)) { + double E1 = mom1.E(); + double E2 = mom2.E(); + double theta = angle(mom1.p3(), mom2.p3()); + Mg = sqrt(2.*E1*E2*(1.-cos(theta))); + return pjet; + } + if (pt1 > pt2) { + pjet = parent1; + } + else { + pjet = parent2; + } + Nsteptot += 1; + if (zg > 0.) Nstep += 1; + } + //cout<<"mySoftDrop failed.\n"; + return PseudoJet(); + } + + + /// Constituent subtraction + PseudoJets ConstSubJet(PseudoJet pjet, map & scmap) { + // sort constituents into two vectors: particles and ghosts (subtraction momenta) + vector parts; + vector ghosts; + foreach (PseudoJet pj, pjet.constituents()) { + if (fuzzyEquals(pj.E(),1e-6,1e-4)) { + map::iterator mapit; + mapit = scmap.find(pj); + if (mapit == scmap.end()) cout<<"Error: did not find matching scattering centre in map.\n"; + else { + FourMomentum submom = Get4Mom(mapit->second); + My4Mom ghost; + ghost.pt = submom.pT(); + ghost.mdelta = sqrt(submom.mass2() + sqr(submom.pT())) - submom.pT(); + ghost.phi = submom.p3().phi(); + ghost.y = submom.rapidity(); + ghosts.push_back(ghost); + } + } + else { + FourMomentum partmom = Get4Mom(pj); + My4Mom part; + part.pt = partmom.pT(); + part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT(); + part.phi = partmom.p3().phi(); + part.y = partmom.rapidity(); + parts.push_back(part); + } + } + //cout<<"number of particles: "< dists; + for (size_t i=0; i < parts.size(); ++i) { + for (size_t j=0; j < ghosts.size(); ++j) { + //cout<<"i, j = "<::iterator liter=dists.begin(); liter != dists.end(); ++liter) { + //cout<<"dealing with dist "<dR<ipart; + size_t gnum = liter->ighost; + double ptp = parts[pnum].pt; + double ptg = ghosts[gnum].pt; + //cout<<"pts: "< ptg) { + parts[pnum].pt -= ptg; + ghosts[gnum].pt = 0.; + } + else { + ghosts[gnum].pt -= ptp; + parts[pnum].pt = 0.; + } + double mdp = parts[pnum].mdelta; + double mdg = ghosts[gnum].mdelta; + //cout<<"masses: "< mdg) { + parts[pnum].mdelta -= mdg; + ghosts[gnum].mdelta = 0.; + } + else { + ghosts[gnum].mdelta -= mdp; + parts[pnum].mdelta = 0.; + } + } + + //sum up resulting 4-momenta to get subtracted jet momentum + FourMomentum partmom, ghostmom, jetmom; + PseudoJets result; + for (size_t i=0; i < parts.size(); ++i) { + if (parts[i].pt > 0.) { + FourMomentum mom((parts[i].pt+parts[i].mdelta)*cosh(parts[i].y), + parts[i].pt*cos(parts[i].phi), + parts[i].pt*sin(parts[i].phi), + (parts[i].pt+parts[i].mdelta)*sinh(parts[i].y)); + partmom += mom; + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + result.push_back(pjet); + } + } + for (size_t i=0; i < ghosts.size(); ++i) { + if (ghosts[i].pt > 0.) { + FourMomentum mom((ghosts[i].pt+ghosts[i].mdelta)*cosh(ghosts[i].y), + ghosts[i].pt*cos(ghosts[i].phi), + ghosts[i].pt*sin(ghosts[i].phi), + (ghosts[i].pt+ghosts[i].mdelta)*sinh(ghosts[i].y)); + ghostmom += mom; + //cout<<"ghost rapidity: "< ptsubs; +//cout<<"==============================\n"; + PseudoJets scatcens; + + vector scafter; + foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { + FourMomentum mom(p->momentum()); + if (fabs(mom.eta()) < 10.) { + if (p->status() == 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + scatcens.push_back(pjet); + } + if(p->status() != 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + pevent.push_back(pjet); + } + } + } + + map scmap; + size_t j(0); + for (size_t i=0; i= _ptedges[l] && ptsub < _ptedges[l+1]) ptbin=l; + } + + int ptbin2(-1); + for (size_t l = 0; l < _nptbins2; ++l) { + if (ptsub >= _ptedges2[l] && ptsub < _ptedges2[l+1]) ptbin2=l; + } + + if (sd_pjet_0 != 0) { + if (ptbin2 >= 0) _h_dR[ptbin2]->fill(dR0,weight); + if (ptbin2 >= 0) _h_nstep[ptbin2]->fill(Nstep0,weight); + if (dR0 > 0.1) { + if (ptbin2 >= 0) _h_zg[ptbin2]->fill(zg0,weight); + if (ptbin >= 0) _h_Mgopt_1[ptbin]->fill(Mg0/ptsub,weight); + } + } + if (sd_pjet_1 != 0) { + if (dR1 > 0.1) { + if (ptbin >= 0) _h_Mgopt_2[ptbin]->fill(Mg1/ptsub,weight); + } + } + } + } + + } + + + /// Normalise histograms etc., after the run + void finalize() { + + for (size_t k = 0; k < _nptbins; ++k) { + normalize(_h_Mgopt_1[k]); + normalize(_h_Mgopt_2[k]); + } + for (size_t k = 0; k < _nptbins2; ++k) { + normalize(_h_zg[k]); + normalize(_h_dR[k]); + normalize(_h_nstep[k]); + } + + } + + + + private: + + vector _ptedges, _ptedges2, _zcuts, _betas; + size_t _nptbins, _nptbins2; + double _jetR; + + Histo1DPtr _h_Mgopt_1[4]; + Histo1DPtr _h_Mgopt_2[4]; + Histo1DPtr _h_zg[6]; + Histo1DPtr _h_dR[6]; + Histo1DPtr _h_nstep[6]; + + + }; + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(CMS_GROOMEDMASS_CONSTSUB); +} Index: trunk/code/CMS_2014_I1299142.plot =================================================================== --- trunk/code/CMS_2014_I1299142.plot (revision 466) +++ trunk/code/CMS_2014_I1299142.plot (revision 467) @@ -1,8 +1,171 @@ # BEGIN PLOT /CMS_2014_I1299142/d01-x01-y01 -#Title=[Uncomment and insert title for histogram d01-x01-y01 here] -#XLabel=[Uncomment and insert x-axis label for histogram d01-x01-y01 here] +Title=p+p ; $100\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$\xi$ #YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] # + any additional plot settings you might like, see make-plots documentation # END PLOT -# ... add more histograms as you need them ... +# BEGIN PLOT /CMS_2014_I1299142/d01-x02-y01 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$\xi$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d01-x03-y01 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$\xi$ +YLabel=PbPb / pp +LogY=0 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d02-x01-y01 +Title=p+p ; $100\,\text{GeV} < p_\perp^\text{jet} < 120\,\text{GeV}$ +XLabel=$\xi$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d02-x02-y01 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 120\,\text{GeV}$ +XLabel=$\xi$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d02-x03-y01 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 120\,\text{GeV}$ +XLabel=$\xi$ +YLabel=PbPb / pp +LogY=0 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d03-x01-y01 +Title=p+p ; $120\,\text{GeV} < p_\perp^\text{jet} < 150\,\text{GeV}$ +XLabel=$\xi$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d03-x02-y01 +Title=Pb+Pb 0-10\% ; $120\,\text{GeV} < p_\perp^\text{jet} < 150\,\text{GeV}$ +XLabel=$\xi$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d03-x03-y01 +Title=Pb+Pb 0-10\% ; $120\,\text{GeV} < p_\perp^\text{jet} < 150\,\text{GeV}$ +XLabel=$\xi$ +YLabel=PbPb / pp +LogY=0 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d04-x01-y01 +Title=p+p ; $150\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$\xi$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d04-x02-y01 +Title=Pb+Pb 0-10\% ; $150\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$\xi$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d04-x03-y01 +Title=Pb+Pb 0-10\% ; $150\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$\xi$ +YLabel=PbPb / pp +LogY=0 +# END PLOT + + +# BEGIN PLOT /CMS_2014_I1299142/d01-x01-y02 +Title=p+p ; $100\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +LogX=1 +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d01-x02-y02 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +LogX=1 +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d01-x03-y02 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +YLabel=PbPb - pp +LogX=1 +LogY=0 +RatioPlotMode=deviation +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d02-x01-y02 +Title=p+p ; $100\,\text{GeV} < p_\perp^\text{jet} < 120\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +LogX=1 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d02-x02-y02 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 120\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +LogX=1 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d02-x03-y02 +Title=Pb+Pb 0-10\% ; $100\,\text{GeV} < p_\perp^\text{jet} < 120\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +YLabel=PbPb - pp +LogY=0 +LogX=1 +RatioPlotMode=deviation +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d03-x01-y02 +Title=p+p ; $120\,\text{GeV} < p_\perp^\text{jet} < 150\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +LogX=1 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d03-x02-y02 +Title=Pb+Pb 0-10\% ; $120\,\text{GeV} < p_\perp^\text{jet} < 150\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +LogX=1 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d03-x03-y02 +Title=Pb+Pb 0-10\% ; $120\,\text{GeV} < p_\perp^\text{jet} < 150\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +YLabel=PbPb - pp +LogY=0 +LogX=1 +RatioPlotMode=deviation +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d04-x01-y02 +Title=p+p ; $150\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +LogX=1 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d04-x02-y02 +Title=Pb+Pb 0-10\% ; $150\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here] +LogX=1 +# END PLOT + +# BEGIN PLOT /CMS_2014_I1299142/d04-x03-y02 +Title=Pb+Pb 0-10\% ; $150\,\text{GeV} < p_\perp^\text{jet} < 300\,\text{GeV}$ +XLabel=$p_\perp\ [\text{GeV}]$ +YLabel=PbPb - pp +LogY=0 +LogX=1 +RatioPlotMode=deviation +# END PLOT + + Index: trunk/code/CMS_2014_I1299142.cc =================================================================== --- trunk/code/CMS_2014_I1299142.cc (revision 466) +++ trunk/code/CMS_2014_I1299142.cc (revision 467) @@ -1,317 +1,602 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { using namespace fastjet; struct FourMomComp { bool operator() (const FourMomentum& p1, const FourMomentum& p2) const {return p1.pT() 0.5) n=0.001; if (cent > 0.3 && cent < 0.5) n=3.88; if (cent > 0.1 && cent < 0.3) n=5.10; if (cent > 0. && cent < 0.1) n=5.23; n=5.23; double sigma(sqrt(c*c+s*s/pttrue+n*n/(pttrue*pttrue))); double r1(1.0*rand()/RAND_MAX), r2(1.0*rand()/RAND_MAX); double fac(max(sqrt(-2.*log(r1))*cos(2.*M_PI*r2)*sigma+1.,0.)); return fac*pttrue; } /// extract thermal momenta, that should be subtracted, from HepMC event vector extractThermalMomenta(const Event & event) { vector thermom; foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { FourMomentum mom(p->momentum()); if (p->status() == 3) { thermom.push_back(mom); } } return thermom; } /// build map from dummy particles to thermal momenta to speed up subtraction map buildThMomMap(const Event & event, const vector * thermom) { map thmap; + size_t j(0); foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { FourMomentum mom(p->momentum()); - if (mom.E() < 1e-5 && mom.E() > 1e-7) { - FourMomentum dummy(mom); - dummy *= 10000.; - double dRmin(1.); - FourMomentum match; - foreach (FourMomentum tm, *thermom) { - double dR(deltaR(dummy, tm)); - if (dR < dRmin) { - dRmin = dR; - match = tm; - } - } - if (dRmin < 1e-5) { - map::iterator mapit; - mapit = thmap.find(mom); - if (mapit != thmap.end()) { - cout<<"Error: dummy is already in map.\n"; - cout<first<second<size()) { + cout<<"Error: number of dummies does not match number of subtractions, will veto Event"< * tmmap, double fac = 1.) { + FourMomentum SubtractJetMom(Jet jet, map * tmmap) { if (tmmap->empty()) return jet.momentum(); FourMomentum sub(0.,0.,0.,0.); foreach (Particle part, jet.constituents()) { - if (part.E() < 1e-5 && part.E() > 1e-7) { + if (fuzzyEquals(part.E(), 1e-6)) { map::iterator mapit; mapit = (*tmmap).find(part.momentum()); if (mapit == (*tmmap).end()) cout<<"Error: did not find matching scattering centre in map.\n"<second; } } FourMomentum jetmom(jet.momentum()); //cout<<"Original momentum: "< > fillGrid(const Event & event, std::set > & gridevent) { // initialise grid vector > grid; grid.resize(_nphibins); for (size_t i = 0; i < _nphibins; ++i) { grid[i].resize(_netabins); for (size_t k = 0; k < _nphibins; ++k) { grid[i][k] = FourMomentum(0., 0., 0., 0.); } } // fill grid foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { FourMomentum mom(p->momentum()); if (fabs(mom.eta()) < 4.) { size_t phibin, etabin; phibin = int(mom.phi(ZERO_2PI)*_nphibins/(2.*M_PI)); etabin = int((mom.eta() + 4.)*_netabins/8.); if (phibin < 0 || phibin > _nphibins-1) {std::cout<<"Error: "< "< _netabins-1) {std::cout<<"Error: "< "<status() == 3) { grid[phibin][etabin].setE(grid[phibin][etabin].E() - mom.E()); grid[phibin][etabin].setPx(grid[phibin][etabin].px() - mom.px()); grid[phibin][etabin].setPy(grid[phibin][etabin].py() - mom.py()); grid[phibin][etabin].setPz(grid[phibin][etabin].pz() - mom.pz()); } else { grid[phibin][etabin].setE(grid[phibin][etabin].E() + mom.E()); grid[phibin][etabin].setPx(grid[phibin][etabin].px() + mom.px()); grid[phibin][etabin].setPy(grid[phibin][etabin].py() + mom.py()); grid[phibin][etabin].setPz(grid[phibin][etabin].pz() + mom.pz()); } /*double cellphi((phibin+0.5)*2*M_PI/_nphibins), celleta((etabin+0.5)*8./_netabins-4.); double celltheta(2.*atan(exp(-celleta))); if (p->status() == 3) { double rho(grid[phibin][etabin].E() - mom.E()); grid[phibin][etabin].setE(rho); grid[phibin][etabin].setPx(fabs(rho)*sin(celltheta)*cos(cellphi)); grid[phibin][etabin].setPy(fabs(rho)*sin(celltheta)*sin(cellphi)); grid[phibin][etabin].setPz(fabs(rho)*cos(celltheta)); } else { // use this for standard recoil treatment double rho(grid[phibin][etabin].E() + mom.E()); grid[phibin][etabin].setE(rho); grid[phibin][etabin].setPx(rho*sin(celltheta)*cos(cellphi)); grid[phibin][etabin].setPy(rho*sin(celltheta)*sin(cellphi)); grid[phibin][etabin].setPz(rho*cos(celltheta)); }*/ gridevent.insert(std::pair(phibin,etabin)); } } return grid; } + /// event-wise Constituent subtraction based on Particles + Particles ConstSubPartEvent(const Particles particles, PseudoJets pjghosts) { + // sort constituents into two vectors: particles and ghosts (subtraction momenta) + Particles subevent; + vector parts; + vector ghosts; + + foreach (Particle p, particles) { + FourMomentum partmom = p.momentum(); + if (fuzzyEquals(partmom.E(),1e-6)) continue; + if (partmom.E() - fabs(partmom.pz()) > 0.) { + MyPart part; + part.id = p.pid(); + part.pt = partmom.pT(); + part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT(); + part.phi = partmom.p3().phi(); + part.y = partmom.rapidity(); + //if (abs(part.y) > 5.) { + // cout<<"particle rapidity: "< 0.) { + My4Mom ghost; + ghost.pt = submom.pT(); + ghost.mdelta = sqrt(submom.mass2() + sqr(submom.pT())) - submom.pT(); + ghost.phi = submom.p3().phi(); + ghost.y = submom.rapidity(); + //cout<<"ghost rapidity: "< dists; + for (size_t i=0; i < parts.size(); ++i) { + for (size_t j=0; j < ghosts.size(); ++j) { + //cout<<"i, j = "< M_PI) Deltaphi = 2.*M_PI - Deltaphi; + dist.dR = sqrt(sqr(Deltaphi) + sqr(parts[i].y-ghosts[j].y)); + dist.ipart = i; + dist.ighost = j; + //cout<<"distance: "<::iterator liter=dists.begin(); liter != dists.end(); ++liter) { + //cout<<"dealing with dist "<dR<dR > _dRloose) break; + size_t pnum = liter->ipart; + size_t gnum = liter->ighost; + double ptp = parts[pnum].pt; + double ptg = ghosts[gnum].pt; + //cout<<"pts: "< ptg) { + parts[pnum].pt -= ptg; + ghosts[gnum].pt = 0.; + } + else { + ghosts[gnum].pt -= ptp; + parts[pnum].pt = 0.; + } + double mdp = parts[pnum].mdelta; + double mdg = ghosts[gnum].mdelta; + //cout<<"masses: "< mdg) { + parts[pnum].mdelta -= mdg; + ghosts[gnum].mdelta = 0.; + } + else { + ghosts[gnum].mdelta -= mdp; + parts[pnum].mdelta = 0.; + } + + } + + //sum up resulting 4-momenta to get subtracted jet momentum + int nparts(0); + for (size_t i=0; i < parts.size(); ++i) { + if (parts[i].pt > 0.) { + FourMomentum mom((parts[i].pt+parts[i].mdelta)*cosh(parts[i].y), + parts[i].pt*cos(parts[i].phi), + parts[i].pt*sin(parts[i].phi), + (parts[i].pt+parts[i].mdelta)*sinh(parts[i].y)); + //cout<<"particle rapidity: "< 0.) { + FourMomentum mom((ghosts[i].pt+ghosts[i].mdelta)*cosh(ghosts[i].y), + ghosts[i].pt*cos(ghosts[i].phi), + ghosts[i].pt*sin(ghosts[i].phi), + (ghosts[i].pt+ghosts[i].mdelta)*sinh(ghosts[i].y)); + //cout<<"ghost rapidity: "<heavy_ion()?event.genEvent()->heavy_ion()->impact_parameter():-1.); if (cent > 0.1) vetoEvent; const ChargedFinalState cfs = applyProjection(event, "CFS"); - + Particles tracks; + Jets jetssub; #if MODE==0 + tracks = cfs.particles(); const Jets alljets = applyProjection(event, "Jets").jetsByPt(3.*GeV); vector thermom = extractThermalMomenta(event); map thmommap = buildThMomMap(event, &thermom); foreach (Jet jet, alljets) { - FourMomentum jetsub = SubtractJetMom(jet, &thmommap, 1.); + FourMomentum jetsub = SubtractJetMom(jet, &thmommap); jetssub.push_back(Jet(jetsub)); } #elif MODE==1 + tracks = cfs.particles(); PseudoJets pevent, pjets; std::set > gridevent; vector > grid = fillGrid(event, gridevent); for (set >::iterator it=gridevent.begin(); it!=gridevent.end(); ++it) { size_t i,k; i = it->first; k = it->second; if (grid[i][k].E() > 0. ) { PseudoJet part(grid[i][k].px(),grid[i][k].py(),grid[i][k].pz(),grid[i][k].E()); pevent.push_back(part); } } JetDefinition jet_def(antikt_algorithm, _jetR); ClusterSequence cs(pevent, jet_def); pjets = sorted_by_pt(cs.inclusive_jets(3.*GeV)); foreach (PseudoJet pjet, pjets) { jetssub.push_back(Jet(pjet)); } +#elif MODE==2 + const FinalState fs = applyProjection(event, "FS"); + const Particles origevent = fs.particles(); + + PseudoJets ghosts; + + foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { + FourMomentum mom(p->momentum()); + if (fabs(mom.eta()) < _etaloose) { + if (p->status() == 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + ghosts.push_back(pjet); + } + } + } + + Particles subevent = ConstSubPartEvent(origevent, ghosts); + + PseudoJets newevent; + foreach (Particle part, subevent) { + newevent.push_back(PseudoJet(part.px(), part.py(), part.pz(), part.E())); + if (isCharged(part) && part.pT() > _trackptmin) tracks.push_back(part); + } + + JetDefinition jet_def(antikt_algorithm, _jetR); + ClusterSequence cs(newevent, jet_def); + PseudoJets pjets; + pjets = cs.inclusive_jets(10.); + + foreach (PseudoJet pj, pjets) { + jetssub.push_back(Jet(pj)); + } + +#elif MODE==3 + const FinalState fs = applyProjection(event, "FS"); + const Jets loosejets = applyProjection(event, "Jets").jetsByPt(_ptminloose); + + Particles origevent; + PseudoJets ghosts; + + foreach (Particle part, fs.particles()) { + foreach (Jet jet, loosejets) { + if (deltaR(part.momentum(),jet.momentum()) < _dRloose) { + origevent.push_back(part); + break; + } + } + } + + foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { + foreach (Jet jet, loosejets) { + FourMomentum mom(p->momentum()); + if (deltaR(mom,jet.momentum()) < _dRloose && p->status() == 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + ghosts.push_back(pjet); + break; + } + } + } + + Particles subevent = ConstSubPartEvent(origevent, ghosts); + + PseudoJets newevent; + foreach (Particle part, subevent) { + newevent.push_back(PseudoJet(part.px(), part.py(), part.pz(), part.E())); + if (isCharged(part) && (part.pT() > _trackptmin)) tracks.push_back(part); + } + + JetDefinition jet_def(antikt_algorithm, _jetR); + ClusterSequence cs(newevent, jet_def); + PseudoJets pjets; + pjets = cs.inclusive_jets(_ptminloose); + + foreach (PseudoJet pj, pjets) { + jetssub.push_back(Jet(pj)); + } + #endif map jetmap; foreach (Jet jet, jetssub) { - double ptreco, ptgen; - ptgen = jet.momentum().pT(); - ptreco = ptsmear(ptgen,cent); - jetmap[jet] = ptreco; + double ptreco, ptgen; + ptgen = jet.momentum().pT(); + ptreco = ptsmear(ptgen,cent); + jetmap[jet] = ptreco; } Jets jets; foreach (Jet jet, jetssub) { if (fabs(jet.momentum().eta()) > _etamin && fabs(jet.momentum().eta()) < _etamax && - jetmap[jet] > _ptmin && jetmap[jet] > _ptmax) + jetmap[jet] > _ptmin && jetmap[jet] < _ptmax) jets.push_back(jet); } foreach (Jet jet, jets) { - _Njets+=weight; + _Njets[0]+=weight; double jetpt = jetmap[jet]; - foreach (Particle part, cfs.particles()){ + int ptbin(-1); + for (size_t l = 0; l < _Nptbins; ++l) { + if (jetpt >= _ptedges[l] && jetpt < _ptedges[l+1]) ptbin=l; + } + //shift: first histo is inclusive + if (ptbin >= 0) ptbin+=1; + _Njets[ptbin]+=weight; + foreach (Particle part, tracks){ double dR(deltaR(part,jet)); if (dR >= _jetR) continue; double z,xi; - z = part.momentum().p3().dot(jet.momentum().p3())/jet.momentum().p3().mod()/jet.momentum().p3().mod(); - z*=jet.momentum().pT()/jetpt; + z = part.momentum().p3().dot(jet.momentum().p3())/sqr(jet.momentum().p3().mod()); + //z*=jet.momentum().pT()/jetpt; xi = log(1./z); - _h_xi->fill(xi, weight); + if (cent < 0.) { + if (ptbin > 0) { + _h_xi_pp[ptbin]->fill(xi, weight); + _h_xi_pp[0]->fill(xi, weight); + _h_pt_pp[ptbin]->fill(part.pT(), weight); + _h_pt_pp[0]->fill(part.pT(), weight); + } + } + else { + if (ptbin > 0) { + _h_xi_PbPb[ptbin]->fill(xi, weight); + _h_xi_PbPb[0]->fill(xi, weight); + _h_pt_PbPb[ptbin]->fill(part.pT(), weight); + _h_pt_PbPb[0]->fill(part.pT(), weight); + } + } + _h_xi_ratio[ptbin]->fill(xi, weight); + _h_xi_ratio[0]->fill(xi, weight); + _h_pt_diff[ptbin]->fill(part.pT(), weight); + _h_pt_diff[0]->fill(part.pT(), weight); } } } /// Normalise histograms etc., after the run void finalize() { - if (_Njets==0.) _Njets=1.; - scale(_h_xi,1./_Njets); + for (size_t i = 0; i < _Nptbins+1; ++i) { + if (_Njets[i]==0.) _Njets[i]=1.; + scale(_h_xi_pp[i],1./_Njets[i]); + scale(_h_xi_PbPb[i],1./_Njets[i]); + scale(_h_xi_ratio[i],1./_Njets[i]); + scale(_h_pt_pp[i],1./_Njets[i]); + scale(_h_pt_PbPb[i],1./_Njets[i]); + scale(_h_pt_diff[i],1./_Njets[i]); + } } private: - double _Njets; + vector _Njets, _ptedges; double _jetR, _trackptmin, _etamin, _etamax, _ptmin, _ptmax; + double _etaloose, _ptminloose, _dRloose; size_t _netabins, _nphibins; + size_t _Nptbins; - Histo1DPtr _h_xi; + Histo1DPtr _h_xi_pp[4], _h_xi_PbPb[4], _h_xi_ratio[4]; + Histo1DPtr _h_pt_pp[4], _h_pt_PbPb[4], _h_pt_diff[4]; }; // The hook for the plugin system -#if MODE==0 +/*#if MODE==0 DECLARE_RIVET_PLUGIN(CMS_2014_I1299142_4MOMSUB); #elif MODE==1 DECLARE_RIVET_PLUGIN(CMS_2014_I1299142_GRIDSUB2); -#endif +#elif MODE==3 + DECLARE_RIVET_PLUGIN(CMS_2014_I1299142_CONSTSUB); +#endif*/ + DECLARE_RIVET_PLUGIN(CMS_2014_I1299142); } Index: trunk/code/CMS_SOFTDROP_4MOMSUB.cc =================================================================== --- trunk/code/CMS_SOFTDROP_4MOMSUB.cc (revision 0) +++ trunk/code/CMS_SOFTDROP_4MOMSUB.cc (revision 467) @@ -0,0 +1,570 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/FastJets.hh" +//#include "Rivet/Tools/Logging.hh" +#include "Rivet/Tools/ParticleIdUtils.hh" +//#include "fastjet/AreaDefinition.hh" +#include "fastjet/ClusterSequence.hh" +//#include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/contrib/SoftDrop.hh" +#include "fastjet/contrib/Recluster.hh" +#include "HepMC/GenParticle.h" +#include "HepMC/GenEvent.h" + +#ifndef MODE +#define MODE 0 +#endif + +namespace Rivet { + + using namespace fastjet; + + struct PJetComp { + bool operator() (const PseudoJet& pj1, const PseudoJet& pj2) const + {return pj1.pt() * scpoint, map & scmap, + double & zg, double & dR, int & Nstep, double & gsj) { + int Nsteptot(0); + PseudoJet pjet, parent1, parent2; + pjet = contrib::Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(pj); + while (pjet.has_parents(parent1,parent2)) { + //cout<<"unsubtracted pts: "< 1.) pt1 = 0.; + else pt1 = mom1.pT(); + if (mom2.E() < 0. || deltaR(mom2.eta(), mom2.phi(), Get4Mom(parent2).eta(), Get4Mom(parent2).phi()) > 1.) pt2 = 0.; + else pt2 = mom2.pT(); + //cout<<"subtracted pts: "< _zcut*pow(dR,_beta)) { + double deltar; + FourMomentum momsum(mom1+mom2); + deltar = deltaR(momsum.eta(),momsum.phi(),mom1.eta(), mom1.phi()); + gsj = mom1.pT()*deltar; + deltar = deltaR(momsum.eta(),momsum.phi(),mom2.eta(), mom2.phi()); + gsj += mom2.pT()*deltar; + gsj /= momsum.pt(); + return pjet; + } + if (pt1 > pt2) { + pjet = parent1; + } + else { + pjet = parent2; + } + Nsteptot += 1; + if (zg > 0.) Nstep += 1; + } + //cout<<"mySoftDrop failed.\n"; + return PseudoJet(); + } + + /// 4-momentum subtraction + FourMomentum SubtractJetMom(PseudoJet pjet, const vector * scatcenters, map & scmap) { + if (scatcenters->empty()) return Get4Mom(pjet); + FourMomentum sub(0.,0.,0.,0.); + bool fillmap(scmap.empty()); + //cout<<"--------------------------------------------------\n"; + //cout<<"size of map: "<::iterator mapit; + mapit = scmap.find(pj); + /*if (mapit != scmap.end()) { + cout<<"Error: dummy is already in map.\n"; + cout<first)<second)<::iterator mapit; + mapit = scmap.find(pj); + if (mapit == scmap.end()) cout<<"Error: did not find matchings scattering centre in map.\n"; + else sub += Get4Mom(mapit->second); + } + } + } + //cout<<"size of map: "< * scatcenters, map & scmap) { + double firstmom(0.); + foreach (PseudoJet part, pjet.constituents()) { + if (part.E() < 1e-5) { + map::iterator mapit; + mapit = scmap.find(part); + if (mapit != scmap.end()) { + double dR(deltaR(jetmom.eta(), jetmom.phi(), Get4Mom(mapit->second).eta(), Get4Mom(mapit->second).phi())); + firstmom -= (mapit->second).pt()*dR; + } + } + else { + double dR(deltaR(jetmom.eta(), jetmom.phi(), part.eta(), part.phi())); + firstmom += part.pt()*dR; + } + } + return firstmom/jetmom.pT(); + } + + + /// recoil energy fraction + double RecEnFrac(PseudoJet pjet, FourMomentum jetmom, vector scafter, map & scmap){ + double recE(0.), subE(0.); + foreach (PseudoJet pj, pjet.constituents()) { + if (pj.E() < 1e-5) { + map::iterator mapit; + mapit = scmap.find(pj); + if (mapit != scmap.end()) { + subE += (mapit->second).E(); + } + } + else { + double dRmin(10.), Emin(0.); + foreach (FourMomentum fmom, scafter) { + double dR(deltaR(Get4Mom(pj).eta(), Get4Mom(pj).phi(), fmom.eta(), fmom.phi())); + if (dR < dRmin) { + dRmin = dR; + Emin = fmom.E(); + } + } + if (dRmin < 1e-5) { + recE += Emin; + } + } + } + return (recE-subE)/jetmom.E(); + } + + + double RecPtFrac(PseudoJet pjet, FourMomentum jetmom, vector scafter, map & scmap){ + FourMomentum recmom(0.,0.,0.,0.), submom(0.,0.,0.,0.); + foreach (PseudoJet pj, pjet.constituents()) { + if (pj.E() < 1e-5) { + map::iterator mapit; + mapit = scmap.find(pj); + if (mapit != scmap.end()) { + submom += Get4Mom(mapit->second); + } + } + else { + double dRmin(10.); + FourMomentum momatmin; + foreach (FourMomentum fmom, scafter) { + double dR(deltaR(Get4Mom(pj).eta(), Get4Mom(pj).phi(), fmom.eta(), fmom.phi())); + if (dR < dRmin) { + dRmin = dR; + momatmin = fmom; + } + } + if (dRmin < 1e-5) { + recmom += momatmin; + } + } + } + recmom -= submom; + double ptfrac(recmom.pT()/jetmom.pt()); + double dR(deltaR(recmom.eta(), recmom.phi(), jetmom.eta(), jetmom.phi())); + if (dR > 1.) ptfrac *= -1.; + return ptfrac; + } + + + + /// 4-momentum adapter + FourMomentum Get4Mom(PseudoJet pjet){ + FourMomentum mom(pjet.E(), pjet.px(), pjet.py(), pjet.pz()); + return mom; + } + + /// Perform the per-event analysis + void analyze(const Event& event) { + const double weight = event.weight(); + + PseudoJets pevent, pjets; + +// PseudoJets evnorec; + list ptsubs; +//cout<<"==============================\n"; + PseudoJets scatcens; + vector scafter; + foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { + FourMomentum mom(p->momentum()); + if (fabs(mom.eta()) < 4.) { + if (p->status() == 3) { +#if MODE==0 + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + scatcens.push_back(pjet); +#endif + _h_pt_sc_bef->fill(mom.pT(),weight); + } +#if MODE==0 + if(p->status() != 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + pevent.push_back(pjet); + } +#else + if(p->status() == 1) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + pevent.push_back(pjet); + } +#endif + if (p->status() == 4) { + _h_pt_sc_aft->fill(mom.pT(),weight); +#if MODE==0 + scafter.push_back(mom); +#endif + } + /*if (p->status() == 1) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + evnorec.push_back(pjet); + }*/ + } + } + + for (size_t k = 0; k < _njetRbins; ++k) { + + JetDefinition jet_def(antikt_algorithm, _jetRs[k]); + ClusterSequence cs(pevent, jet_def); + pjets = sorted_by_pt(cs.inclusive_jets(50.0)); + + foreach (PseudoJet pjet, pjets) { + if (fabs(pjet.eta()) < 2.5) { + + map scmap; + FourMomentum jetsub = SubtractJetMom(pjet, &scatcens, scmap); + double ptsub; + if (jetsub.E() < 0. || deltaR(jetsub.eta(), jetsub.phi(), Get4Mom(pjet).eta(), Get4Mom(pjet).phi()) > 1.) ptsub = 0.; + else ptsub = jetsub.pT(); + + if (ptsub < 50.) continue; + _h_jetpt[k]->fill(ptsub, weight); + //ptsubs.push_back(ptsub); + + int ptbin(-1); + for (size_t l = 0; l < _nptbins; ++l) { + if (ptsub >= _ptedges[l] && ptsub < _ptedges[l+1]) ptbin=l; + } + if (ptbin < 0) continue; + + double g; + g = Girth(pjet, jetsub, &scatcens, scmap); + _h_girth[k]->fill(g,weight); + + int gbin(-1); + for (size_t j = 0; j < _ngirthbins; ++j) { + if (g >= _girthedges[j]*_jetRs[k] && g < _girthedges[j+1]*_jetRs[k]) gbin=j; + } + + double zg(0.), dR(0.), gsj(0.); + int Nstep(0); + PseudoJet sd_pjet = SoftDrop(pjet, &scatcens, scmap,zg,dR,Nstep,gsj); + + if (sd_pjet != 0) { + _h_grptratio[k]->fill(sd_pjet.pt()/pjet.pt(),weight); + _h_zg[k]->fill(zg,weight); + _h_MassoverpT[k]->fill(sd_pjet.m()/pjet.pt(), weight); + _h_deltar[k]->fill(dR,weight); + _h_nstep[k]->fill(Nstep,weight); + if (dR > 0.1) { + _h_zgdrcut[k]->fill(zg,weight); + _h_MassoverpTdrcut[k]->fill(sd_pjet.m()/pjet.pt(), weight); + if(ptbin == 0){ + _h_zgdrcut_ptbin0[k]->fill(zg,weight); + _h_MassoverpTdrcut_ptbin0[k]->fill(sd_pjet.m()/pjet.pt(), weight); + } + if(ptbin == 4){ + _h_zgdrcut_ptbin4[k]->fill(zg,weight); + _h_MassoverpTdrcut_ptbin4[k]->fill(sd_pjet.m()/pjet.pt(), weight); + } + _h_nstepdrcut[k]->fill(Nstep,weight); + if (Nstep == 0) _h_zgdrcutnstep0[k]->fill(zg,weight); + if (gbin >= 0) _h_zg_girth[gbin][k]->fill(zg,weight); + + PseudoJet subjet1, subjet2; + if (sd_pjet.has_parents(subjet1,subjet2)){ + FourMomentum mom1, mom2; + mom1 = SubtractJetMom(subjet1, &scatcens, scmap); + mom2 = SubtractJetMom(subjet2, &scatcens, scmap); + double pt1, pt2, g1, g2; + pt1 = mom1.pT(); + pt2 = mom2.pT(); + g1 = Girth(subjet1, mom1, &scatcens, scmap); + g2 = Girth(subjet2, mom2, &scatcens, scmap); + if (pt1 > pt2) { + _p_girth_lead_zg[k]->fill(zg,g1,weight); + _p_girth_subl_zg[k]->fill(zg,g2,weight); + } + else { + _p_girth_lead_zg[k]->fill(zg,g2,weight); + _p_girth_subl_zg[k]->fill(zg,g1,weight); + } +#if MODE==0 + double frac1, frac2,ptfrac1,ptfrac2; + frac1 = RecEnFrac(subjet1, mom1, scafter, scmap); + frac2 = RecEnFrac(subjet2, mom2, scafter, scmap); + ptfrac1 = RecPtFrac(subjet1, mom1, scafter, scmap); + ptfrac2 = RecPtFrac(subjet2, mom2, scafter, scmap); + if (pt1 > pt2) { + _p_recen_lead_zg[k]->fill(zg,frac1,weight); + _p_recen_sublead_zg[k]->fill(zg,frac2,weight); + _p_recpt_lead_zg[k]->fill(zg,ptfrac1,weight); + _p_recpt_sublead_zg[k]->fill(zg,ptfrac2,weight); + _p_recpt_abs_lead_zg[k]->fill(zg,ptfrac1*pt1,weight); + _p_recpt_abs_sublead_zg[k]->fill(zg,ptfrac2*pt2,weight); + } + else { + _p_recen_lead_zg[k]->fill(zg,frac2,weight); + _p_recen_sublead_zg[k]->fill(zg,frac1,weight); + _p_recpt_lead_zg[k]->fill(zg,ptfrac2,weight); + _p_recpt_sublead_zg[k]->fill(zg,ptfrac1,weight); + _p_recpt_abs_lead_zg[k]->fill(zg,ptfrac2*pt2,weight); + _p_recpt_abs_sublead_zg[k]->fill(zg,ptfrac1*pt1,weight); + } +#endif + } + } + } + + } + } + + + /*ptsubs.sort(); + ptsubs.reverse(); + ClusterSequence cs2(evnorec, jet_def); + PseudoJets pjets2 = sorted_by_pt(cs2.inclusive_jets(50.0)); + + list::iterator listit(ptsubs.begin()); + for (size_t i=0; ifill(deltapt,weight); + listit++; + }*/ + + } + + } + + + /// Normalise histograms etc., after the run + void finalize() { + + normalize(_h_pt_sc_bef); + normalize(_h_pt_sc_aft); + for (size_t k = 0; k < _njetRbins; ++k) { + scale(_h_jetpt[k], crossSection()/sumOfWeights()); + normalize(_h_zg[k]); + normalize(_h_MassoverpT[k]); + normalize(_h_deltar[k]); + normalize(_h_nstep[k]); + normalize(_h_zgdrcut[k]); + normalize(_h_MassoverpTdrcut[k]); + normalize(_h_zgdrcut_ptbin0[k]); + normalize(_h_MassoverpTdrcut_ptbin0[k]); + normalize(_h_zgdrcut_ptbin4[k]); + normalize(_h_MassoverpTdrcut_ptbin4[k]); + normalize(_h_nstepdrcut[k]); + normalize(_h_zgdrcutnstep0[k]); + normalize(_h_grptratio[k]); + normalize(_h_girth[k]); + for (size_t j = 0; j < _ngirthbins; ++j) { + normalize(_h_zg_girth[j][k]); + } + //normalize(_h_deltapt[k]); + /*for (size_t j = 0; j < _nptbins; ++j) { + normalize(_h_zg[j]); + normalize(_h_nstep[j]); + }*/ + } + + } + + + + private: + + vector _ptedges, _jetRs, _girthedges; + size_t _nptbins, _njetRbins, _ngirthbins; + double _zcut, _beta; + + size_t _nphibins, _netabins; + + Histo1DPtr _h_pt_sc_bef; + Histo1DPtr _h_pt_sc_aft; + Histo1DPtr _h_jetpt[1]; + //Histo1DPtr _h_zg[7]; + //Histo1DPtr _h_nstep[7]; + Histo1DPtr _h_zg[1]; + Histo1DPtr _h_MassoverpT[1]; + Histo1DPtr _h_deltar[1]; + Histo1DPtr _h_nstep[1]; + Histo1DPtr _h_zgdrcut[1]; + Histo1DPtr _h_MassoverpTdrcut[1]; + Histo1DPtr _h_zgdrcut_ptbin0[1]; + Histo1DPtr _h_MassoverpTdrcut_ptbin0[1]; + Histo1DPtr _h_zgdrcut_ptbin4[1]; + Histo1DPtr _h_MassoverpTdrcut_ptbin4[1]; + Histo1DPtr _h_nstepdrcut[1]; + Histo1DPtr _h_zgdrcutnstep0[1]; + Histo1DPtr _h_grptratio[1]; + Histo1DPtr _h_girth[1]; + Histo1DPtr _h_zg_girth[4][1]; + Profile1DPtr _p_girth_lead_zg[1]; + Profile1DPtr _p_girth_subl_zg[1]; +#if MODE==0 + Profile1DPtr _p_recen_lead_zg[1]; + Profile1DPtr _p_recen_sublead_zg[1]; + Profile1DPtr _p_recpt_lead_zg[1]; + Profile1DPtr _p_recpt_sublead_zg[1]; + Profile1DPtr _p_recpt_abs_lead_zg[1]; + Profile1DPtr _p_recpt_abs_sublead_zg[1]; +#endif + //Histo1DPtr _h_deltapt[2]; + + + }; + + + + // The hook for the plugin system +#if MODE==0 + DECLARE_RIVET_PLUGIN(CMS_SOFTDROP_4MOMSUB); +#else + DECLARE_RIVET_PLUGIN(CMS_SOFTDROP_NOREC); +#endif +} Index: trunk/code/CMS_HIN_12_004.cc =================================================================== --- trunk/code/CMS_HIN_12_004.cc (revision 466) +++ trunk/code/CMS_HIN_12_004.cc (revision 467) @@ -1,109 +1,109 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/FastJets.hh" /// @todo Include more projections as required, e.g. ChargedFinalState, FastJets, ZFinder... namespace Rivet { class CMS_HIN_12_004 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor CMS_HIN_12_004() : Analysis("CMS_HIN_12_004") { - //setNeedsCrossSection(true); + setNeedsCrossSection(true); } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { _centedges += 0.0, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9; _ncentbins = _centedges.size()-1; FinalState fs(-5.0, 5.0, 0.*GeV); addProjection(fs, "FS"); FastJets fj(fs, FastJets::ANTIKT, 0.3); fj.useInvisibles(); addProjection(fj, "jets03"); for (size_t i = 0; i < _ncentbins; ++i) { _sumwtcentbins[i]=0.; stringstream ss; ss << "centbin_" << i; const string pname = ss.str(); _histos[i] = bookHisto1D(pname,refData(1, 1, i+1)); _ratios[i] = bookHisto1D(1, 1, i+1); } } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); const double cent = (event.genEvent()->heavy_ion()?event.genEvent()->heavy_ion()->impact_parameter():-1.); Cut cuts = Cuts::abseta < 2.0 && Cuts::pT > 100.*GeV && Cuts::pT < 300.*GeV; const Jets jets = applyProjection(event, "jets03").jetsByPt(cuts); for (size_t i = 0; i < _ncentbins; ++i) { if (cent < 0. || (cent >= _centedges[i] && cent < _centedges[i+1])) { _sumwtcentbins[i]+=weight; // std::cout<fill(jet.momentum().pT(),weight); _ratios[i]->fill(jet.momentum().pT(),weight); } } } } /// Normalise histograms etc., after the run void finalize() { for (size_t i = 0; i < _ncentbins; ++i) { scale(_histos[i],crossSection()/(_sumwtcentbins[i]>0.?_sumwtcentbins[i]:1.)); scale(_ratios[i],crossSection()/(_sumwtcentbins[i]>0.?_sumwtcentbins[i]:1.)); } } //@} private: vector _centedges; double _sumwtcentbins[6]; size_t _ncentbins; private: Histo1DPtr _histos[6]; Histo1DPtr _ratios[6]; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_HIN_12_004); } Index: trunk/code/CMS_GROOMEDMASS_4MOMSUB.cc =================================================================== --- trunk/code/CMS_GROOMEDMASS_4MOMSUB.cc (revision 0) +++ trunk/code/CMS_GROOMEDMASS_4MOMSUB.cc (revision 467) @@ -0,0 +1,273 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/FastJets.hh" +//#include "Rivet/Tools/Logging.hh" +#include "Rivet/Tools/ParticleIdUtils.hh" +//#include "fastjet/AreaDefinition.hh" +#include "fastjet/ClusterSequence.hh" +//#include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/contrib/SoftDrop.hh" +#include "fastjet/contrib/Recluster.hh" +#include "HepMC/GenParticle.h" +#include "HepMC/GenEvent.h" + +namespace Rivet { + + using namespace fastjet; + + struct PJetComp { + bool operator() (const PseudoJet& pj1, const PseudoJet& pj2) const + {return pj1.pt() * scpoint, map & scmap, + double & zcut, double & beta, double & zg, double & dR, int & Nstep, double & Mg) { + int Nsteptot(0); + PseudoJet pjet, parent1, parent2; + pjet = contrib::Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(pj); + while (pjet.has_parents(parent1,parent2)) { + //cout<<"unsubtracted pts: "< 1.) pt1 = 0.; + else pt1 = mom1.pT(); + if (mom2.E() < 0. || deltaR(mom2.eta(), mom2.phi(), Get4Mom(parent2).eta(), Get4Mom(parent2).phi()) > 1.) pt2 = 0.; + else pt2 = mom2.pT(); + //cout<<"subtracted pts: "< zcut*pow(dR/_jetR,beta)) { + double E1 = mom1.E(); + double E2 = mom2.E(); + double theta = angle(mom1.p3(), mom2.p3()); + Mg = sqrt(2.*E1*E2*(1.-cos(theta))); + return pjet; + } + if (pt1 > pt2) { + pjet = parent1; + } + else { + pjet = parent2; + } + Nsteptot += 1; + if (zg > 0.) Nstep += 1; + } + //cout<<"mySoftDrop failed.\n"; + return PseudoJet(); + } + + /// 4-momentum subtraction + FourMomentum SubtractJetMom(PseudoJet pjet, const vector * scatcenters, map & scmap) { + if (scatcenters->empty()) return Get4Mom(pjet); + FourMomentum sub(0.,0.,0.,0.); + bool fillmap(scmap.empty()); + //cout<<"--------------------------------------------------\n"; + //cout<<"size of map: "< 0.) sub += submom; + else sub -= submom; + map::iterator mapit; + mapit = scmap.find(pj); + /*if (mapit != scmap.end()) { + cout<<"Error: dummy is already in map.\n"; + cout<first)<second)<::iterator mapit; + mapit = scmap.find(pj); + if (mapit == scmap.end()) cout<<"Error: did not find matchings scattering centre in map.\n"; + else { + FourMomentum submom = Get4Mom(mapit->second); + if (submom.mass2() > 0.) sub += submom; + else sub -= submom; + } + } + } + } + //cout<<"size of map: "< ptsubs; +//cout<<"==============================\n"; + PseudoJets scatcens; + + vector scafter; + foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { + FourMomentum mom(p->momentum()); + if (fabs(mom.eta()) < 10.) { + if (p->status() == 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + scatcens.push_back(pjet); + } + if(p->status() != 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + pevent.push_back(pjet); + } + } + } + + map scmap; + size_t j(0); + for (size_t i=0; i 1.) ptsub = 0.; + else ptsub = jetsub.pT(); + + int ptbin(-1); + for (size_t l = 0; l < _nptbins; ++l) { + if (ptsub >= _ptedges[l] && ptsub < _ptedges[l+1]) ptbin=l; + } + if (ptbin < 0) continue; + + double zg0(0.), dR0(0.), Mg0(0.); + int Nstep0(0); + double zg1(0.), dR1(0.), Mg1(0.); + int Nstep1(0); + PseudoJet sd_pjet_0 = SoftDrop(pjet, &scatcens, scmap, _zcuts[0], _betas[0], zg0, dR0, Nstep0, Mg0); + PseudoJet sd_pjet_1 = SoftDrop(pjet, &scatcens, scmap, _zcuts[1], _betas[1], zg1, dR1, Nstep1, Mg1); + + if (sd_pjet_0 != 0) { + if (dR0 > 0.1) { + _h_Mgopt_1[ptbin]->fill(Mg0/ptsub,weight); + } + } + if (sd_pjet_1 != 0) { + if (dR1 > 0.1) { + _h_Mgopt_2[ptbin]->fill(Mg1/ptsub,weight); + } + } + } + } + + } + + + /// Normalise histograms etc., after the run + void finalize() { + + for (size_t k = 0; k < _nptbins; ++k) { + normalize(_h_Mgopt_1[k]); + normalize(_h_Mgopt_2[k]); + } + + } + + + + private: + + vector _ptedges, _zcuts, _betas; + size_t _nptbins; + double _jetR; + + Histo1DPtr _h_Mgopt_1[4]; + Histo1DPtr _h_Mgopt_2[4]; + + + }; + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(CMS_GROOMEDMASS_4MOMSUB); +} Index: trunk/code/CMS_2013_I1256590_4MOMSUB.yoda =================================================================== --- trunk/code/CMS_2013_I1256590_4MOMSUB.yoda (revision 466) +++ trunk/code/CMS_2013_I1256590_4MOMSUB.yoda (revision 467) @@ -1,206 +1,206 @@ # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d01-x01-y03 Path=/REF/CMS_2013_I1256590_4MOMSUB/d01-x01-y03 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ -2.500000e-02 2.500000e-02 2.500000e-02 9.996849e-01 1.534981e-02 1.496040e-02 -7.500000e-02 2.500000e-02 2.500000e-02 9.914790e-01 4.180134e-02 4.529186e-02 -1.250000e-01 2.500000e-02 2.500000e-02 9.037203e-01 4.181550e-02 3.772316e-02 -1.750000e-01 2.500000e-02 2.500000e-02 1.016727e+00 5.312961e-02 5.671926e-02 -2.250000e-01 2.500000e-02 2.500000e-02 1.273673e+00 1.439611e-01 1.400741e-01 -2.750000e-01 2.500000e-02 2.500000e-02 1.382885e+00 1.060964e-01 1.097498e-01 +2.500000e-02 2.500000e-02 2.500000e-02 1.02200e+00 4.60000e-02 4.60000e-02 +7.500000e-02 2.500000e-02 2.500000e-02 9.59000e-01 5.30000e-02 5.30000e-02 +1.250000e-01 2.500000e-02 2.500000e-02 8.62000e-01 6.40000e-02 6.50000e-02 +1.750000e-01 2.500000e-02 2.500000e-02 9.60000e-01 9.20000e-02 9.10000e-02 +2.250000e-01 2.500000e-02 2.500000e-02 1.20900e+00 1.37000e-01 1.38000e-01 +2.750000e-01 2.500000e-02 2.500000e-02 1.35700e+00 2.01000e-01 2.02000e-01 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d01-x01-y02 Path=/REF/CMS_2013_I1256590_4MOMSUB/d01-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ -2.500000e-02 2.500000e-02 2.500000e-02 1.266784e+01 0.000000e+00 0.000000e+00 -7.500000e-02 2.500000e-02 2.500000e-02 4.565196e+00 0.000000e+00 0.000000e+00 -1.250000e-01 2.500000e-02 2.500000e-02 1.470847e+00 0.000000e+00 0.000000e+00 -1.750000e-01 2.500000e-02 2.500000e-02 6.714706e-01 0.000000e+00 0.000000e+00 -2.250000e-01 2.500000e-02 2.500000e-02 3.694601e-01 0.000000e+00 0.000000e+00 -2.750000e-01 2.500000e-02 2.500000e-02 2.007718e-01 0.000000e+00 0.000000e+00 +2.500000e-02 2.500000e-02 2.500000e-02 1.23630e+01 2.33000e-02 2.33000e-02 +7.500000e-02 2.500000e-02 2.500000e-02 4.75400e+00 1.34000e-2 1.34000e-02 +1.250000e-01 2.500000e-02 2.500000e-02 1.52800e+00 6.00000e-03 6.00000e-03 +1.750000e-01 2.500000e-02 2.500000e-02 7.36000e-01 4.00000e-03 4.00000e-03 +2.250000e-01 2.500000e-02 2.500000e-02 4.02000e-01 2.00000e-03 2.00000e-03 +2.750000e-01 2.500000e-02 2.500000e-02 2.17000e-01 2.00000e-03 2.00000e-03 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d01-x01-y01 Path=/REF/CMS_2013_I1256590_4MOMSUB/d01-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ -2.500000e-02 2.500000e-02 2.500000e-02 1.266784e+01 0.000000e+00 0.000000e+00 -7.500000e-02 2.500000e-02 2.500000e-02 4.565196e+00 0.000000e+00 0.000000e+00 -1.250000e-01 2.500000e-02 2.500000e-02 1.298714e+00 0.000000e+00 0.000000e+00 -1.750000e-01 2.500000e-02 2.500000e-02 6.883952e-01 0.000000e+00 0.000000e+00 -2.250000e-01 2.500000e-02 2.500000e-02 4.738880e-01 0.000000e+00 0.000000e+00 -2.750000e-01 2.500000e-02 2.500000e-02 2.774874e-01 0.000000e+00 0.000000e+00 +2.500000e-02 2.500000e-02 2.500000e-02 1.26360e+01 9.84000e-01 9.83000e-01 +7.500000e-02 2.500000e-02 2.500000e-02 4.56000e+00 3.99000e-01 3.97000e-01 +1.250000e-01 2.500000e-02 2.500000e-02 1.31800e+00 1.41000e-01 1.40000e-01 +1.750000e-01 2.500000e-02 2.500000e-02 7.07000e-01 8.70000e-02 8.60000e-02 +2.250000e-01 2.500000e-02 2.500000e-02 4.86000e-01 6.50000e-02 6.60000e-02 +2.750000e-01 2.500000e-02 2.500000e-02 2.94000e-01 4.90000e-02 5.10000e-02 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d03-x01-y01 Path=/REF/CMS_2013_I1256590_4MOMSUB/d03-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-01 2.500000e-01 2.500000e-01 4.589493e-02 2.640026e-02 2.607364e-02 7.500000e-01 2.500000e-01 2.500000e-01 3.302992e-01 1.372295e-01 1.308755e-01 1.250000e+00 2.500000e-01 2.500000e-01 7.553553e-01 2.064520e-01 2.113501e-01 1.750000e+00 2.500000e-01 2.500000e-01 1.219367e+00 1.948854e-01 1.904400e-01 2.250000e+00 2.500000e-01 2.500000e-01 1.560549e+00 1.710532e-01 1.668595e-01 2.750000e+00 2.500000e-01 2.500000e-01 1.804274e+00 2.437255e-01 2.817904e-01 3.250000e+00 2.500000e-01 2.500000e-01 2.086064e+00 4.080724e-01 3.258002e-01 3.750000e+00 2.500000e-01 2.500000e-01 2.447121e+00 4.787017e-01 4.655116e-01 4.250000e+00 2.500000e-01 2.500000e-01 2.556013e+00 3.771229e-01 3.566198e-01 4.750000e+00 2.500000e-01 2.500000e-01 1.273626e+00 4.731263e-01 4.790329e-01 5.250000e+00 2.500000e-01 2.500000e-01 1.213483e-01 5.143852e-02 1.293548e-01 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d03-x01-y03 Path=/REF/CMS_2013_I1256590_4MOMSUB/d03-x01-y03 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-01 2.500000e-01 2.500000e-01 9.511278e-01 2.506266e-01 2.631579e-01 7.500000e-01 2.500000e-01 2.500000e-01 9.949875e-01 1.817043e-01 1.879699e-01 1.250000e+00 2.500000e-01 2.500000e-01 8.508772e-01 1.190476e-01 1.190476e-01 1.750000e+00 2.500000e-01 2.500000e-01 8.132832e-01 1.002506e-01 9.398496e-02 2.250000e+00 2.500000e-01 2.500000e-01 7.380952e-01 8.145363e-02 8.145363e-02 2.750000e+00 2.500000e-01 2.500000e-01 7.506266e-01 8.145363e-02 8.145363e-02 3.250000e+00 2.500000e-01 2.500000e-01 8.508772e-01 8.771930e-02 8.771930e-02 3.750000e+00 2.500000e-01 2.500000e-01 1.095238e+00 1.190476e-01 1.127820e-01 4.250000e+00 2.500000e-01 2.500000e-01 1.508772e+00 1.879699e-01 1.879699e-01 4.750000e+00 2.500000e-01 2.500000e-01 1.765664e+00 3.320802e-01 3.320802e-01 5.250000e+00 2.500000e-01 2.500000e-01 2.223058e+00 6.892231e-01 6.892231e-01 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d03-x01-y02 Path=/REF/CMS_2013_I1256590_4MOMSUB/d03-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-01 2.500000e-01 2.500000e-01 4.793716e-02 0.000000e+00 0.000000e+00 7.500000e-01 2.500000e-01 2.500000e-01 3.302992e-01 0.000000e+00 0.000000e+00 1.250000e+00 2.500000e-01 2.500000e-01 8.990452e-01 0.000000e+00 0.000000e+00 1.750000e+00 2.500000e-01 2.500000e-01 1.494066e+00 0.000000e+00 0.000000e+00 2.250000e+00 2.500000e-01 2.500000e-01 2.147498e+00 0.000000e+00 0.000000e+00 2.750000e+00 2.500000e-01 2.500000e-01 2.447121e+00 0.000000e+00 0.000000e+00 3.250000e+00 2.500000e-01 2.500000e-01 2.447121e+00 0.000000e+00 0.000000e+00 3.750000e+00 2.500000e-01 2.500000e-01 2.243057e+00 0.000000e+00 0.000000e+00 4.250000e+00 2.500000e-01 2.500000e-01 1.702521e+00 0.000000e+00 0.000000e+00 4.750000e+00 2.500000e-01 2.500000e-01 7.231755e-01 0.000000e+00 0.000000e+00 5.250000e+00 2.500000e-01 2.500000e-01 5.462544e-02 0.000000e+00 0.000000e+00 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d04-x01-y02 Path=/REF/CMS_2013_I1256590_4MOMSUB/d04-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 1.097304e+00 9.730405e-02 9.730405e-02 9.843526e-01 0.000000e+00 0.000000e+00 1.310850e+00 1.162410e-01 1.162410e-01 1.032045e+00 0.000000e+00 0.000000e+00 1.565954e+00 1.388630e-01 1.388630e-01 1.032045e+00 0.000000e+00 0.000000e+00 1.870703e+00 1.658865e-01 1.658865e-01 9.689500e-01 0.000000e+00 0.000000e+00 2.234760e+00 1.981700e-01 1.981700e-01 9.537884e-01 0.000000e+00 0.000000e+00 2.669665e+00 2.367350e-01 2.367350e-01 8.407319e-01 0.000000e+00 0.000000e+00 3.189208e+00 2.828065e-01 2.828065e-01 7.648242e-01 0.000000e+00 0.000000e+00 3.809856e+00 3.378430e-01 3.378430e-01 6.132976e-01 0.000000e+00 0.000000e+00 4.551290e+00 4.035905e-01 4.035905e-01 5.238152e-01 0.000000e+00 0.000000e+00 5.437013e+00 4.821335e-01 4.821335e-01 4.473886e-01 0.000000e+00 0.000000e+00 6.495107e+00 5.759600e-01 5.759600e-01 3.644549e-01 0.000000e+00 0.000000e+00 7.759116e+00 6.880480e-01 6.880480e-01 2.922493e-01 0.000000e+00 0.000000e+00 9.269112e+00 8.219480e-01 8.219480e-01 2.418587e-01 0.000000e+00 0.000000e+00 1.107298e+01 9.819200e-01 9.819200e-01 1.764312e-01 0.000000e+00 0.000000e+00 1.322790e+01 1.173000e+00 1.173000e+00 1.370839e-01 0.000000e+00 0.000000e+00 1.580216e+01 1.401255e+00 1.401255e+00 1.015896e-01 0.000000e+00 0.000000e+00 1.887738e+01 1.673970e+00 1.673970e+00 6.957701e-02 0.000000e+00 0.000000e+00 2.255110e+01 1.999745e+00 1.999745e+00 4.917913e-02 0.000000e+00 0.000000e+00 2.693975e+01 2.388910e+00 2.388910e+00 3.368195e-02 0.000000e+00 0.000000e+00 3.218248e+01 2.853820e+00 2.853820e+00 2.306820e-02 0.000000e+00 0.000000e+00 3.844549e+01 3.409190e+00 3.409190e+00 1.287031e-02 0.000000e+00 0.000000e+00 4.592735e+01 4.072655e+00 4.072655e+00 6.532336e-03 0.000000e+00 0.000000e+00 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d04-x01-y03 Path=/REF/CMS_2013_I1256590_4MOMSUB/d04-x01-y03 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 1.097304e+00 9.730405e-02 9.730405e-02 1.050239e+00 5.598086e-01 5.645933e-01 1.310850e+00 1.162410e-01 1.162410e-01 6.339713e-01 4.354067e-01 4.354067e-01 1.565954e+00 1.388630e-01 1.388630e-01 6.818182e-01 3.636364e-01 3.684211e-01 1.870703e+00 1.658865e-01 1.658865e-01 4.043062e-01 2.918660e-01 2.918660e-01 2.234760e+00 1.981700e-01 1.981700e-01 2.224880e-01 2.344498e-01 2.392344e-01 2.669665e+00 2.367350e-01 2.367350e-01 5.502392e-02 1.818182e-01 1.866029e-01 3.189208e+00 2.828065e-01 2.828065e-01 -5.502392e-02 1.387560e-01 1.435407e-01 3.809856e+00 3.378430e-01 3.378430e-01 -6.937799e-02 1.052632e-01 1.004785e-01 4.551290e+00 4.035905e-01 4.035905e-01 -9.808612e-02 7.177033e-02 7.177033e-02 5.437013e+00 4.821335e-01 4.821335e-01 -1.028708e-01 5.263158e-02 5.263158e-02 6.495107e+00 5.759600e-01 5.759600e-01 -8.851675e-02 3.827751e-02 3.349282e-02 7.759116e+00 6.880480e-01 6.880480e-01 -8.373206e-02 2.870813e-02 2.392344e-02 9.269112e+00 8.219480e-01 8.219480e-01 -6.937799e-02 1.435407e-02 2.392344e-02 1.107298e+01 9.819200e-01 9.819200e-01 -4.066986e-02 1.913876e-02 1.435407e-02 1.322790e+01 1.173000e+00 1.173000e+00 -3.588517e-02 1.435407e-02 1.435407e-02 1.580216e+01 1.401255e+00 1.401255e+00 -2.631579e-02 9.569378e-03 1.435407e-02 1.887738e+01 1.673970e+00 1.673970e+00 -7.177033e-03 1.435407e-02 4.784689e-03 2.255110e+01 1.999745e+00 1.999745e+00 -7.177033e-03 9.569378e-03 4.784689e-03 2.693975e+01 2.388910e+00 2.388910e+00 -2.392344e-03 0.000000e+00 0.000000e+00 3.218248e+01 2.853820e+00 2.853820e+00 2.392344e-03 0.000000e+00 0.000000e+00 3.844549e+01 3.409190e+00 3.409190e+00 2.392344e-03 0.000000e+00 0.000000e+00 4.592735e+01 4.072655e+00 4.072655e+00 -2.392344e-03 0.000000e+00 0.000000e+00 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d04-x01-y01 Path=/REF/CMS_2013_I1256590_4MOMSUB/d04-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 1.097304e+00 9.730405e-02 9.730405e-02 2.033383e+00 4.782024e-01 4.627072e-01 1.310850e+00 1.162410e-01 1.162410e-01 1.656449e+00 3.489591e-01 3.451167e-01 1.565954e+00 1.388630e-01 1.388630e-01 1.709530e+00 2.947629e-01 3.238530e-01 1.870703e+00 1.658865e-01 1.658865e-01 1.370839e+00 2.541162e-01 2.341778e-01 2.234760e+00 1.981700e-01 1.981700e-01 1.189440e+00 2.204898e-01 1.813987e-01 2.669665e+00 2.367350e-01 2.367350e-01 8.954777e-01 1.774117e-01 1.696391e-01 3.189208e+00 2.828065e-01 2.828065e-01 7.068302e-01 1.400369e-01 1.339018e-01 3.809856e+00 3.378430e-01 3.378430e-01 5.491943e-01 1.088062e-01 1.040393e-01 4.551290e+00 4.035905e-01 4.035905e-01 4.267141e-01 8.454043e-02 8.890473e-02 5.437013e+00 4.821335e-01 4.821335e-01 3.421737e-01 6.779131e-02 5.845289e-02 6.495107e+00 5.759600e-01 5.759600e-01 2.700890e-01 4.656965e-02 4.613880e-02 7.759116e+00 6.880480e-01 6.880480e-01 2.065706e-01 3.561760e-02 3.528807e-02 9.269112e+00 8.219480e-01 8.219480e-01 1.709530e-01 2.947629e-02 2.607165e-02 1.107298e+01 9.819200e-01 9.819200e-01 1.328274e-01 2.115517e-02 2.269066e-02 1.322790e+01 1.173000e+00 1.173000e+00 1.015896e-01 1.883195e-02 1.924512e-02 1.580216e+01 1.401255e+00 1.401255e+00 7.648242e-02 1.611232e-02 1.593491e-02 1.887738e+01 1.673970e+00 1.673970e+00 5.667932e-02 1.332960e-02 1.073732e-02 2.255110e+01 1.999745e+00 1.999745e+00 4.006265e-02 1.083773e-02 1.149923e-02 2.693975e+01 2.388910e+00 2.388910e+00 2.831749e-02 8.615028e-03 9.295896e-03 3.218248e+01 2.853820e+00 2.853820e+00 1.909070e-02 6.620015e-03 6.266981e-03 3.844549e+01 3.409190e+00 3.409190e+00 1.170828e-02 4.750580e-03 4.856209e-03 4.592735e+01 4.072655e+00 4.072655e+00 6.636175e-03 2.933691e-03 3.207351e-03 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d02-x01-y01 Path=/REF/CMS_2013_I1256590_4MOMSUB/d02-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-02 2.500000e-02 2.500000e-02 6.375449e-01 0.000000e+00 0.000000e+00 7.500000e-02 2.500000e-02 2.500000e-02 8.643001e-01 0.000000e+00 0.000000e+00 1.250000e-01 2.500000e-02 2.500000e-02 9.278651e-01 0.000000e+00 0.000000e+00 1.750000e-01 2.500000e-02 2.500000e-02 9.612714e-01 0.000000e+00 0.000000e+00 2.250000e-01 2.500000e-02 2.500000e-02 9.858084e-01 0.000000e+00 0.000000e+00 2.750000e-01 2.500000e-02 2.500000e-02 9.979320e-01 0.000000e+00 0.000000e+00 # END YODA_SCATTER2D # BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590_4MOMSUB/d02-x01-y02 Path=/REF/CMS_2013_I1256590_4MOMSUB/d02-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-02 2.500000e-02 2.500000e-02 6.357711e-01 0.000000e+00 0.000000e+00 7.500000e-02 2.500000e-02 2.500000e-02 8.678478e-01 0.000000e+00 0.000000e+00 1.250000e-01 2.500000e-02 2.500000e-02 9.366820e-01 0.000000e+00 0.000000e+00 1.750000e-01 2.500000e-02 2.500000e-02 9.700882e-01 0.000000e+00 0.000000e+00 2.250000e-01 2.500000e-02 2.500000e-02 9.892897e-01 0.000000e+00 0.000000e+00 2.750000e-01 2.500000e-02 2.500000e-02 9.978656e-01 0.000000e+00 0.000000e+00 # END YODA_SCATTER2D Index: trunk/code/CMS_2012_I1090064.plot =================================================================== --- trunk/code/CMS_2012_I1090064.plot (revision 466) +++ trunk/code/CMS_2012_I1090064.plot (revision 467) @@ -1,215 +1,231 @@ # BEGIN PLOT /CMS_2012_I1090064/* FullRange=1 # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d01-x01-y0* +#NormalizeToIntegral=1 Scale=0.10472 # END PLOT +# BEGIN PLOT /CMS_2012_I1090064/d03-x01-y0* +#NormalizeToIntegral=1 +Scale=0.06 +# END PLOT + +# BEGIN PLOT /CMS_2012_I1090064/d04-x01-y0* +#NormalizeToIntegral=1 +Scale=0.06 +# END PLOT + +# BEGIN PLOT /CMS_2012_I1090064/d05-x01-y0* +#NormalizeToIntegral=1 +Scale=0.1 +# END PLOT + # BEGIN PLOT /CMS_2012_I1090064/d01-x01-y01 Title=PbPb 0-20, $120\ \mathrm{GeV} < p_{\perp,1} < 150\ \mathrm{GeV}$ XLabel=$\Delta \phi$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d01-x01-y02 Title=PbPb 0-20, $150\ \mathrm{GeV} < p_{\perp,1} < 180\ \mathrm{GeV}$ XLabel=$\Delta \phi$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d01-x01-y03 Title=PbPb 0-20, $180\ \mathrm{GeV} < p_{\perp,1} < 220\ \mathrm{GeV}$ XLabel=$\Delta \phi$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d01-x01-y04 Title=PbPb 0-20, $220\ \mathrm{GeV} < p_{\perp,1} < 260\ \mathrm{GeV}$ XLabel=$\Delta \phi$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d01-x01-y05 #Title=PbPb 0-20, $260\ \mathrm{GeV} < p_{\perp,1} < 300\ \mathrm{GeV}$ XLabel=$\Delta \phi$ YLabel=event fraction END PLOT # BEGIN PLOT /CMS_2012_I1090064/d01-x01-y06 Title=PbPb 0-20, $300\ \mathrm{GeV} < p_{\perp,1} < 500\ \mathrm{GeV}$ XLabel=$\Delta \phi$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d02-x01-y01 LogY=0 XMax=220. YMin=0.7 Title=di-jet fraction in pp XLabel=$p_{\perp,1}$ [GeV] YLabel=$N_\mathrm{dijet}/N_\mathrm{leading jet}$ # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d02-x01-y02 LogY=0 XMax=370. YMin=0.7 Title=di-jet fraction in PbPb 0-20 XLabel=$p_{\perp,1}$ [GeV] YLabel=$N_\mathrm{dijet}/N_\mathrm{leading jet}$ # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d03-x01-y0* Scale=0.06 LogY=0 # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d03-x01-y01 Title=di-jet asymmetry in PbPb 0-20, $120\ \mathrm{GeV} < p_{\perp,1} < 150\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d03-x01-y02 Title=di-jet asymmetry in PbPb 0-20, $150\ \mathrm{GeV} < p_{\perp,1} < 180\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d03-x01-y03 Title=di-jet asymmetry in PbPb 0-20, $180\ \mathrm{GeV} < p_{\perp,1} < 220\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d03-x01-y04 Title=di-jet asymmetry in PbPb 0-20, $220\ \mathrm{GeV} < p_{\perp,1} < 260\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d03-x01-y05 Title=di-jet asymmetry in PbPb 0-20, $260\ \mathrm{GeV} < p_{\perp,1} < 300\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d03-x01-y06 Title=di-jet asymmetry in PbPb 0-20, $300\ \mathrm{GeV} < p_{\perp,1} < 500\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d04-x01-y0* Scale=0.06 LogY=0 # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d04-x01-y01 Title=di-jet asymmetry in PbPb 0-10, $p_{\perp,1} > 120\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d04-x01-y02 Title=di-jet asymmetry in PbPb 10-20, $p_{\perp,1} > 120\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d04-x01-y03 Title=di-jet asymmetry in PbPb 20-30, $p_{\perp,1} > 120\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d04-x01-y04 Title=di-jet asymmetry in PbPb 30-50, $p_{\perp,1} > 120\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d04-x01-y05 Title=di-jet asymmetry in PbPb 50-70, $p_{\perp,1} > 120\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d04-x01-y06 Title=di-jet asymmetry in PbPb 70-100, $p_{\perp,1} > 120\ \mathrm{GeV}$ XLabel=$A_J$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d05-x01-y0* Scale=0.1 LogY=0 # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d05-x01-y01 Title=di-jet $p_\perp$ ratio in PbPb 0-20, $120\ \mathrm{GeV} < p_{\perp,1} < 150\ \mathrm{GeV}$ XLabel=$p_{\perp,2}/p_{\perp,1}$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d05-x01-y02 Title=di-jet $p_\perp$ ratio in PbPb 0-20, $150\ \mathrm{GeV} < p_{\perp,1} < 180\ \mathrm{GeV}$ XLabel=$p_{\perp,2}/p_{\perp,1}$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d05-x01-y03 Title=di-jet $p_\perp$ ratio in PbPb 0-20, $180\ \mathrm{GeV} < p_{\perp,1} < 220\ \mathrm{GeV}$ XLabel=$p_{\perp,2}/p_{\perp,1}$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d05-x01-y04 Title=di-jet $p_\perp$ ratio in PbPb 0-20, $220\ \mathrm{GeV} < p_{\perp,1} < 260\ \mathrm{GeV}$ XLabel=$p_{\perp,2}/p_{\perp,1}$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d05-x01-y05 Title=di-jet $p_\perp$ ratio in PbPb 0-20, $260\ \mathrm{GeV} < p_{\perp,1} < 300\ \mathrm{GeV}$ XLabel=$p_{\perp,2}/p_{\perp,1}$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d05-x01-y06 Title=di-jet $p_\perp$ ratio in PbPb 0-20, $300\ \mathrm{GeV} < p_{\perp,1} < 500\ \mathrm{GeV}$ XLabel=$p_{\perp,2}/p_{\perp,1}$ YLabel=event fraction # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d06-x01-y0* YMin=0.5 XMax=300. LogY=0 # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d06-x01-y01 Title=mean di-jet $p_\perp$ ratio in pp XLabel=$p_{\perp,1}$ [GeV] YLabel=$\langle p_{\perp,2}/p_{\perp,1}\rangle$ # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d06-x01-y02 Title=mean di-jet $p_\perp$ ratio in PbPb 0-20 XLabel=$p_{\perp,1}$ [GeV] YLabel=$\langle p_{\perp,2}/p_{\perp,1}\rangle$ # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d06-x01-y03 Title=mean di-jet $p_\perp$ ratio in PbPb 20-50 XLabel=$p_{\perp,1}$ [GeV] YLabel=$\langle p_{\perp,2}/p_{\perp,1}\rangle$ # END PLOT # BEGIN PLOT /CMS_2012_I1090064/d06-x01-y04 Title=mean di-jet $p_\perp$ ratio in PbPb 50-100 XLabel=$p_{\perp,1}$ [GeV] YLabel=$\langle p_{\perp,2}/p_{\perp,1}\rangle$ # END PLOT Index: trunk/code/CMS_2013_I1256590_CONSTSUB.plot =================================================================== --- trunk/code/CMS_2013_I1256590_CONSTSUB.plot (revision 0) +++ trunk/code/CMS_2013_I1256590_CONSTSUB.plot (revision 467) @@ -0,0 +1,81 @@ +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d01-x01-y01 +Title=PbPb 0-10 +XLabel=$r$ +YLabel=$\rho(r)$ +LogY=1 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d01-x01-y02 +Title=pp smeared to PbPb resolution +XLabel=$r$ +YLabel=$\rho(r)$ +LogY=1 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d01-x01-y03 +Title=PbPb(0-10)/pp +XLabel=$r$ +YLabel=$\rho_\text{PbPb}/\rho_\text{pp}$ +LogY=0 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d02-x01-y01 +Title=PbPb 0-10 +XLabel=$r$ +YLabel=$\Psi(r)$ +LogY=0 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d02-x01-y02 +Title=pp smeared to PbPb resolution +XLabel=$r$ +YLabel=$\Psi(r)$ +LogY=0 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d03-x01-y01 +Title=PbPb 0-10 +XLabel=$\xi = \ln(1/z)$ +YLabel=$1/N_\text{jet}\, \text{d}N_\text{track}/\text{d}\xi$ +LogY=1 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d03-x01-y02 +Title=pp smeared to PbPb resolution +XLabel=$\xi = \ln(1/z)$ +YLabel=$1/N_\text{jet}\, \text{d}N_\text{track}/\text{d}\xi$ +LogY=1 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d03-x01-y03 +Title=PbPb(0-10)/pp +XLabel=$\xi = \ln(1/z)$ +YLabel=$D_\text{PbPb}(\xi)/D_\text{pp}(\xi)$ +LogY=0 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d04-x01-y01 +Title=PbPb 0-10 +XLabel=$p_\perp\, [\text{GeV}]$ +YLabel=$1/N_\text{jet}\, \text{d}N_\text{track}/\text{d}p_\perp\, [\text{GeV}^{-1}]$ +LogX=1 +LogY=1 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d04-x01-y02 +Title=pp smeared to PbPb resolution +XLabel=$p_\perp\, [\text{GeV}]$ +YLabel=$1/N_\text{jet}\, \text{d}N_\text{track}/\text{d}p_\perp\, [\text{GeV}^{-1}]$ +LogX=1 +LogY=1 +# END PLOT + +# BEGIN PLOT /CMS_2013_I1256590_CONSTSUB/d04-x01-y03 +Title=PbPb(0-10)-pp +XLabel=$p_\perp\, [\text{GeV}]$ +YLabel=$D_\text{PbPb}(p_\perp)-D_\text{pp}(p_\perp)$ +LogX=1 +LogY=0 +# END PLOT + + Index: trunk/code/CMS_2013_I1256590.cc =================================================================== --- trunk/code/CMS_2013_I1256590.cc (revision 466) +++ trunk/code/CMS_2013_I1256590.cc (revision 467) @@ -1,444 +1,607 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include #include #include /// @todo Include more projections as required, e.g. ChargedFinalState, FastJets, ZFinder... namespace Rivet { using namespace fastjet; struct FourMomComp { bool operator() (const FourMomentum& p1, const FourMomentum& p2) const {return p1.pT() 0.5) n=0.001; if (cent > 0.3 && cent < 0.5) n=3.88; if (cent > 0.1 && cent < 0.3) n=5.10; if (cent > 0. && cent < 0.1) n=5.23; n=5.23; double sigma(sqrt(c*c+s*s/pttrue+n*n/(pttrue*pttrue))); double r1(1.0*rand()/RAND_MAX), r2(1.0*rand()/RAND_MAX); double fac(max(sqrt(-2.*log(r1))*cos(2.*M_PI*r2)*sigma+1.,0.)); return fac*pttrue; } /// extract thermal momenta, that should be subtracted, from HepMC event vector extractThermalMomenta(const Event & event) { vector thermom; foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { FourMomentum mom(p->momentum()); if (p->status() == 3) { thermom.push_back(mom); } } return thermom; } /// build map from dummy particles to thermal momenta to speed up subtraction map buildThMomMap(const Event & event, const vector * thermom) { map thmap; foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { FourMomentum mom(p->momentum()); if (mom.E() < 1e-5 && mom.E() > 1e-7) { FourMomentum dummy(mom); dummy *= 10000.; double dRmin(1.); FourMomentum match; foreach (FourMomentum tm, *thermom) { double dR(deltaR(dummy, tm)); if (dR < dRmin) { dRmin = dR; match = tm; } } if (dRmin < 1e-5) { map::iterator mapit; mapit = thmap.find(mom); if (mapit != thmap.end()) { cout<<"Error: dummy is already in map.\n"; cout<first<second< buildThMomMap2(const Event & event, const vector * thermom) { + map thmap; + size_t j(0); + foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { + FourMomentum mom(p->momentum()); + if (mom.E() == 1e-6) { + if (j == thermom->size()) { + cout<<"Error: number of dummies does not match number of subtractions, will veto Event"< 1e-4) continue; + thmap[mom]=sc; + j++; + } + } + return thmap; + } + /// 4-momentum subtraction FourMomentum SubtractJetMom(Jet jet, map * tmmap, double fac = 1.) { if (tmmap->empty()) return jet.momentum(); FourMomentum sub(0.,0.,0.,0.); foreach (Particle part, jet.constituents()) { if (part.E() < 1e-5 && part.E() > 1e-7) { map::iterator mapit; mapit = (*tmmap).find(part.momentum()); if (mapit == (*tmmap).end()) cout<<"Error: did not find matching scattering centre in map.\n"<second; - } + else { + FourMomentum submom(mapit->second); + if (submom.mass2() < 0.) sub -= submom; + /*{ + if ((submom.E() < 0. && submom.eta() == -part.momentum().eta()) || (submom.E() >= 0. && submom.eta() == part.momentum().eta())) { + sub -= submom; + } + else { + sub += submom; + } + }*/ + else { + sub += submom; + } + } + } } FourMomentum jetmom(jet.momentum()); //cout<<"Original momentum: "< ConstSubJetMom(Jet jet, map & scmap) { + // sort constituents into two vectors: particles and ghosts (subtraction momenta) + vector parts; + vector ghosts; + foreach (PseudoJet pj, pjet.constituents()) { + if (fuzzyEquals(pj.E(),1e-6,1e-4)) { + map::iterator mapit; + mapit = scmap.find(pj); + if (mapit == scmap.end()) cout<<"Error: did not find matching scattering centre in map.\n"; + else { + FourMomentum submom = Get4Mom(mapit->second); + My4Mom ghost; + ghost.pt = submom.pT(); + ghost.mdelta = sqrt(submom.mass2() + sqr(submom.pT())) - submom.pT(); + ghost.phi = submom.p3().phi(); + ghost.y = submom.rapidity(); + ghosts.push_back(ghost); + } + } + else { + FourMomentum partmom = Get4Mom(pj); + My4Mom part; + part.pt = partmom.pT(); + part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT(); + part.phi = partmom.p3().phi(); + part.y = partmom.rapidity(); + parts.push_back(part); + } + } + //cout<<"number of particles: "< dists; + for (size_t i=0; i < parts.size(); ++i) { + for (size_t j=0; j < ghosts.size(); ++j) { + //cout<<"i, j = "<::iterator liter=dists.begin(); liter != dists.end(); ++liter) { + //cout<<"dealing with dist "<dR<ipart; + size_t gnum = liter->ighost; + double ptp = parts[pnum].pt; + double ptg = ghosts[gnum].pt; + //cout<<"pts: "< ptg) { + parts[pnum].pt -= ptg; + ghosts[gnum].pt = 0.; + } + else { + ghosts[gnum].pt -= ptp; + parts[pnum].pt = 0.; + } + double mdp = parts[pnum].mdelta; + double mdg = ghosts[gnum].mdelta; + //cout<<"masses: "< mdg) { + parts[pnum].mdelta -= mdg; + ghosts[gnum].mdelta = 0.; + } + else { + ghosts[gnum].mdelta -= mdp; + parts[pnum].mdelta = 0.; + } + } + + //sum up resulting 4-momenta to get subtracted jet momentum + FourMomentum partmom, ghostmom, jetmom; + for (size_t i=0; i < parts.size(); ++i) { + if (parts[i].pt > 0.) partmom += FourMomentum((parts[i].pt+parts[i].mdelta)*cosh(parts[i].y), + parts[i].pt*cos(parts[i].phi), + parts[i].pt*sin(parts[i].phi), + (parts[i].pt+parts[i].mdelta)*sinh(parts[i].y)); + } + for (size_t i=0; i < ghosts.size(); ++i) { + if (ghosts[i].pt > 0.) ghostmom += FourMomentum((ghosts[i].pt+ghosts[i].mdelta)*cosh(ghosts[i].y), + ghosts[i].pt*cos(ghosts[i].phi), + ghosts[i].pt*sin(ghosts[i].phi), + (ghosts[i].pt+ghosts[i].mdelta)*sinh(ghosts[i].y)); + } + //cout<<"particle component: "<particles(); else parts = jet.constituents(); double scalefac; double ntracks(0.), nparts(0.); FourMomentum chmom(0.,0.,0.,0.), allmom(0.,0.,0.,0.); foreach (Particle part, parts) { if (deltaR(part,jet) < _jetR && part.momentum().E() > 1e-5) { nparts += 1.; allmom += part.momentum(); if (PID::charge(part.pdgId()) != 0 && part.pT() > ptcut) { ntracks += 1.; chmom += part.momentum(); } } } scalefac = (nparts>0.?ntracks/nparts:1); //scalefac = (allmom.pT()>0.?chmom.pT()/allmom.pT():1.); return scalefac; } /// initialise and fill grid for grid subtration vector > fillGrid(const Event & event, std::set > & gridevent) { // initialise grid vector > grid; grid.resize(_nphibins); for (size_t i = 0; i < _nphibins; ++i) { grid[i].resize(_netabins); for (size_t k = 0; k < _nphibins; ++k) { grid[i][k] = FourMomentum(0., 0., 0., 0.); } } // fill grid foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { FourMomentum mom(p->momentum()); if (fabs(mom.eta()) < 4.) { size_t phibin, etabin; phibin = int(mom.phi(ZERO_2PI)*_nphibins/(2.*M_PI)); etabin = int((mom.eta() + 4.)*_netabins/8.); if (phibin < 0 || phibin > _nphibins-1) {std::cout<<"Error: "< "< _netabins-1) {std::cout<<"Error: "< "<status() == 3) { grid[phibin][etabin].setE(grid[phibin][etabin].E() - mom.E()); grid[phibin][etabin].setPx(grid[phibin][etabin].px() - mom.px()); grid[phibin][etabin].setPy(grid[phibin][etabin].py() - mom.py()); grid[phibin][etabin].setPz(grid[phibin][etabin].pz() - mom.pz()); } else { grid[phibin][etabin].setE(grid[phibin][etabin].E() + mom.E()); grid[phibin][etabin].setPx(grid[phibin][etabin].px() + mom.px()); grid[phibin][etabin].setPy(grid[phibin][etabin].py() + mom.py()); grid[phibin][etabin].setPz(grid[phibin][etabin].pz() + mom.pz()); } /*double cellphi((phibin+0.5)*2*M_PI/_nphibins), celleta((etabin+0.5)*8./_netabins-4.); double celltheta(2.*atan(exp(-celleta))); if (p->status() == 3) { double rho(grid[phibin][etabin].E() - mom.E()); grid[phibin][etabin].setE(rho); grid[phibin][etabin].setPx(fabs(rho)*sin(celltheta)*cos(cellphi)); grid[phibin][etabin].setPy(fabs(rho)*sin(celltheta)*sin(cellphi)); grid[phibin][etabin].setPz(fabs(rho)*cos(celltheta)); } else { // use this for standard recoil treatment double rho(grid[phibin][etabin].E() + mom.E()); grid[phibin][etabin].setE(rho); grid[phibin][etabin].setPx(rho*sin(celltheta)*cos(cellphi)); grid[phibin][etabin].setPy(rho*sin(celltheta)*sin(cellphi)); grid[phibin][etabin].setPz(rho*cos(celltheta)); }*/ gridevent.insert(std::pair(phibin,etabin)); } } return grid; } /// 4-momentum adapter FourMomentum Get4Mom(PseudoJet pjet){ FourMomentum mom(pjet.E(), pjet.px(), pjet.py(), pjet.pz()); return mom; } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); const double cent = (event.genEvent()->heavy_ion()?event.genEvent()->heavy_ion()->impact_parameter():-1.); const Jets jets03 = applyProjection(event, "Jets03").jetsByPt(3.*GeV); const FinalState fs = applyProjection(event, "FS"); Jets jetssub; #if MODE!=1 vector thermom = extractThermalMomenta(event); - map thmommap = buildThMomMap(event, &thermom); + map thmommap = buildThMomMap2(event, &thermom); + if (thmommap.size() != thermom.size()) { + cout<<"Error: size of map does not match number of momenta: "< > gridevent; vector > grid = fillGrid(event, gridevent); for (set >::iterator it=gridevent.begin(); it!=gridevent.end(); ++it) { size_t i,k; i = it->first; k = it->second; if (grid[i][k].E() > 0. ) { PseudoJet part(grid[i][k].px(),grid[i][k].py(),grid[i][k].pz(),grid[i][k].E()); pevent.push_back(part); } } JetDefinition jet_def(antikt_algorithm, _jetR); ClusterSequence cs(pevent, jet_def); pjets = sorted_by_pt(cs.inclusive_jets(3.*GeV)); foreach (PseudoJet pjet, pjets) { jetssub.push_back(Jet(pjet)); } #endif map jetmap; foreach (Jet jet, jetssub) { double ptreco, ptgen; ptgen = jet.momentum().pT(); - ptreco = ptsmear(ptgen,cent); + //ptreco = ptsmear(ptgen,cent); + ptreco = ptsmear(ptgen,0.05); jetmap[jet] = ptreco; } Jets jets; foreach (Jet jet, jetssub) { if (fabs(jet.momentum().eta()) > 0.3 && fabs(jet.momentum().eta()) < 2.0 && jetmap[jet] > 100.*GeV) jets.push_back(jet); } // if (jets.size() == 0) vetoEvent; foreach (Jet jet, jets) { vector difvec(_nbins, 0.), difvecvac(_nbins, 0.); vector intvec(_nbins, 0.), intvecvac(_nbins, 0.); _Njets+=weight; double jetpt = jetmap[jet]; -#if MODE==0 +#if (MODE==0 || MODE==4) double scalefac=1.; #elif MODE==2 double scalefac = getScaleFac(jet, _trackptmin); _h_scalefac->fill(scalefac, weight); #endif int bin; foreach (Particle part, fs.particles()){ double dR(deltaR(part,jet)); if (dR >= _jetR) continue; bin = int(dR*_nbins/_jetR); #if MODE!=1 if (part.E() < 1e-5 && part.E() > 1e-7) { double thpt(0.); if (thmommap.find(part.momentum()) == thmommap.end()) cout<<"Error: dummy not found in map.\n"; else { thpt = thmommap[part.momentum()].pT(); - if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]-=scalefac*thpt/jetpt; - if (cent<0.) difvecvac[bin]-=scalefac*thpt/jetpt; - for (size_t j=bin; j<_nbins; j++) { - if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]-=scalefac*thpt/jetpt; - if (cent<0.) intvecvac[j]-=scalefac*thpt/jetpt; + double m2 = thmommap[part.momentum()].mass2(); + if (m2 > 0.) { + if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]-=scalefac*thpt/jetpt; + if (cent<0.) difvecvac[bin]-=scalefac*thpt/jetpt; + for (size_t j=bin; j<_nbins; j++) { + if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]-=scalefac*thpt/jetpt; + if (cent<0.) intvecvac[j]-=scalefac*thpt/jetpt; + } + } + else { + if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]+=scalefac*thpt/jetpt; + if (cent<0.) difvecvac[bin]+=scalefac*thpt/jetpt; + for (size_t j=bin; j<_nbins; j++) { + if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]+=scalefac*thpt/jetpt; + if (cent<0.) intvecvac[j]+=scalefac*thpt/jetpt; + } } } } #endif -#if MODE==0 +#if (MODE==0 || MODE==4) else { if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]+=part.momentum().pT()/jetpt; if (cent<0.) difvecvac[bin]+=part.momentum().pT()/jetpt; for (size_t j=bin; j<_nbins; j++) { if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]+=part.momentum().pT()/jetpt; if (cent<0.) intvecvac[j]+=part.momentum().pT()/jetpt; } } #elif MODE==2 else if ((PID::charge(part.pdgId()) != 0) && (part.pT() > _trackptmin)) { //else if ((PID::charge(part.pdgId()) != 0)) { if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]+=part.momentum().pT()/jetpt; if (cent<0.) difvecvac[bin]+=part.momentum().pT()/jetpt; for (size_t j=bin; j<_nbins; j++) { if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]+=part.momentum().pT()/jetpt; if (cent<0.) intvecvac[j]+=part.momentum().pT()/jetpt; } } #endif if ((PID::charge(part.pdgId()) != 0) && (part.pT() > _trackptmin)) { double z,xi; z = part.momentum().p3().dot(jet.momentum().p3())/jet.momentum().p3().mod()/jet.momentum().p3().mod(); z*=jet.momentum().pT()/jetpt; xi = log(1/z); if (cent<0. || (cent>=0.0 && cent<0.1)) _h_ffxi->fill(xi,weight); if (cent<0.) _h_ffxiv->fill(xi,weight); if (cent<0. || (cent>=0.0 && cent<0.1)) _h_ffpt->fill(part.momentum().pT(),weight); if (cent<0.) _h_ffptv->fill(part.momentum().pT(),weight); } } #if MODE==1 foreach (PseudoJet pjet, pevent) { FourMomentum cellmom = Get4Mom(pjet); double dR(deltaR(cellmom,jet)); if (dR >= _jetR) continue; bin = int(dR*_nbins/_jetR); if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]+=cellmom.pT()/jetpt; if (cent<0.) difvecvac[bin]+=cellmom.pT()/jetpt; for (size_t j=bin; j<_nbins; j++) { if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]+=cellmom.pT()/jetpt; if (cent<0.) intvecvac[j]+=cellmom.pT()/jetpt; } } #endif for (size_t i=0; i < _nbins; ++i){ _p_difjetshape->fill(_delR_center[i], difvec[i]/_deltar, weight); _p_difjetshapev->fill(_delR_center[i], difvecvac[i]/_deltar, weight); + if (cent < 0.) _p_difjetshaper->fill(_delR_center[i], difvecvac[i]/_deltar, weight); + else _p_difjetshaper->fill(_delR_center[i], difvec[i]/_deltar, weight); _p_intjetshape->fill(_delR_center[i], intvec[i], weight); _p_intjetshapev->fill(_delR_center[i], intvecvac[i], weight); } } } /// Normalise histograms etc., after the run void finalize() { if (_Njets==0.) _Njets=1.; scale(_h_ffxi,1./_Njets); scale(_h_ffxiv,1./_Njets); scale(_h_ffpt,1./_Njets); scale(_h_ffptv,1./_Njets); #if MODE==2 normalize(_h_scalefac); #endif } private: double _Njets; double _jetR, _trackptmin; double _deltar; vector _delR_center; size_t _nbins; size_t _netabins, _nphibins; private: Profile1DPtr _p_intjetshape; Profile1DPtr _p_intjetshapev; Profile1DPtr _p_difjetshape; Profile1DPtr _p_difjetshapev; + Profile1DPtr _p_difjetshaper; Histo1DPtr _h_ffxi; Histo1DPtr _h_ffxiv; Histo1DPtr _h_ffpt; Histo1DPtr _h_ffptv; #if MODE==2 Histo1DPtr _h_scalefac; #endif }; // The hook for the plugin system #if MODE==0 DECLARE_RIVET_PLUGIN(CMS_2013_I1256590_4MOMSUB); #elif MODE==1 DECLARE_RIVET_PLUGIN(CMS_2013_I1256590_GRIDSUB2); #elif MODE==2 DECLARE_RIVET_PLUGIN(CMS_2013_I1256590_CHARGED); +#elif MODE==4 + DECLARE_RIVET_PLUGIN(CMS_2013_I1256590_QSUB); #endif } Index: trunk/code/CMS_GROOMEDMASS_CONSTSUB.plot =================================================================== --- trunk/code/CMS_GROOMEDMASS_CONSTSUB.plot (revision 0) +++ trunk/code/CMS_GROOMEDMASS_CONSTSUB.plot (revision 467) @@ -0,0 +1,38 @@ +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/* +LogY=0 +XLabel=$M_g/p_\perp^\text{jet}$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_1_140_160 +Title=$140\,\text{GeV} < p_\perp^\text{jet} < $140\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_1_160_180 +Title=$160\,\text{GeV} < p_\perp^\text{jet} < $180\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_1_180_200 +Title=$180\,\text{GeV} < p_\perp^\text{jet} < $200\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_1_200_300 +Title=$200\,\text{GeV} < p_\perp^\text{jet} < $300\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_2_140_160 +Title=$140\,\text{GeV} < p_\perp^\text{jet} < $140\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_2_160_180 +Title=$160\,\text{GeV} < p_\perp^\text{jet} < $180\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_2_180_200 +Title=$180\,\text{GeV} < p_\perp^\text{jet} < $200\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_CONSTSUB/Mgoverpt_2_200_300 +Title=$200\,\text{GeV} < p_\perp^\text{jet} < $300\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + + Index: trunk/code/CMS_2013_I1256590_CONSTSUB.cc =================================================================== --- trunk/code/CMS_2013_I1256590_CONSTSUB.cc (revision 0) +++ trunk/code/CMS_2013_I1256590_CONSTSUB.cc (revision 467) @@ -0,0 +1,476 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Tools/Logging.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/FastJets.hh" +#include "Rivet/Projections/ChargedFinalState.hh" +#include +#include +#include +/// @todo Include more projections as required, e.g. ChargedFinalState, FastJets, ZFinder... + +namespace Rivet { + +using namespace fastjet; + + struct FourMomComp { + bool operator() (const FourMomentum& p1, const FourMomentum& p2) const + {return p1.pT() 0.5) n=0.001; + if (cent > 0.3 && cent < 0.5) n=3.88; + if (cent > 0.1 && cent < 0.3) n=5.10; + if (cent > 0. && cent < 0.1) n=5.23; + n=5.23; + double sigma(sqrt(c*c+s*s/pttrue+n*n/(pttrue*pttrue))); + double r1(1.0*rand()/RAND_MAX), r2(1.0*rand()/RAND_MAX); + double fac(max(sqrt(-2.*log(r1))*cos(2.*M_PI*r2)*sigma+1.,0.)); + return fac*pttrue; + } + + + + PseudoJets ConstSubEvent(PseudoJets pjparts, PseudoJets pjghosts) { + // sort constituents into two vectors: particles and ghosts (subtraction momenta) + PseudoJets subevent; + FourMomentum totpartmom, totghostmom; + vector parts; + vector ghosts; + + foreach (PseudoJet pj, pjparts) { + FourMomentum partmom = Get4Mom(pj); + if (partmom.E() - fabs(partmom.pz()) > 0.) { + My4Mom part; + part.pt = partmom.pT(); + part.mdelta = sqrt(partmom.mass2() + sqr(partmom.pT())) - partmom.pT(); + part.phi = partmom.p3().phi(); + part.y = partmom.rapidity(); + //if (abs(part.y) > 5.) { + // cout<<"particle rapidity: "< 0.) { + My4Mom ghost; + ghost.pt = submom.pT(); + ghost.mdelta = sqrt(submom.mass2() + sqr(submom.pT())) - submom.pT(); + ghost.phi = submom.p3().phi(); + ghost.y = submom.rapidity(); + //cout<<"ghost rapidity: "< dists; + for (size_t i=0; i < parts.size(); ++i) { + for (size_t j=0; j < ghosts.size(); ++j) { + //cout<<"i, j = "<::iterator liter=dists.begin(); liter != dists.end(); ++liter) { + //cout<<"dealing with dist "<dR<dR > 0.3) break; + size_t pnum = liter->ipart; + size_t gnum = liter->ighost; + double ptp = parts[pnum].pt; + double ptg = ghosts[gnum].pt; + //cout<<"pts: "< ptg) { + parts[pnum].pt -= ptg; + ghosts[gnum].pt = 0.; + } + else { + ghosts[gnum].pt -= ptp; + parts[pnum].pt = 0.; + } + double mdp = parts[pnum].mdelta; + double mdg = ghosts[gnum].mdelta; + //cout<<"masses: "< mdg) { + parts[pnum].mdelta -= mdg; + ghosts[gnum].mdelta = 0.; + } + else { + ghosts[gnum].mdelta -= mdp; + parts[pnum].mdelta = 0.; + } + + } + + for (size_t i=0; i < parts.size(); ++i) { + if (parts[i].pt > 0.) { + FourMomentum mom((parts[i].pt+parts[i].mdelta)*cosh(parts[i].y), + parts[i].pt*cos(parts[i].phi), + parts[i].pt*sin(parts[i].phi), + (parts[i].pt+parts[i].mdelta)*sinh(parts[i].y)); + totpartmom += mom; + //cout<<"particle rapidity: "< 0.) { + FourMomentum mom((ghosts[i].pt+ghosts[i].mdelta)*cosh(ghosts[i].y), + ghosts[i].pt*cos(ghosts[i].phi), + ghosts[i].pt*sin(ghosts[i].phi), + (ghosts[i].pt+ghosts[i].mdelta)*sinh(ghosts[i].y)); + totghostmom += mom; + //cout<<"ghost rapidity: "<particles(); + else parts = jet.constituents(); + double scalefac; + double ntracks(0.), nparts(0.); + FourMomentum chmom(0.,0.,0.,0.), allmom(0.,0.,0.,0.); + foreach (Particle part, parts) { + if (deltaR(part,jet) < _jetR && part.momentum().E() > 1e-5) { + nparts += 1.; + allmom += part.momentum(); + if (PID::charge(part.pdgId()) != 0 && part.pT() > ptcut) { + ntracks += 1.; + chmom += part.momentum(); + } + } + } + scalefac = (nparts>0.?ntracks/nparts:1); + //scalefac = (allmom.pT()>0.?chmom.pT()/allmom.pT():1.); + return scalefac; + } + + +/// 4-momentum adapter + FourMomentum Get4Mom(PseudoJet pjet){ + FourMomentum mom(pjet.E(), pjet.px(), pjet.py(), pjet.pz()); + return mom; + } + + + + /// Perform the per-event analysis + void analyze(const Event& event) { + const double weight = event.weight(); + const double cent = (event.genEvent()->heavy_ion()?event.genEvent()->heavy_ion()->impact_parameter():-1.); + + + PseudoJets parts, ghosts; + + foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { + FourMomentum mom(p->momentum()); + if(p->status() == 1) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + if (mom.E() > 1e-6) parts.push_back(pjet); + } + else if (p->status() == 3) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + ghosts.push_back(pjet); + } + else if(p->status() == 4) { + PseudoJet pjet(mom.px(),mom.py(),mom.pz(),mom.E()); + parts.push_back(pjet); + } + } + + PseudoJets newevent; + newevent = ConstSubEvent(parts, ghosts); + + + JetDefinition jet_def(antikt_algorithm, 0.3); + ClusterSequence cs(newevent, jet_def); + PseudoJets pjets; + pjets = sorted_by_pt(cs.inclusive_jets(10.0)); + + + foreach (PseudoJet pj, pjets) { + FourMomentum jet(Get4Mom(pj)); + double ptreco, ptgen; + ptgen = pj.pt(); + ptreco = ptsmear(ptgen,cent); + + if (fabs(pj.eta()) < 0.3 || fabs(pj.eta()) > 2.0 || ptreco < 100.*GeV) continue; + + vector difvec(_nbins, 0.), difvecvac(_nbins, 0.); + vector intvec(_nbins, 0.), intvecvac(_nbins, 0.); + + _Njets+=weight; + double jetpt = ptreco; + double scalefac=1.; + int bin; + foreach (PseudoJet pjpart, newevent){ + FourMomentum part(Get4Mom(pjpart)); + double dR(deltaR(part,jet)); + if (dR >= _jetR) continue; + bin = int(dR*_nbins/_jetR); + if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]+=part.pT()/jetpt; + if (cent<0.) difvecvac[bin]+=part.pT()/jetpt; + for (size_t j=bin; j<_nbins; j++) { + if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]+=part.pT()/jetpt; + if (cent<0.) intvecvac[j]+=part.pT()/jetpt; + } + if (part.pT() > _trackptmin) { + double z,xi; + z = part.p3().dot(jet.p3())/jet.p3().mod()/jet.p3().mod(); + z*=jet.pT()/jetpt; + xi = log(1/z); + if (cent<0. || (cent>=0.0 && cent<0.1)) _h_ffxi->fill(xi,weight); + if (cent<0.) _h_ffxiv->fill(xi,weight); + if (cent<0. || (cent>=0.0 && cent<0.1)) _h_ffpt->fill(part.pT(),weight); + if (cent<0.) _h_ffptv->fill(part.pT(),weight); + } + } + + for (size_t i=0; i < _nbins; ++i){ + _p_difjetshape->fill(_delR_center[i], difvec[i]/_deltar, weight); + _p_difjetshapev->fill(_delR_center[i], difvecvac[i]/_deltar, weight); + if (cent < 0.) _p_difjetshaper->fill(_delR_center[i], difvecvac[i]/_deltar, weight); + else _p_difjetshaper->fill(_delR_center[i], difvec[i]/_deltar, weight); + _p_intjetshape->fill(_delR_center[i], intvec[i], weight); + _p_intjetshapev->fill(_delR_center[i], intvecvac[i], weight); + } + } + + + /*const Jets jets03 = applyProjection(event, "Jets03").jetsByPt(3.*GeV); + const FinalState fs = applyProjection(event, "FS"); + + + Jets jetssub; + vector thermom = extractThermalMomenta(event); + map thmommap = buildThMomMap2(event, &thermom); + if (thmommap.size() != thermom.size()) { + cout<<"Error: size of map does not match number of momenta: "< jetmap; + foreach (Jet jet, jetssub) { + double ptreco, ptgen; + ptgen = jet.momentum().pT(); + ptreco = ptsmear(ptgen,cent); + //ptreco = ptsmear(ptgen,0.05); + jetmap[jet] = ptreco; + } + + Jets jets; + foreach (Jet jet, jetssub) { + if (fabs(jet.momentum().eta()) > 0.3 && fabs(jet.momentum().eta()) < 2.0 && jetmap[jet] > 100.*GeV) + jets.push_back(jet); + } + + // if (jets.size() == 0) vetoEvent; + + foreach (Jet jet, jets) { + vector difvec(_nbins, 0.), difvecvac(_nbins, 0.); + vector intvec(_nbins, 0.), intvecvac(_nbins, 0.); + _Njets+=weight; + double jetpt = jetmap[jet]; + double scalefac=1.; + int bin; + foreach (Particle part, fs.particles()){ + double dR(deltaR(part,jet)); + if (dR >= _jetR) continue; + bin = int(dR*_nbins/_jetR); + if (part.E() == 1e-6) { + double thpt(0.); + if (thmommap.find(part.momentum()) == thmommap.end()) cout<<"Error: dummy not found in map.\n"; + else { + thpt = thmommap[part.momentum()].pT(); + double m2 = thmommap[part.momentum()].mass2(); + if (m2 > 0.) { + if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]-=scalefac*thpt/jetpt; + if (cent<0.) difvecvac[bin]-=scalefac*thpt/jetpt; + for (size_t j=bin; j<_nbins; j++) { + if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]-=scalefac*thpt/jetpt; + if (cent<0.) intvecvac[j]-=scalefac*thpt/jetpt; + } + } + else { + if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]+=scalefac*thpt/jetpt; + if (cent<0.) difvecvac[bin]+=scalefac*thpt/jetpt; + for (size_t j=bin; j<_nbins; j++) { + if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]+=scalefac*thpt/jetpt; + if (cent<0.) intvecvac[j]+=scalefac*thpt/jetpt; + } + } + } + } + else { + if (cent<0. || (cent>=0.0 && cent<0.1)) difvec[bin]+=part.momentum().pT()/jetpt; + if (cent<0.) difvecvac[bin]+=part.momentum().pT()/jetpt; + for (size_t j=bin; j<_nbins; j++) { + if (cent<0. || (cent>=0.0 && cent<0.1)) intvec[j]+=part.momentum().pT()/jetpt; + if (cent<0.) intvecvac[j]+=part.momentum().pT()/jetpt; + } + } + if ((PID::charge(part.pdgId()) != 0) && (part.pT() > _trackptmin)) { + double z,xi; + z = part.momentum().p3().dot(jet.momentum().p3())/jet.momentum().p3().mod()/jet.momentum().p3().mod(); + z*=jet.momentum().pT()/jetpt; + xi = log(1/z); + if (cent<0. || (cent>=0.0 && cent<0.1)) _h_ffxi->fill(xi,weight); + if (cent<0.) _h_ffxiv->fill(xi,weight); + if (cent<0. || (cent>=0.0 && cent<0.1)) _h_ffpt->fill(part.momentum().pT(),weight); + if (cent<0.) _h_ffptv->fill(part.momentum().pT(),weight); + } + } + for (size_t i=0; i < _nbins; ++i){ + _p_difjetshape->fill(_delR_center[i], difvec[i]/_deltar, weight); + _p_difjetshapev->fill(_delR_center[i], difvecvac[i]/_deltar, weight); + if (cent < 0.) _p_difjetshaper->fill(_delR_center[i], difvecvac[i]/_deltar, weight); + else _p_difjetshaper->fill(_delR_center[i], difvec[i]/_deltar, weight); + _p_intjetshape->fill(_delR_center[i], intvec[i], weight); + _p_intjetshapev->fill(_delR_center[i], intvecvac[i], weight); + } + }*/ + } + + + /// Normalise histograms etc., after the run + void finalize() { + + if (_Njets==0.) _Njets=1.; + scale(_h_ffxi,1./_Njets); + scale(_h_ffxiv,1./_Njets); + scale(_h_ffpt,1./_Njets); + scale(_h_ffptv,1./_Njets); + } + + + private: + + double _Njets; + double _jetR, _trackptmin; + double _deltar; + vector _delR_center; + size_t _nbins; + size_t _netabins, _nphibins; + + private: + + Profile1DPtr _p_intjetshape; + Profile1DPtr _p_intjetshapev; + Profile1DPtr _p_difjetshape; + Profile1DPtr _p_difjetshapev; + Profile1DPtr _p_difjetshaper; + Histo1DPtr _h_ffxi; + Histo1DPtr _h_ffxiv; + Histo1DPtr _h_ffpt; + Histo1DPtr _h_ffptv; + + }; + + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(CMS_2013_I1256590_CONSTSUB); + +} Index: trunk/code/CMS_GROOMEDMASS_4MOMSUB.plot =================================================================== --- trunk/code/CMS_GROOMEDMASS_4MOMSUB.plot (revision 0) +++ trunk/code/CMS_GROOMEDMASS_4MOMSUB.plot (revision 467) @@ -0,0 +1,38 @@ +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/* +LogY=0 +XLabel=$M_g/p_\perp^\text{jet}$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_1_140_160 +Title=$140\,\text{GeV} < p_\perp^\text{jet} < $140\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_1_160_180 +Title=$160\,\text{GeV} < p_\perp^\text{jet} < $180\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_1_180_200 +Title=$180\,\text{GeV} < p_\perp^\text{jet} < $200\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_1_200_300 +Title=$200\,\text{GeV} < p_\perp^\text{jet} < $300\,\text{GeV} \qquad z_\text{cut}=0.1\ \beta = 0$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_2_140_160 +Title=$140\,\text{GeV} < p_\perp^\text{jet} < $140\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_2_160_180 +Title=$160\,\text{GeV} < p_\perp^\text{jet} < $180\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_2_180_200 +Title=$180\,\text{GeV} < p_\perp^\text{jet} < $200\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + +# BEGIN PLOT /CMS_GROOMEDMASS_4MOMSUB/Mgoverpt_2_200_300 +Title=$200\,\text{GeV} < p_\perp^\text{jet} < $300\,\text{GeV} \qquad z_\text{cut}=0.5\ \beta = 1.5$ +# END PLOT + + Index: trunk/code/CMS_2013_I1256590.yoda =================================================================== --- trunk/code/CMS_2013_I1256590.yoda (revision 466) +++ trunk/code/CMS_2013_I1256590.yoda (revision 467) @@ -1,206 +1,206 @@ -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d01-x01-y03 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d01-x01-y03 Path=/REF/CMS_2013_I1256590/d01-x01-y03 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-02 2.500000e-02 2.500000e-02 9.996849e-01 1.534981e-02 1.496040e-02 7.500000e-02 2.500000e-02 2.500000e-02 9.914790e-01 4.180134e-02 4.529186e-02 1.250000e-01 2.500000e-02 2.500000e-02 9.037203e-01 4.181550e-02 3.772316e-02 1.750000e-01 2.500000e-02 2.500000e-02 1.016727e+00 5.312961e-02 5.671926e-02 2.250000e-01 2.500000e-02 2.500000e-02 1.273673e+00 1.439611e-01 1.400741e-01 2.750000e-01 2.500000e-02 2.500000e-02 1.382885e+00 1.060964e-01 1.097498e-01 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d01-x01-y02 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d01-x01-y02 Path=/REF/CMS_2013_I1256590/d01-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-02 2.500000e-02 2.500000e-02 1.266784e+01 0.000000e+00 0.000000e+00 7.500000e-02 2.500000e-02 2.500000e-02 4.565196e+00 0.000000e+00 0.000000e+00 1.250000e-01 2.500000e-02 2.500000e-02 1.470847e+00 0.000000e+00 0.000000e+00 1.750000e-01 2.500000e-02 2.500000e-02 6.714706e-01 0.000000e+00 0.000000e+00 2.250000e-01 2.500000e-02 2.500000e-02 3.694601e-01 0.000000e+00 0.000000e+00 2.750000e-01 2.500000e-02 2.500000e-02 2.007718e-01 0.000000e+00 0.000000e+00 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d01-x01-y01 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d01-x01-y01 Path=/REF/CMS_2013_I1256590/d01-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-02 2.500000e-02 2.500000e-02 1.266784e+01 0.000000e+00 0.000000e+00 7.500000e-02 2.500000e-02 2.500000e-02 4.565196e+00 0.000000e+00 0.000000e+00 1.250000e-01 2.500000e-02 2.500000e-02 1.298714e+00 0.000000e+00 0.000000e+00 1.750000e-01 2.500000e-02 2.500000e-02 6.883952e-01 0.000000e+00 0.000000e+00 2.250000e-01 2.500000e-02 2.500000e-02 4.738880e-01 0.000000e+00 0.000000e+00 2.750000e-01 2.500000e-02 2.500000e-02 2.774874e-01 0.000000e+00 0.000000e+00 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d03-x01-y01 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d03-x01-y01 Path=/REF/CMS_2013_I1256590/d03-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-01 2.500000e-01 2.500000e-01 4.589493e-02 2.640026e-02 2.607364e-02 7.500000e-01 2.500000e-01 2.500000e-01 3.302992e-01 1.372295e-01 1.308755e-01 1.250000e+00 2.500000e-01 2.500000e-01 7.553553e-01 2.064520e-01 2.113501e-01 1.750000e+00 2.500000e-01 2.500000e-01 1.219367e+00 1.948854e-01 1.904400e-01 2.250000e+00 2.500000e-01 2.500000e-01 1.560549e+00 1.710532e-01 1.668595e-01 2.750000e+00 2.500000e-01 2.500000e-01 1.804274e+00 2.437255e-01 2.817904e-01 3.250000e+00 2.500000e-01 2.500000e-01 2.086064e+00 4.080724e-01 3.258002e-01 3.750000e+00 2.500000e-01 2.500000e-01 2.447121e+00 4.787017e-01 4.655116e-01 4.250000e+00 2.500000e-01 2.500000e-01 2.556013e+00 3.771229e-01 3.566198e-01 4.750000e+00 2.500000e-01 2.500000e-01 1.273626e+00 4.731263e-01 4.790329e-01 5.250000e+00 2.500000e-01 2.500000e-01 1.213483e-01 5.143852e-02 1.293548e-01 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d03-x01-y03 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d03-x01-y03 Path=/REF/CMS_2013_I1256590/d03-x01-y03 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-01 2.500000e-01 2.500000e-01 9.511278e-01 2.506266e-01 2.631579e-01 7.500000e-01 2.500000e-01 2.500000e-01 9.949875e-01 1.817043e-01 1.879699e-01 1.250000e+00 2.500000e-01 2.500000e-01 8.508772e-01 1.190476e-01 1.190476e-01 1.750000e+00 2.500000e-01 2.500000e-01 8.132832e-01 1.002506e-01 9.398496e-02 2.250000e+00 2.500000e-01 2.500000e-01 7.380952e-01 8.145363e-02 8.145363e-02 2.750000e+00 2.500000e-01 2.500000e-01 7.506266e-01 8.145363e-02 8.145363e-02 3.250000e+00 2.500000e-01 2.500000e-01 8.508772e-01 8.771930e-02 8.771930e-02 3.750000e+00 2.500000e-01 2.500000e-01 1.095238e+00 1.190476e-01 1.127820e-01 4.250000e+00 2.500000e-01 2.500000e-01 1.508772e+00 1.879699e-01 1.879699e-01 4.750000e+00 2.500000e-01 2.500000e-01 1.765664e+00 3.320802e-01 3.320802e-01 5.250000e+00 2.500000e-01 2.500000e-01 2.223058e+00 6.892231e-01 6.892231e-01 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d03-x01-y02 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d03-x01-y02 Path=/REF/CMS_2013_I1256590/d03-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-01 2.500000e-01 2.500000e-01 4.793716e-02 0.000000e+00 0.000000e+00 7.500000e-01 2.500000e-01 2.500000e-01 3.302992e-01 0.000000e+00 0.000000e+00 1.250000e+00 2.500000e-01 2.500000e-01 8.990452e-01 0.000000e+00 0.000000e+00 1.750000e+00 2.500000e-01 2.500000e-01 1.494066e+00 0.000000e+00 0.000000e+00 2.250000e+00 2.500000e-01 2.500000e-01 2.147498e+00 0.000000e+00 0.000000e+00 2.750000e+00 2.500000e-01 2.500000e-01 2.447121e+00 0.000000e+00 0.000000e+00 3.250000e+00 2.500000e-01 2.500000e-01 2.447121e+00 0.000000e+00 0.000000e+00 3.750000e+00 2.500000e-01 2.500000e-01 2.243057e+00 0.000000e+00 0.000000e+00 4.250000e+00 2.500000e-01 2.500000e-01 1.702521e+00 0.000000e+00 0.000000e+00 4.750000e+00 2.500000e-01 2.500000e-01 7.231755e-01 0.000000e+00 0.000000e+00 5.250000e+00 2.500000e-01 2.500000e-01 5.462544e-02 0.000000e+00 0.000000e+00 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d04-x01-y02 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d04-x01-y02 Path=/REF/CMS_2013_I1256590/d04-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 1.097304e+00 9.730405e-02 9.730405e-02 9.843526e-01 0.000000e+00 0.000000e+00 1.310850e+00 1.162410e-01 1.162410e-01 1.032045e+00 0.000000e+00 0.000000e+00 1.565954e+00 1.388630e-01 1.388630e-01 1.032045e+00 0.000000e+00 0.000000e+00 1.870703e+00 1.658865e-01 1.658865e-01 9.689500e-01 0.000000e+00 0.000000e+00 2.234760e+00 1.981700e-01 1.981700e-01 9.537884e-01 0.000000e+00 0.000000e+00 2.669665e+00 2.367350e-01 2.367350e-01 8.407319e-01 0.000000e+00 0.000000e+00 3.189208e+00 2.828065e-01 2.828065e-01 7.648242e-01 0.000000e+00 0.000000e+00 3.809856e+00 3.378430e-01 3.378430e-01 6.132976e-01 0.000000e+00 0.000000e+00 4.551290e+00 4.035905e-01 4.035905e-01 5.238152e-01 0.000000e+00 0.000000e+00 5.437013e+00 4.821335e-01 4.821335e-01 4.473886e-01 0.000000e+00 0.000000e+00 6.495107e+00 5.759600e-01 5.759600e-01 3.644549e-01 0.000000e+00 0.000000e+00 7.759116e+00 6.880480e-01 6.880480e-01 2.922493e-01 0.000000e+00 0.000000e+00 9.269112e+00 8.219480e-01 8.219480e-01 2.418587e-01 0.000000e+00 0.000000e+00 1.107298e+01 9.819200e-01 9.819200e-01 1.764312e-01 0.000000e+00 0.000000e+00 1.322790e+01 1.173000e+00 1.173000e+00 1.370839e-01 0.000000e+00 0.000000e+00 1.580216e+01 1.401255e+00 1.401255e+00 1.015896e-01 0.000000e+00 0.000000e+00 1.887738e+01 1.673970e+00 1.673970e+00 6.957701e-02 0.000000e+00 0.000000e+00 2.255110e+01 1.999745e+00 1.999745e+00 4.917913e-02 0.000000e+00 0.000000e+00 2.693975e+01 2.388910e+00 2.388910e+00 3.368195e-02 0.000000e+00 0.000000e+00 3.218248e+01 2.853820e+00 2.853820e+00 2.306820e-02 0.000000e+00 0.000000e+00 3.844549e+01 3.409190e+00 3.409190e+00 1.287031e-02 0.000000e+00 0.000000e+00 4.592735e+01 4.072655e+00 4.072655e+00 6.532336e-03 0.000000e+00 0.000000e+00 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d04-x01-y03 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d04-x01-y03 Path=/REF/CMS_2013_I1256590/d04-x01-y03 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 1.097304e+00 9.730405e-02 9.730405e-02 1.050239e+00 5.598086e-01 5.645933e-01 1.310850e+00 1.162410e-01 1.162410e-01 6.339713e-01 4.354067e-01 4.354067e-01 1.565954e+00 1.388630e-01 1.388630e-01 6.818182e-01 3.636364e-01 3.684211e-01 1.870703e+00 1.658865e-01 1.658865e-01 4.043062e-01 2.918660e-01 2.918660e-01 2.234760e+00 1.981700e-01 1.981700e-01 2.224880e-01 2.344498e-01 2.392344e-01 2.669665e+00 2.367350e-01 2.367350e-01 5.502392e-02 1.818182e-01 1.866029e-01 3.189208e+00 2.828065e-01 2.828065e-01 -5.502392e-02 1.387560e-01 1.435407e-01 3.809856e+00 3.378430e-01 3.378430e-01 -6.937799e-02 1.052632e-01 1.004785e-01 4.551290e+00 4.035905e-01 4.035905e-01 -9.808612e-02 7.177033e-02 7.177033e-02 5.437013e+00 4.821335e-01 4.821335e-01 -1.028708e-01 5.263158e-02 5.263158e-02 6.495107e+00 5.759600e-01 5.759600e-01 -8.851675e-02 3.827751e-02 3.349282e-02 7.759116e+00 6.880480e-01 6.880480e-01 -8.373206e-02 2.870813e-02 2.392344e-02 9.269112e+00 8.219480e-01 8.219480e-01 -6.937799e-02 1.435407e-02 2.392344e-02 1.107298e+01 9.819200e-01 9.819200e-01 -4.066986e-02 1.913876e-02 1.435407e-02 1.322790e+01 1.173000e+00 1.173000e+00 -3.588517e-02 1.435407e-02 1.435407e-02 1.580216e+01 1.401255e+00 1.401255e+00 -2.631579e-02 9.569378e-03 1.435407e-02 1.887738e+01 1.673970e+00 1.673970e+00 -7.177033e-03 1.435407e-02 4.784689e-03 2.255110e+01 1.999745e+00 1.999745e+00 -7.177033e-03 9.569378e-03 4.784689e-03 2.693975e+01 2.388910e+00 2.388910e+00 -2.392344e-03 0.000000e+00 0.000000e+00 3.218248e+01 2.853820e+00 2.853820e+00 2.392344e-03 0.000000e+00 0.000000e+00 3.844549e+01 3.409190e+00 3.409190e+00 2.392344e-03 0.000000e+00 0.000000e+00 4.592735e+01 4.072655e+00 4.072655e+00 -2.392344e-03 0.000000e+00 0.000000e+00 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d04-x01-y01 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d04-x01-y01 Path=/REF/CMS_2013_I1256590/d04-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 1.097304e+00 9.730405e-02 9.730405e-02 2.033383e+00 4.782024e-01 4.627072e-01 1.310850e+00 1.162410e-01 1.162410e-01 1.656449e+00 3.489591e-01 3.451167e-01 1.565954e+00 1.388630e-01 1.388630e-01 1.709530e+00 2.947629e-01 3.238530e-01 1.870703e+00 1.658865e-01 1.658865e-01 1.370839e+00 2.541162e-01 2.341778e-01 2.234760e+00 1.981700e-01 1.981700e-01 1.189440e+00 2.204898e-01 1.813987e-01 2.669665e+00 2.367350e-01 2.367350e-01 8.954777e-01 1.774117e-01 1.696391e-01 3.189208e+00 2.828065e-01 2.828065e-01 7.068302e-01 1.400369e-01 1.339018e-01 3.809856e+00 3.378430e-01 3.378430e-01 5.491943e-01 1.088062e-01 1.040393e-01 4.551290e+00 4.035905e-01 4.035905e-01 4.267141e-01 8.454043e-02 8.890473e-02 5.437013e+00 4.821335e-01 4.821335e-01 3.421737e-01 6.779131e-02 5.845289e-02 6.495107e+00 5.759600e-01 5.759600e-01 2.700890e-01 4.656965e-02 4.613880e-02 7.759116e+00 6.880480e-01 6.880480e-01 2.065706e-01 3.561760e-02 3.528807e-02 9.269112e+00 8.219480e-01 8.219480e-01 1.709530e-01 2.947629e-02 2.607165e-02 1.107298e+01 9.819200e-01 9.819200e-01 1.328274e-01 2.115517e-02 2.269066e-02 1.322790e+01 1.173000e+00 1.173000e+00 1.015896e-01 1.883195e-02 1.924512e-02 1.580216e+01 1.401255e+00 1.401255e+00 7.648242e-02 1.611232e-02 1.593491e-02 1.887738e+01 1.673970e+00 1.673970e+00 5.667932e-02 1.332960e-02 1.073732e-02 2.255110e+01 1.999745e+00 1.999745e+00 4.006265e-02 1.083773e-02 1.149923e-02 2.693975e+01 2.388910e+00 2.388910e+00 2.831749e-02 8.615028e-03 9.295896e-03 3.218248e+01 2.853820e+00 2.853820e+00 1.909070e-02 6.620015e-03 6.266981e-03 3.844549e+01 3.409190e+00 3.409190e+00 1.170828e-02 4.750580e-03 4.856209e-03 4.592735e+01 4.072655e+00 4.072655e+00 6.636175e-03 2.933691e-03 3.207351e-03 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d02-x01-y01 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d02-x01-y01 Path=/REF/CMS_2013_I1256590/d02-x01-y01 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-02 2.500000e-02 2.500000e-02 6.375449e-01 0.000000e+00 0.000000e+00 7.500000e-02 2.500000e-02 2.500000e-02 8.643001e-01 0.000000e+00 0.000000e+00 1.250000e-01 2.500000e-02 2.500000e-02 9.278651e-01 0.000000e+00 0.000000e+00 1.750000e-01 2.500000e-02 2.500000e-02 9.612714e-01 0.000000e+00 0.000000e+00 2.250000e-01 2.500000e-02 2.500000e-02 9.858084e-01 0.000000e+00 0.000000e+00 2.750000e-01 2.500000e-02 2.500000e-02 9.979320e-01 0.000000e+00 0.000000e+00 -END YODA_SCATTER2D +# END YODA_SCATTER2D -BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d02-x01-y02 +# BEGIN YODA_SCATTER2D /REF/CMS_2013_I1256590/d02-x01-y02 Path=/REF/CMS_2013_I1256590/d02-x01-y02 Title= Type=Scatter2D # xval xerr- xerr+ yval yerr- yerr+ 2.500000e-02 2.500000e-02 2.500000e-02 6.357711e-01 0.000000e+00 0.000000e+00 7.500000e-02 2.500000e-02 2.500000e-02 8.678478e-01 0.000000e+00 0.000000e+00 1.250000e-01 2.500000e-02 2.500000e-02 9.366820e-01 0.000000e+00 0.000000e+00 1.750000e-01 2.500000e-02 2.500000e-02 9.700882e-01 0.000000e+00 0.000000e+00 2.250000e-01 2.500000e-02 2.500000e-02 9.892897e-01 0.000000e+00 0.000000e+00 2.750000e-01 2.500000e-02 2.500000e-02 9.978656e-01 0.000000e+00 0.000000e+00 -END YODA_SCATTER2D +# END YODA_SCATTER2D