← Back to branch summary

~ubuntu-branches/ubuntu/saucy/rivet/saucy-proposed

« back to all changes in this revision

Viewing changes to src/Analyses/BABAR_2005_S6181155.cc

  • Committer: Package Import Robot
  • Author(s): Lifeng Sun
  • Date: 2013-05-07 09:24:27 UTC
  • mfrom: (1.2.1) (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130507092427-rpfijl4mn3y87ek7
Tags: 1.8.3-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Analyses/BABAR_2005_S6181155.cc
 
1
// -*- C++ -*-
 
2
#include <iostream>
 
3
#include "Rivet/Analysis.hh"
 
4
#include "Rivet/RivetAIDA.hh"
 
5
#include "Rivet/Tools/ParticleIdUtils.hh"
 
6
#include "Rivet/Projections/Beam.hh"
 
7
#include "Rivet/Projections/UnstableFinalState.hh"
 
8
 
 
9
namespace Rivet {
 
10
 
 
11
 
 
12
  /// @brief BABAR Xi_c baryons from fragmentation
 
13
  /// @author Peter Richardson
 
14
  class BABAR_2005_S6181155 : public Analysis {
 
15
  public:
 
16
 
 
17
    BABAR_2005_S6181155()
 
18
      : Analysis("BABAR_2005_S6181155")
 
19
    { }
 
20
 
 
21
 
 
22
    void analyze(const Event& e) {
 
23
      const double weight = e.weight();
 
24
 
 
25
      // Loop through unstable FS particles and look for charmed mesons/baryons
 
26
      const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS");
 
27
 
 
28
      const Beam beamproj = applyProjection<Beam>(e, "Beams");
 
29
      const ParticlePair& beams = beamproj.beams();
 
30
      FourMomentum mom_tot = beams.first.momentum() + beams.second.momentum();
 
31
      LorentzTransform cms_boost(-mom_tot.boostVector());
 
32
      const double s = sqr(beamproj.sqrtS());
 
33
 
 
34
      const bool onresonance = fuzzyEquals(beamproj.sqrtS(), 10.58, 2E-3);
 
35
 
 
36
      foreach (const Particle& p, ufs.particles()) {
 
37
        // 3-momentum in CMS frame
 
38
        const double mom = cms_boost.transform(p.momentum()).vector3().mod();
 
39
        // only looking at Xi_c^0
 
40
        if(abs(p.pdgId()) != 4132 ) continue;
 
41
        if (onresonance) {
 
42
          _histOnResonanceA_norm->fill(mom,weight);
 
43
          _histOnResonanceB_norm->fill(mom,weight);
 
44
        }
 
45
        else {
 
46
          _histOffResonance_norm->fill(mom,s/sqr(10.58)*weight);
 
47
        }
 
48
        MSG_DEBUG("mom = " << mom);
 
49
        // off-resonance cross section
 
50
        if(checkDecay(p.genParticle())) {
 
51
          if (onresonance) {
 
52
            _histOnResonanceA->fill(mom,weight);
 
53
            _histOnResonanceB->fill(mom,weight);
 
54
          }
 
55
          else {
 
56
            _histOffResonance->fill(mom,s/sqr(10.58)*weight);
 
57
            _sigma->fill(10.6,weight);
 
58
          }
 
59
        }
 
60
      }
 
61
    } // analyze
 
62
 
 
63
    void finalize() {
 
64
      scale(_histOnResonanceA, crossSection()/femtobarn/sumOfWeights());
 
65
      scale(_histOnResonanceB, crossSection()/femtobarn/sumOfWeights());
 
66
      scale(_histOffResonance, crossSection()/femtobarn/sumOfWeights());
 
67
      scale(_sigma           , crossSection()/femtobarn/sumOfWeights());
 
68
      normalize(_histOnResonanceA_norm);
 
69
      normalize(_histOnResonanceB_norm);
 
70
      normalize(_histOffResonance_norm);
 
71
    } // finalize
 
72
 
 
73
 
 
74
    void init() {
 
75
      addProjection(Beam(), "Beams");
 
76
      addProjection(UnstableFinalState(), "UFS");
 
77
 
 
78
      _histOnResonanceA = bookHistogram1D(1,1,1);
 
79
      _histOnResonanceB = bookHistogram1D(2,1,1);
 
80
      _histOffResonance = bookHistogram1D(2,1,2);
 
81
      _sigma            = bookHistogram1D(3,1,1);
 
82
      _histOnResonanceA_norm = bookHistogram1D(4,1,1);
 
83
      _histOnResonanceB_norm = bookHistogram1D(5,1,1);
 
84
      _histOffResonance_norm = bookHistogram1D(5,1,2);
 
85
 
 
86
    } // init
 
87
 
 
88
  private:
 
89
 
 
90
    //@{
 
91
    /// Histograms
 
92
    AIDA::IHistogram1D *_histOnResonanceA;
 
93
    AIDA::IHistogram1D *_histOnResonanceB;
 
94
    AIDA::IHistogram1D *_histOffResonance;
 
95
    AIDA::IHistogram1D *_sigma           ;
 
96
    AIDA::IHistogram1D *_histOnResonanceA_norm;
 
97
    AIDA::IHistogram1D *_histOnResonanceB_norm;
 
98
    AIDA::IHistogram1D *_histOffResonance_norm;
 
99
    //@}
 
100
 
 
101
    bool checkDecay(const GenParticle & p) {
 
102
      unsigned int nstable=0,npip=0,npim=0;
 
103
      unsigned int nXim=0,nXip=0;
 
104
      findDecayProducts(p,nstable,npip,npim,
 
105
                        nXip,nXim);
 
106
      int id = p.pdg_id();
 
107
      // Xi_c
 
108
      if(id==4132) {
 
109
        if(nstable==2&&nXim==1&&npip==1) return true;
 
110
      }
 
111
      else if(id==-4132) {
 
112
        if(nstable==2&&nXip==1&&npim==1) return true;
 
113
      }
 
114
      return false;
 
115
    }
 
116
 
 
117
    void findDecayProducts(const GenParticle & p,
 
118
                           unsigned int & nstable, 
 
119
                           unsigned int & npip, unsigned int & npim,
 
120
                           unsigned int & nXip, unsigned int & nXim) {
 
121
      const GenVertex* dv = p.end_vertex();
 
122
      for (GenVertex::particles_out_const_iterator
 
123
             pp = dv->particles_out_const_begin();
 
124
           pp != dv->particles_out_const_end(); ++pp) {
 
125
        int id = (*pp)->pdg_id();
 
126
        if(id==3312) { 
 
127
          ++nXim;
 
128
          ++nstable;
 
129
        }
 
130
        else if(id==-3312) { 
 
131
          ++nXip;
 
132
          ++nstable;
 
133
        }
 
134
        else if(id==111||id==221)
 
135
          ++nstable;
 
136
        else if((*pp)->end_vertex())
 
137
          findDecayProducts(**pp,nstable,npip,npim,nXip,nXim);
 
138
        else {
 
139
          if(id!=22) ++nstable;
 
140
          if     (id ==   211) ++npip;
 
141
          else if(id ==  -211) ++npim;
 
142
        }
 
143
      }
 
144
    }
 
145
  };
 
146
 
 
147
  // The hook for the plugin system
 
148
  DECLARE_RIVET_PLUGIN(BABAR_2005_S6181155);
 
149
 
 
150
}