00001
00002 #include "Rivet/Rivet.hh"
00003 #include "Rivet/Analyses/HepEx9506012.hh"
00004 #include "Rivet/RivetCLHEP.hh"
00005 #include "Rivet/RivetAIDA.hh"
00006 using namespace CLHEP;
00007
00008
00009 namespace Rivet {
00010
00011 const double HepEx9506012::_xmin = -6.0;
00012 const double HepEx9506012::_xmax = 6.0;
00013
00014
00015 void HepEx9506012::init() {
00016 _hEtFlow = vector<AIDA::IHistogram1D *>(_nbin);
00017 _hEtFlowStat = vector<AIDA::IHistogram1D *>(_nbin);
00018 _nev = vector<double>(_nbin);
00019 for (int i = 0; i < _nbin; ++i) {
00020 string istr(1, char('1' + i));
00021 _hEtFlow[i] = bookHistogram1D(istr, "dEt/d[c] CMS bin=" + istr, _nb, _xmin, _xmax);
00022 _hEtFlowStat[i] = bookHistogram1D(istr, "stat dEt/d[c] CMS bin=1" + istr, _nb, _xmin, _xmax);
00023 }
00024 _hAvEt = bookHistogram1D("21tmp", "<Et> vs kin. bin", _nbin, 1.0, 10.0);
00025 _hAvX = bookHistogram1D("22tmp", "<x> vs kin. bin", _nbin, 1.0, 10.0);
00026 _hAvQ2 = bookHistogram1D("23tmp", "<Q2> vs kin. bin", _nbin, 1.0, 10.0);
00027 _hN = bookHistogram1D("24", "# events vs kin. bin", _nbin, 1.0, 10.0);
00028 }
00029
00030
00031 int HepEx9506012::getbin(const DISKinematics& dk) {
00032 const double GeV2 = GeV*GeV;
00033
00034 if ( dk.Q2() > 5.0*GeV2 && dk.Q2() <= 10.0*GeV2 ) {
00035 if ( dk.x() > 0.0001 && dk.x() <= 0.0002 )
00036 return 0;
00037 else if ( dk.x() > 0.0002 && dk.x() <= 0.0005 && dk.Q2() > 6.0*GeV2 )
00038 return 1;
00039 }
00040 else if ( dk.Q2() > 10.0*GeV2 && dk.Q2() <= 20.0*GeV2 ){
00041 if ( dk.x() > 0.0002 && dk.x() <= 0.0005 )
00042 return 2;
00043 else if ( dk.x() > 0.0005 && dk.x() <= 0.0008 )
00044 return 3;
00045 else if ( dk.x() > 0.0008 && dk.x() <= 0.0015 )
00046 return 4;
00047 else if ( dk.x() > 0.0015 && dk.x() <= 0.0040 )
00048 return 5;
00049 }
00050 else if ( dk.Q2() > 20.0*GeV2 && dk.Q2() <= 50.0*GeV2 ){
00051 if ( dk.x() > 0.0005 && dk.x() <= 0.0014 )
00052 return 6;
00053 else if ( dk.x() > 0.0014 && dk.x() <= 0.0030 )
00054 return 7;
00055 else if ( dk.x() > 0.0030 && dk.x() <= 0.0100 )
00056 return 8;
00057 }
00058 return -1;
00059 }
00060
00061
00062 void HepEx9506012::analyze(const Event& event) {
00063 const FinalStateHCM& fs = event.applyProjection(_fshcmproj);
00064 const DISKinematics& dk = event.applyProjection(_diskinproj);
00065 const CentralEtHCM y1 = event.applyProjection(_y1hcmproj);
00066
00067 int ibin = getbin(dk);
00068 if (ibin < 0) return;
00069
00070 const double weight = event.weight();
00071
00072 for ( int i = 0, N = fs.particles().size(); i < N; ++i ) {
00073 double rap = fs.particles()[i].getMomentum().rapidity();
00074 double et = fs.particles()[i].getMomentum().et();
00075 _hEtFlow[ibin]->fill(rap, weight*et/GeV);
00076 _hEtFlowStat[ibin]->fill(rap, weight*et/GeV);
00077 }
00078
00079 _nev[ibin] += weight;
00080 _hAvEt->fill(ibin + 1.5, weight*y1.sumEt()/GeV);
00081 _hAvX->fill(ibin + 1.5, weight*dk.x());
00082 _hAvQ2->fill(ibin + 1.5, weight*dk.Q2()/(GeV*GeV));
00083 _hN->fill(ibin + 1.5, weight);
00084 }
00085
00086
00087 void HepEx9506012::finalize() {
00088 for ( int ibin = 0; ibin < _nbin; ++ibin ) {
00089 _hEtFlow[ibin]->scale(1.0/(_nev[ibin]*double(_nb)/(_xmax-_xmin)));
00090 _hEtFlowStat[ibin]->scale(1.0/(_nev[ibin]*double(_nb)/(_xmax-_xmin)));
00091 }
00092 AIDA::IHistogram1D* h = 0;
00093 h = histogramFactory().divide("/HepEx9506012/21", *_hAvEt, *_hN);
00094 h->setTitle(_hAvEt->title());
00095 histogramFactory().destroy(_hAvEt);
00096 h = histogramFactory().divide("/HepEx9506012/22", *_hAvX, *_hN);
00097 h->setTitle(_hAvX->title());
00098 histogramFactory().destroy(_hAvX);
00099 h = histogramFactory().divide("/HepEx9506012/23", *_hAvQ2, *_hN);
00100 h->setTitle(_hAvQ2->title());
00101 histogramFactory().destroy(_hAvQ2);
00102 }
00103
00104 }