File: 06-area.cc

package info (click to toggle)
fastjet 3.0.6%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 9,404 kB
  • ctags: 3,770
  • sloc: cpp: 21,498; sh: 10,546; fortran: 673; makefile: 527; ansic: 131
file content (136 lines) | stat: -rw-r--r-- 5,397 bytes parent folder | download | duplicates (3)
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
127
128
129
130
131
132
133
134
135
136
//----------------------------------------------------------------------
/// \file
/// \page Example06 06 - using jet areas
///
/// fastjet example program for jet areas
/// It mostly illustrates the usage of the 
/// fastjet::AreaDefinition and fastjet::ClusterSequenceArea classes
///
/// run it with    : ./06-area < data/single-event.dat
///
/// Source code: 06-area.cc
//----------------------------------------------------------------------

//STARTHEADER
// $Id: 06-area.cc 2684 2011-11-14 07:41:44Z soyez $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
//  FastJet is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  The algorithms that underlie FastJet have required considerable
//  development and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  FastJet is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER

#include "fastjet/ClusterSequenceArea.hh"  // use this instead of the "usual" ClusterSequence to get area support
#include <iostream> // needed for io
#include <cstdio>   // needed for io

using namespace std;

/// an example program showing how to use fastjet
int main(){
  
  // read in input particles
  //----------------------------------------------------------
  vector<fastjet::PseudoJet> input_particles;
  
  double px, py , pz, E;
  while (cin >> px >> py >> pz >> E) {
    // create a fastjet::PseudoJet with these components and put it onto
    // back of the input_particles vector
    input_particles.push_back(fastjet::PseudoJet(px,py,pz,E)); 
  }
  

  // create a jet definition: 
  // a jet algorithm with a given radius parameter
  //----------------------------------------------------------
  double R = 0.6;
  fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R);


  // Now we also need an AreaDefinition to define the properties of the 
  // area we want
  //
  // This is made of 2 building blocks:
  //  - the area type:
  //    passive, active, active with explicit ghosts, or Voronoi area
  //  - the specifications:
  //    a VoronoiSpec or a GhostedAreaSpec for the 3 ghost-bases ones
  // 
  //---------------------------------------------------------- For
  // GhostedAreaSpec (as below), the minimal info you have to provide
  // is up to what rapidity ghosts are placed. 
  // Other commonm parameters (that mostly have an impact on the
  // precision on the area) include the number of repetitions
  // (i.e. the number of different sets of ghosts that are used) and
  // the ghost density (controlled through the ghost_area).
  // Other, more exotic, parameters (not shown here) control how ghosts
  // are placed.
  //
  // The ghost rapidity interval should be large enough to cover the
  // jets for which you want to calculate. E.g. if you want to
  // calculate the area of jets up to |y|=4, you need to put ghosts up
  // to at least 4+R (or, optionally, up to the largest particle
  // rapidity if this is smaller).
  double maxrap = 5.0;
  unsigned int n_repeat = 3; // default is 1
  double ghost_area = 0.01; // this is the default
  fastjet::GhostedAreaSpec area_spec(maxrap, n_repeat, ghost_area);

  fastjet::AreaDefinition area_def(fastjet::active_area, area_spec);

  // run the jet clustering with the above jet and area definitions
  //
  // The only change is the usage of a ClusterSequenceArea rather than
  //a ClusterSequence
  //----------------------------------------------------------
  fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def);


  // get the resulting jets ordered in pt
  //----------------------------------------------------------
  double ptmin = 5.0;
  vector<fastjet::PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(ptmin));


  // tell the user what was done
  //  - the description of the algorithm and area used
  //  - extract the inclusive jets with pt > 5 GeV
  //    show the output as 
  //      {index, rap, phi, pt, number of constituents}
  //----------------------------------------------------------
  cout << endl;
  cout << "Ran " << jet_def.description() << endl;
  cout << "Area: " << area_def.description() << endl << endl;

  // label the columns
  printf("%5s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area", "area error");
 
  // print out the details for each jet
  for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
    printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
	   inclusive_jets[i].rap(), inclusive_jets[i].phi(), inclusive_jets[i].perp(),
	   inclusive_jets[i].area(), inclusive_jets[i].area_error());
  }

  return 0;
}