Page MenuHomeHEPForge
Files F10881269

No OneTemporary
Actions

View Options

diff --git a/src/CMSlrnsac.cc b/src/CMSlrnsac.cc
--- a/src/CMSlrnsac.cc
+++ b/src/CMSlrnsac.cc
@@ -1,379 +1,392 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the CMSlrnsac class.
//
#include "CMSlrnsac.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace ThePEG;
CMSlrnsac::CMSlrnsac(): doRweight(0) {}
CMSlrnsac::~CMSlrnsac() {}
void CMSlrnsac::analyze(const tPVector & particles, double weight) {
vector<int> nch(5);
vector< vector<double> > thisphi(5);
vector< vector<double> > thiseta(5);
vector<double> allphi;
vector<double> alleta;
vector<double> phi13;
vector<double> eta13;
int Nch04 = 0;
for ( int i = 0, N = particles.size(); i < N; ++i ) {
if ( particles[i]->data().iCharge() == 0 ) continue;
double eta =particles[i]->eta();
if ( abs(eta) >= etamax ) continue;
Energy pt = particles[i]->momentum().perp();
if ( pt < 0.1*GeV ) continue;
if ( pt >= 0.4*GeV ) ++Nch04;
int ipt = ptbin(pt);
double phi = particles[i]->momentum().phi();
thisphi[ipt].push_back(phi);
thiseta[ipt].push_back(eta);
allphi.push_back(phi);
alleta.push_back(eta);
if ( ipt == 1 || ipt == 2 ) {
phi13.push_back(phi);
eta13.push_back(eta);
}
}
if ( Nch04 < 2 ) return;
int inch = nchbin(Nch04);
if ( refNch04[inch] == 0 ) {
refNch04[inch] = Nch04;
refphi[inch] = thisphi;
refeta[inch] = thiseta;
refweight[inch] = weight;
return;
}
mult->fill(Nch04, weight);
- double wt = weight/double(Nch04*(Nch04 - 1));
- sumweight[inch] += weight;
- double rwt = weight/double(Nch04*refNch04[inch]);
+
+ double rweight = weight;
switch ( doRweight ) {
case 1:
- rwt = sqrt(weight*refweight[inch])/double(Nch04*refNch04[inch]);
- sumrefweight[inch] += sqrt(weight*refweight[inch]);
+ rweight = sqrt(weight*refweight[inch]);
break;
case 2:
- rwt = weight*refweight[inch]/double(Nch04*refNch04[inch]);
- sumrefweight[inch] += weight*refweight[inch];
+ rweight = weight*refweight[inch];
break;
- default:
- sumrefweight[inch] += weight;
}
+ sumweight[inch] += weight;
+ sumrefweight[inch] += rweight;
+
+ double wt = weight/double(Nch04*(Nch04 - 1));
+ double rwt = rweight/double(Nch04*refNch04[inch]);
// First do phi-projections.
for ( int ipt = 0; ipt < 4; ++ipt ) {
for ( int i = 0, N = thisphi[ipt].size(); i < N; ++i ) {
// Start with signal sample.
for ( int j = i + 1; j < N; ++ j ) {
double deta = abs(thiseta[ipt][i] - thiseta[ipt][j]);
if ( deta < 2.0 ) continue;
double dphi = abs(thisphi[ipt][i] - thisphi[ipt][j]);
if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
Sphi[inch][ipt]->fill(dphi, wt);
}
// Then background sample.
for ( int j = 0, M = refphi[inch][ipt].size(); j < M; ++j ) {
double deta = abs(thiseta[ipt][i] - refeta[inch][ipt][j]);
if ( deta < 2.0 ) continue;
double dphi = abs(thisphi[ipt][i] - refphi[inch][ipt][j]);
if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
Bphi[inch][ipt]->fill(dphi, rwt);
}
}
}
// Now do 2d histos. Begin with 1<pt<3 bins
vector<double> refallphi = refphi[inch][1];
refallphi.insert(refallphi.end(),
refphi[inch][2].begin(), refphi[inch][2].end());
vector<double> refalleta = refeta[inch][1];
refalleta.insert(refalleta.end(),
refeta[inch][2].begin(), refeta[inch][2].end());
-
+
+ if ( phi13.size() > 1 )
+ wt = weight/double(phi13.size()*(phi13.size() - 1));
+ if ( phi13.size() > 0 && refallphi.size() > 0 )
+ rwt = rweight/double(phi13.size()*refallphi.size());
+
for ( int i = 0, N = phi13.size(); i < N; ++i ) {
// First Signal.
for ( int j = i + 1; j < N; ++j ) {
double deta = abs(eta13[i] - eta13[j]);
double dphi = abs(phi13[i] - phi13[j]);
if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
SMB13->fill(deta, dphi, wt);
SMB13->fill(-deta, dphi, wt);
if ( inch == 3 ) {
SHN13->fill(deta, dphi, wt);
SHN13->fill(-deta, dphi, wt);
}
dphi = dphi <= 0.5*Constants::pi? -dphi: 2.0*Constants::pi - dphi;
SMB13->fill(deta, dphi, wt);
SMB13->fill(-deta, dphi, wt);
if ( inch == 3 ) {
SHN13->fill(deta, dphi, wt);
SHN13->fill(-deta, dphi, wt);
}
}
// Then background.
for ( int j = 0, M = refallphi.size(); j < M; ++j ) {
double deta = abs(eta13[i] - refalleta[j]);
double dphi = abs(phi13[i] - refallphi[j]);
if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
BMB13->fill(deta, dphi, rwt);
BMB13->fill(-deta, dphi, rwt);
if ( inch == 3 ) {
BHN13->fill(deta, dphi, rwt);
BHN13->fill(-deta, dphi, rwt);
}
dphi = dphi <= 0.5*Constants::pi? -dphi: 2.0*Constants::pi - dphi;
BMB13->fill(deta, dphi, rwt);
BMB13->fill(-deta, dphi, rwt);
if ( inch == 3 ) {
BHN13->fill(deta, dphi, rwt);
BHN13->fill(-deta, dphi, rwt);
}
}
}
// Continue 2D with 0.1<pt bins
refallphi.insert(refallphi.end(),
refphi[inch][0].begin(), refphi[inch][0].end());
refallphi.insert(refallphi.end(),
refphi[inch][3].begin(), refphi[inch][3].end());
refallphi.insert(refallphi.end(),
refphi[inch][4].begin(), refphi[inch][4].end());
refalleta.insert(refalleta.end(),
refeta[inch][0].begin(), refeta[inch][0].end());
refalleta.insert(refalleta.end(),
refeta[inch][3].begin(), refeta[inch][3].end());
refalleta.insert(refalleta.end(),
refeta[inch][4].begin(), refeta[inch][4].end());
+ if ( allphi.size() > 1 )
+ wt = weight/double(allphi.size()*(allphi.size() - 1));
+ if ( allphi.size() > 0 && refallphi.size() > 0 )
+ rwt = rweight/double(allphi.size()*refallphi.size());
+
for ( int i = 0, N = allphi.size(); i < N; ++i ) {
// First Signal.
for ( int j = i + 1; j < N; ++j ) {
double deta = abs(alleta[i] - alleta[j]);
double dphi = abs(allphi[i] - allphi[j]);
if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
SMB01->fill(deta, dphi, wt);
SMB01->fill(-deta, dphi, wt);
if ( inch == 3 ) {
SHN01->fill(deta, dphi, wt);
SHN01->fill(-deta, dphi, wt);
}
- dphi = dphi <= 0.5*Constants::pi? -dphi: 2.0*Constants::pi - dphi;
+ dphi = ( dphi <= 0.5*Constants::pi? -dphi: 2.0*Constants::pi - dphi );
SMB01->fill(deta, dphi, wt);
SMB01->fill(-deta, dphi, wt);
if ( inch == 3 ) {
SHN01->fill(deta, dphi, wt);
SHN01->fill(-deta, dphi, wt);
}
}
// Then background.
for ( int j = 0, M = refallphi.size(); j < M; ++j ) {
double deta = abs(alleta[i] - refalleta[j]);
double dphi = abs(allphi[i] - refallphi[j]);
if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
BMB01->fill(deta, dphi, rwt);
BMB01->fill(-deta, dphi, rwt);
if ( inch == 3 ) {
BHN01->fill(deta, dphi, rwt);
BHN01->fill(-deta, dphi, rwt);
}
dphi = dphi <= 0.5*Constants::pi? -dphi: 2.0*Constants::pi - dphi;
BMB01->fill(deta, dphi, rwt);
BMB01->fill(-deta, dphi, rwt);
if ( inch == 3 ) {
BHN01->fill(deta, dphi, rwt);
BHN01->fill(-deta, dphi, rwt);
}
}
}
refNch04[inch] = Nch04;
refphi[inch] = thisphi;
refeta[inch] = thiseta;
refweight[inch] = weight;
}
IBPtr CMSlrnsac::clone() const {
return new_ptr(*this);
}
IBPtr CMSlrnsac::fullclone() const {
return new_ptr(*this);
}
void CMSlrnsac::dofinish() {
AnalysisHandler::dofinish();
unitNormalize(mult);
double sumw = 0.0;
double sumrefw = 0.0;
for ( int inch = 0; inch < 4; ++inch ) {
sumw += sumweight[inch];
sumrefw += sumrefweight[inch];
for ( int ipt = 0; ipt < 4; ++ipt ) {
if ( sumweight[inch] > 0.0 && sumrefweight[inch] > 0.0 ) {
Sphi[inch][ipt]->scale(1.0/sumweight[inch]);
Bphi[inch][ipt]->scale(1.0/sumrefweight[inch]);
}
string indx;
indx += '0' + inch + 1;
indx += '0' + ipt + 1;
iHistogramFactory().divide
- ("/CMSlrnsac/Rphi" + indx, *Sphi[inch][ipt], *Bphi[inch][ipt]);
+ (dir + "/Rphi" + indx, *Sphi[inch][ipt], *Bphi[inch][ipt]);
}
}
if ( sumw > 0.0 && sumrefw > 0.0 ) {
SMB01->scale(1.0/sumw);
BMB01->scale(1.0/sumrefw);
SMB13->scale(1.0/sumw);
BMB13->scale(1.0/sumrefw);
}
if ( sumweight[3] > 0.0 && sumrefweight[3] > 0.0 ) {
SHN01->scale(1.0/sumweight[3]);
BHN01->scale(1.0/sumrefweight[3]);
SHN13->scale(1.0/sumweight[3]);
BHN13->scale(1.0/sumrefweight[3]);
}
- iHistogramFactory().divide("/CMSlrnsac/RMB01", *SMB01, *BMB01);
- iHistogramFactory().divide("/CMSlrnsac/RHN01", *SHN01, *BHN01);
- iHistogramFactory().divide("/CMSlrnsac/RMB13", *SMB13, *BMB13);
- iHistogramFactory().divide("/CMSlrnsac/RHN13", *SHN13, *BHN13);
+ iHistogramFactory().divide(dir + "/RMB01", *SMB01, *BMB01);
+ iHistogramFactory().divide(dir + "/RHN01", *SHN01, *BHN01);
+ iHistogramFactory().divide(dir + "/RMB13", *SMB13, *BMB13);
+ iHistogramFactory().divide(dir + "/RHN13", *SHN13, *BHN13);
}
void CMSlrnsac::doinitrun() {
AnalysisHandler::doinitrun();
etamax = 2.4;
histogramFactory().registerClient(this);
- histogramFactory().mkdirs("/CMSlrnsac");
+ dir = "/CMSlrnsac";
+ if ( doRweight == 1 ) dir = "/CMSlrnsac1";
+ if ( doRweight == 2 ) dir = "/CMSlrnsac2";
+ histogramFactory().mkdirs(dir);
Bphi = Sphi = vector< vector<tH1DPtr> >(4, vector<tH1DPtr>(4));
for ( int inch = 0; inch < 4; ++inch )
for ( int ipt = 0; ipt < 4; ++ ipt ) {
string indx;
indx += '0' + inch + 1;
indx += '0' + ipt + 1;
Sphi[inch][ipt] =
- histogramFactory().createHistogram1D("/CMSlrnsac/Sphi" + indx,
+ histogramFactory().createHistogram1D(dir + "/Sphi" + indx,
17, 0.0, Constants::pi);
Bphi[inch][ipt] =
- histogramFactory().createHistogram1D("/CMSlrnsac/Bphi" + indx,
+ histogramFactory().createHistogram1D(dir + "/Bphi" + indx,
17, 0.0, Constants::pi);
}
SMB01 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/SMB01",
+ (dir + "/SMB01",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
BMB01 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/BMB01",
+ (dir + "/BMB01",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
SMB13 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/SMB13",
+ (dir + "/SMB13",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
BMB13 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/BMB13",
+ (dir + "/BMB13",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
SHN01 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/SHN01",
+ (dir + "/SHN01",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
BHN01 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/BHN01",
+ (dir + "/BHN01",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
SHN13 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/SHN13",
+ (dir + "/SHN13",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
BHN13 =
histogramFactory().createHistogram2D
- ("/CMSlrnsac/BHN13",
+ (dir + "/BHN13",
33, -2.0*etamax, 2.0*etamax,
30, -0.5*Constants::pi,1.5*Constants::pi);
- mult = histogramFactory().createHistogram1D("/CMSlrnsac/Nch",
+ mult = histogramFactory().createHistogram1D(dir + "/Nch",
200, -0.5, 199.5);
refphi = vector< vector< vector<double> > >(4);
refeta = vector< vector< vector<double> > >(4);
refNch04 = vector<int>(4);
refweight = vector<double>(4);
sumrefweight = vector<double>(4);
sumweight = vector<double>(4);
}
void CMSlrnsac::persistentOutput(PersistentOStream & os) const {
os << doRweight;
}
void CMSlrnsac::persistentInput(PersistentIStream & is, int) {
is >> doRweight;
}
ClassDescription<CMSlrnsac> CMSlrnsac::initCMSlrnsac;
// Definition of the static class description member.
void CMSlrnsac::Init() {
static ClassDocumentation<CMSlrnsac> documentation
("There is no documentation for the CMSlrnsac class");
static Switch<CMSlrnsac,int> interfaceRweight
("Rweight",
"How to weight the background samples",
&CMSlrnsac::doRweight, 0, true, false);
static SwitchOption interfaceRweightSameAsSignal
(interfaceRweight,
"SameAsSignal",
"Fill with the same weight as the signal irrespectively of the weight "
"of the reference event.",
0);
static SwitchOption interfaceRweightSqrtOfWeights
(interfaceRweight,
"SqrtOfWeights",
"Fill with the square root of the weight of the signal event times the "
"weight of the reference event.",
1);
static SwitchOption interfaceRweightProdOfWeights
(interfaceRweight,
"ProdOfWeights",
"Fill with the weight of the signal event times the weight "
"of the reference event.",
2);
}
diff --git a/src/CMSlrnsac.h b/src/CMSlrnsac.h
--- a/src/CMSlrnsac.h
+++ b/src/CMSlrnsac.h
@@ -1,208 +1,209 @@
// -*- C++ -*-
#ifndef THEPEG_CMSlrnsac_H
#define THEPEG_CMSlrnsac_H
//
// This is the declaration of the CMSlrnsac class.
//
#include "ThePEG/Handlers/AnalysisHandler.h"
namespace ThePEG {
/**
* Here is the documentation of the CMSlrnsac class.
*
* @see \ref CMSlrnsacInterfaces "The interfaces"
* defined for CMSlrnsac.
*/
class CMSlrnsac: public AnalysisHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
CMSlrnsac();
/**
* The destructor.
*/
virtual ~CMSlrnsac();
//@}
public:
/** @name Virtual functions required by the AnalysisHandler class. */
//@{
/**
* Analyze the given vector of particles. The default version calls
* analyze(tPPtr) for each of the particles.
* @param particles the vector of pointers to particles to be analyzed
* @param weight the event weight
*/
virtual void analyze(const tPVector & particles, double weight);
//@}
/**
* Find the correct pt-bin
*/
static int ptbin(Energy pt) {
if ( pt < 1.0*GeV ) return 0;
if ( pt < 2.0*GeV ) return 1;
if ( pt < 3.0*GeV ) return 2;
if ( pt < 4.0*GeV ) return 3;
return 4;
}
/**
* Find the correct multiplicity bin.
*/
static int nchbin(int n) {
if ( n < 35 ) return 0;
if ( n < 90 ) return 1;
if ( n < 110 ) return 2;
return 3;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
vector< vector<tH1DPtr> > Sphi;
vector< vector<tH1DPtr> > Bphi;
tH2DPtr SMB01;
tH2DPtr BMB01;
tH2DPtr SMB13;
tH2DPtr BMB13;
tH2DPtr SHN01;
tH2DPtr BHN01;
tH2DPtr SHN13;
tH2DPtr BHN13;
vector< vector< vector<double> > > refphi;
vector< vector< vector<double> > > refeta;
vector<int> refNch04;
vector<double> refweight;
vector<double> sumrefweight;
vector<double> sumweight;
tH1DPtr mult;
double etamax;
int doRweight;
+ string dir;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<CMSlrnsac> initCMSlrnsac;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
CMSlrnsac & operator=(const CMSlrnsac &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of CMSlrnsac. */
template <>
struct BaseClassTrait<CMSlrnsac,1> {
/** Typedef of the first base class of CMSlrnsac. */
typedef AnalysisHandler NthBase;
};
/** This template specialization informs ThePEG about the name of
* the CMSlrnsac class and the shared object where it is defined. */
template <>
struct ClassTraits<CMSlrnsac>
: public ClassTraitsBase<CMSlrnsac> {
/** Return a platform-independent class name */
static string className() { return "ThePEG::CMSlrnsac"; }
/**
* The name of a file containing the dynamic library where the class
* CMSlrnsac is implemented. It may also include several, space-separated,
* libraries if the class CMSlrnsac depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "CMSlrnsac.so"; }
};
/** @endcond */
}
#endif /* THEPEG_CMSlrnsac_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:05 AM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4982926
Default Alt Text
(18 KB)

Event Timeline