File: GammaGammaAnalysis.cc

package info (click to toggle)
herwig%2B%2B 2.6.0-1
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 27,128 kB
  • ctags: 24,739
  • sloc: cpp: 188,949; fortran: 23,193; sh: 11,365; python: 5,069; ansic: 3,539; makefile: 1,865; perl: 2
file content (126 lines) | stat: -rw-r--r-- 4,547 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
// -*- C++ -*-
//
// GammaGammaAnalysis.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GammaGammaAnalysis class.
//

#include "GammaGammaAnalysis.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"

using namespace Herwig;
GammaGammaAnalysis::GammaGammaAnalysis() :
  _ptharder(0.,250.,100), _ptsofter(0.,250.,100), _ptpair(0.,250.,100), 
  _Eharder(0.,3000.,100), _Esofter(0.,3000.,100), _Epair(0.,6000.,100), 
  _rapharder(-12.,12.,120), _rapsofter(-12.,12.,120), _rappair(-12.,12.,120), 
  _phiharder(-Constants::pi,Constants::pi,50), 
  _phisofter(-Constants::pi,Constants::pi,50), 
  _deltaphi(-Constants::pi,Constants::pi,100), 
  _mpair(0,1000,100)  {}

void GammaGammaAnalysis::dofinish() {
  AnalysisHandler::dofinish();
  string fname = generator()->filename() + string("-") + name() + string(".top");
  ofstream outfile(fname.c_str());
  using namespace HistogramOptions;
  _ptharder.topdrawOutput(outfile,Frame|Ylog,"BLACK","pt harder");
  _ptsofter.topdrawOutput(outfile,Frame|Ylog,"BLACK","pt softer");
  _ptpair.topdrawOutput(outfile,Frame|Ylog,"BLACK","pt pair");
  _Eharder.topdrawOutput(outfile,Frame|Ylog,"BLACK","E harder");
  _Esofter.topdrawOutput(outfile,Frame|Ylog,"BLACK","E softer");
  _Epair.topdrawOutput(outfile,Frame|Ylog,"BLACK","E pair");
  _rapharder.topdrawOutput(outfile,Frame,"BLACK","y harder");
  _rapsofter.topdrawOutput(outfile,Frame,"BLACK","y softer");
  _rappair.topdrawOutput(outfile,Frame,"BLACK","y pair");
  _phiharder.topdrawOutput(outfile,Frame,"BLACK","phi harder");
  _phisofter.topdrawOutput(outfile,Frame,"BLACK","phi softer");
  _deltaphi.topdrawOutput(outfile,Frame,"BLACK","Delta phi");
  _mpair.topdrawOutput(outfile,Frame|Ylog,"BLACK","M pair");
  outfile.close();
}

namespace {
  inline Lorentz5Momentum getMomentum(tcPPtr particle) {
    return particle->momentum();
    //Lorentz5Momentum tmp = particle->children()[0]->next()->momentum();
    //tmp += particle->children()[1]->next()->momentum();
    //tmp.rescaleMass();
    //return tmp;
  } 
}


void GammaGammaAnalysis::analyze(tEventPtr event, long, int, int) {
  //  AnalysisHandler::analyze(event, ieve, loop, state);
  // Rotate to CMS, extract final state particles and call analyze(particles).
  // find the Z
  Lorentz5Momentum p1, p2, ppair;  
  bool foundphotons = false;
  set<tcPPtr> particles;
  event->selectFinalState(inserter(particles));

  for(set<tcPPtr>::const_iterator it = particles.begin(); 
      it != particles.end(); ++it) {
    if((**it).id()==ParticleID::gamma) {
      // find the two hardest photons in the event
      if( getMomentum(*it).perp() > p2.perp() ) {
	if (getMomentum(*it).perp() > p1.perp()) {
	  p2 = p1;
	  p1 = getMomentum(*it);
	} else {
	  p2 = getMomentum(*it);
	}
      }
    }
  }

  //  cerr << "E1 = " <<  p1.e()/GeV << ", E1 = " <<  p1.e()/GeV << "\n";
  if (p1.perp()/GeV > 0 && p2.perp()/GeV > 0) foundphotons = true;

  ppair = p1 + p2;
  if (foundphotons) {
    _ptharder += p1.perp()/GeV;
    _ptsofter += p2.perp()/GeV;
    _ptpair += ppair.perp()/GeV;
    _Eharder += p1.e()/GeV;
    _Esofter += p2.e()/GeV;
    _Epair += ppair.e()/GeV;
    _rapharder += p1.rapidity();
    _rapsofter += p2.rapidity();
    _rappair += ppair.rapidity();
    _phiharder += p1.phi();
    _phisofter += p2.phi();
    _deltaphi += (p2.vect()).deltaPhi(p1.vect());
    _mpair += ppair.m()/GeV;
  } else {
    cerr << "Analysis/GammaGammaAnalysis: Found no hard photon in event " 
	 << event->number()  << ".\n";
    generator()->log() << "Analysis/GammaGammaAnalysis: " 
		       << "Found no hard photon in event " 
		       << event->number()  << ".\n"
		       << *event;    
  }  
}

NoPIOClassDescription<GammaGammaAnalysis> GammaGammaAnalysis::initGammaGammaAnalysis;
// Definition of the static class description member.

void GammaGammaAnalysis::Init() {

  static ClassDocumentation<GammaGammaAnalysis> documentation
    ("There is no documentation for the GammaGammaAnalysis class");

}