00001
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/RivetAIDA.hh"
00004 #include "Rivet/Tools/ParticleIdUtils.hh"
00005 #include "Rivet/Projections/Beam.hh"
00006 #include "Rivet/Projections/FinalState.hh"
00007 #include "Rivet/Projections/ChargedFinalState.hh"
00008 #include "Rivet/Projections/InitialQuarks.hh"
00009
00010 namespace Rivet {
00011
00012
00013
00014
00015 class OPAL_1998_S3780481 : public Analysis {
00016 public:
00017
00018
00019 OPAL_1998_S3780481() : Analysis("OPAL_1998_S3780481")
00020 {
00021 setBeams(ELECTRON, POSITRON);
00022
00023
00024 _weightedTotalPartNum = 0;
00025 _SumOfudsWeights = 0;
00026 _SumOfcWeights = 0;
00027 _SumOfbWeights = 0;
00028 }
00029
00030
00031
00032
00033
00034 void analyze(const Event& e) {
00035
00036 const FinalState& fs = applyProjection<FinalState>(e, "FS");
00037 const size_t numParticles = fs.particles().size();
00038
00039
00040 if (numParticles < 2) {
00041 getLog() << Log::DEBUG << "Failed ncharged cut" << endl;
00042 vetoEvent;
00043 }
00044 getLog() << Log::DEBUG << "Passed ncharged cut" << endl;
00045
00046
00047 const double weight = e.weight();
00048 _weightedTotalPartNum += numParticles * weight;
00049
00050
00051 const ParticlePair& beams = applyProjection<Beam>(e, "Beams").beams();
00052 const double meanBeamMom = ( beams.first.momentum().vector3().mod() +
00053 beams.second.momentum().vector3().mod() ) / 2.0;
00054 getLog() << Log::DEBUG << "Avg beam momentum = " << meanBeamMom << endl;
00055
00056 int flavour = 0;
00057 const InitialQuarks& iqf = applyProjection<InitialQuarks>(e, "IQF");
00058
00059
00060
00061 if (iqf.particles().size() == 2) {
00062 flavour = abs( iqf.particles().front().pdgId() );
00063 }
00064 else {
00065 map<int, double> quarkmap;
00066 foreach (const Particle& p, iqf.particles()) {
00067 if (quarkmap[p.pdgId()] < p.momentum().E()) {
00068 quarkmap[p.pdgId()] = p.momentum().E();
00069 }
00070 }
00071 double maxenergy = 0.;
00072 for (int i = 1; i <= 5; ++i) {
00073 if (quarkmap[i]+quarkmap[-i] > maxenergy) {
00074 flavour = i;
00075 }
00076 }
00077 }
00078
00079 switch (flavour) {
00080 case 1:
00081 case 2:
00082 case 3:
00083 _SumOfudsWeights += weight;
00084 break;
00085 case 4:
00086 _SumOfcWeights += weight;
00087 break;
00088 case 5:
00089 _SumOfbWeights += weight;
00090 break;
00091 }
00092
00093 foreach (const Particle& p, fs.particles()) {
00094 const double xp = p.momentum().vector3().mod()/meanBeamMom;
00095 const double logxp = -std::log(xp);
00096 _histXpall->fill(xp, weight);
00097 _histLogXpall->fill(logxp, weight);
00098 _histMultiChargedall->fill(_histMultiChargedall->binMean(0), weight);
00099 switch (flavour) {
00100
00101 case DQUARK:
00102 case UQUARK:
00103 case SQUARK:
00104 _histXpuds->fill(xp, weight);
00105 _histLogXpuds->fill(logxp, weight);
00106 _histMultiChargeduds->fill(_histMultiChargeduds->binMean(0), weight);
00107 break;
00108 case CQUARK:
00109 _histXpc->fill(xp, weight);
00110 _histLogXpc->fill(logxp, weight);
00111 _histMultiChargedc->fill(_histMultiChargedc->binMean(0), weight);
00112 break;
00113 case BQUARK:
00114 _histXpb->fill(xp, weight);
00115 _histLogXpb->fill(logxp, weight);
00116 _histMultiChargedb->fill(_histMultiChargedb->binMean(0), weight);
00117 break;
00118 }
00119 }
00120
00121 }
00122
00123
00124 void init() {
00125
00126 addProjection(Beam(), "Beams");
00127 addProjection(ChargedFinalState(), "FS");
00128 addProjection(InitialQuarks(), "IQF");
00129
00130
00131 _histXpuds = bookHistogram1D(1, 1, 1);
00132 _histXpc = bookHistogram1D(2, 1, 1);
00133 _histXpb = bookHistogram1D(3, 1, 1);
00134 _histXpall = bookHistogram1D(4, 1, 1);
00135 _histLogXpuds = bookHistogram1D(5, 1, 1);
00136 _histLogXpc = bookHistogram1D(6, 1, 1);
00137 _histLogXpb = bookHistogram1D(7, 1, 1);
00138 _histLogXpall = bookHistogram1D(8, 1, 1);
00139 _histMultiChargeduds = bookHistogram1D(9, 1, 1);
00140 _histMultiChargedc = bookHistogram1D(9, 1, 2);
00141 _histMultiChargedb = bookHistogram1D(9, 1, 3);
00142 _histMultiChargedall = bookHistogram1D(9, 1, 4);
00143 }
00144
00145
00146
00147 void finalize() {
00148 const double avgNumParts = _weightedTotalPartNum / sumOfWeights();
00149 normalize(_histXpuds , avgNumParts);
00150 normalize(_histXpc , avgNumParts);
00151 normalize(_histXpb , avgNumParts);
00152 normalize(_histXpall , avgNumParts);
00153 normalize(_histLogXpuds , avgNumParts);
00154 normalize(_histLogXpc , avgNumParts);
00155 normalize(_histLogXpb , avgNumParts);
00156 normalize(_histLogXpall , avgNumParts);
00157
00158 scale(_histMultiChargeduds, 1.0/_SumOfudsWeights);
00159 scale(_histMultiChargedc , 1.0/_SumOfcWeights);
00160 scale(_histMultiChargedb , 1.0/_SumOfbWeights);
00161 scale(_histMultiChargedall, 1.0/sumOfWeights());
00162 }
00163
00164
00165
00166
00167 private:
00168
00169
00170
00171
00172 double _weightedTotalPartNum;
00173
00174 double _SumOfudsWeights;
00175 double _SumOfcWeights;
00176 double _SumOfbWeights;
00177
00178 AIDA::IHistogram1D *_histXpuds;
00179 AIDA::IHistogram1D *_histXpc;
00180 AIDA::IHistogram1D *_histXpb;
00181 AIDA::IHistogram1D *_histXpall;
00182 AIDA::IHistogram1D *_histLogXpuds;
00183 AIDA::IHistogram1D *_histLogXpc;
00184 AIDA::IHistogram1D *_histLogXpb;
00185 AIDA::IHistogram1D *_histLogXpall;
00186 AIDA::IHistogram1D *_histMultiChargeduds;
00187 AIDA::IHistogram1D *_histMultiChargedc;
00188 AIDA::IHistogram1D *_histMultiChargedb;
00189 AIDA::IHistogram1D *_histMultiChargedall;
00190
00191
00192
00193 };
00194
00195
00196
00197
00198 AnalysisBuilder<OPAL_1998_S3780481> plugin_OPAL_1998_S3780481;
00199
00200 }