File: CDFMidPointPlugin.cc

package info (click to toggle)
fastjet 3.4.0%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 9,552 kB
  • sloc: cpp: 78,628; python: 6,112; sh: 1,038; fortran: 673; makefile: 636; ansic: 161
file content (176 lines) | stat: -rw-r--r-- 6,636 bytes parent folder | download | duplicates (2)
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2005-2021, 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. They are described in the original FastJet paper,
//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
//  FastJet as part of work towards a scientific publication, please
//  quote the version you use and include a citation to the manual and
//  optionally also to hep-ph/0512210.
//
//  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/>.
//----------------------------------------------------------------------
//FJENDHEADER

#include "fastjet/CDFMidPointPlugin.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/Error.hh"
#include <sstream>

// CDF stuff
#include "MidPointAlgorithm.hh"
#include "PhysicsTower.hh"
#include "Cluster.hh"

FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

using namespace std;
using namespace cdf;

thread_safety_helpers::FirstTimeTrue CDFMidPointPlugin::_first_time;

string CDFMidPointPlugin::description () const {
  ostringstream desc;
  
  string sm_scale_string = "split-merge uses ";
  switch(_sm_scale) {
  case SM_pt:
    sm_scale_string += "pt";
    break;
  case SM_Et:
    sm_scale_string += "Et";
    break;
  case SM_mt:
    sm_scale_string += "mt";
    break;
  case SM_pttilde:
    sm_scale_string += "pttilde (scalar sum of pts)";
    break;
  default:
    ostringstream err;
    err << "Unrecognized split-merge scale choice = " << _sm_scale;
    throw Error(err.str());
  }

  
  if (cone_area_fraction() == 1) {
    desc << "CDF MidPoint jet algorithm, with " ;
  } else {
    desc << "CDF MidPoint+Searchcone jet algorithm, with ";
  }
  desc << "seed_threshold = "     << seed_threshold     () << ", "
       << "cone_radius = "        << cone_radius        () << ", "
       << "cone_area_fraction = " << cone_area_fraction () << ", " 
       << "max_pair_size = "      << max_pair_size      () << ", "
       << "max_iterations = "     << max_iterations     () << ", "
       << "overlap_threshold  = " << overlap_threshold  () << ", "
       << sm_scale_string ;

  return desc.str();
}


void CDFMidPointPlugin::run_clustering(ClusterSequence & clust_seq) const {
  // print a banner if we run this for the first time
  _print_banner(clust_seq.fastjet_banner_stream());
 
  // create the physics towers needed by the CDF code
  vector<PhysicsTower> towers;
  towers.reserve(clust_seq.jets().size());
  for (unsigned i = 0; i < clust_seq.jets().size(); i++) {
    LorentzVector fourvect(clust_seq.jets()[i].px(),
			   clust_seq.jets()[i].py(),
			   clust_seq.jets()[i].pz(),
			   clust_seq.jets()[i].E());
    PhysicsTower tower(fourvect);
    // misuse one of the indices for tracking, since the MidPoint
    // implementation doesn't seem to make use of these indices
    tower.calTower.iEta = i;
    towers.push_back(tower);
  }

  // prepare the CDF algorithm
  MidPointAlgorithm m(_seed_threshold,_cone_radius,_cone_area_fraction,
		      _max_pair_size,_max_iterations,_overlap_threshold,
                      MidPointAlgorithm::SplitMergeScale(_sm_scale));
    
  // run the CDF algorithm
  std::vector<Cluster> jets;
  m.run(towers,jets);


  // now transfer the jets back into our own structure -- we will
  // mimic the cone code with a sequential recombination sequence in
  // which the jets are built up by adding one particle at a time
  for(vector<Cluster>::const_iterator jetIter = jets.begin(); 
                                      jetIter != jets.end(); jetIter++) {
    const vector<PhysicsTower> & tower_list = jetIter->towerList;
    int jet_k = tower_list[0].calTower.iEta;
  
    int ntow = int(jetIter->towerList.size());
    for (int itow = 1; itow < ntow; itow++) {
      int jet_i = jet_k;
      // retrieve our misappropriated index for the jet
      int jet_j = tower_list[itow].calTower.iEta;
      // do a fake recombination step with dij=0
      double dij = 0.0;
      clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
    }
  
    // NB: put a sensible looking d_iB just to be nice...
    double d_iB = clust_seq.jets()[jet_k].perp2();
    clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
  }


  // following code is for testing only
  //cout << endl;
  //for(vector<Cluster>::const_iterator jetIter = jets.begin(); 
  //                                    jetIter != jets.end(); jetIter++) {
  //  cout << jetIter->fourVector.pt() << " " << jetIter->fourVector.y() << endl;
  //}
  //cout << "-----------------------------------------------------\n";
  //vector<PseudoJet> ourjets(clust_seq.inclusive_jets());
  //for (vector<PseudoJet>::const_reverse_iterator ourjet = ourjets.rbegin();
  //     ourjet != ourjets.rend(); ourjet++) {
  //  cout << ourjet->perp() << " " << ourjet->rap() << endl;
  //}
  //cout << endl;
}

// print a banner for reference to the 3rd-party code
void CDFMidPointPlugin::_print_banner(ostream *ostr) const{
  if (! _first_time()) return;

  // make sure the user has not set the banner stream to NULL
  if (!ostr) return;  

  (*ostr) << "#-------------------------------------------------------------------------" << endl;
  (*ostr) << "# You are running the CDF MidPoint plugin for FastJet                     " << endl;
  (*ostr) << "# This is based on an implementation provided by Joey Huston.             " << endl;
  (*ostr) << "# If you use this plugin, please cite                                     " << endl;
  (*ostr) << "#   G. C. Blazey et al., hep-ex/0005012.                                  " << endl;
  (*ostr) << "# in addition to the usual FastJet reference.                             " << endl;
  (*ostr) << "#-------------------------------------------------------------------------" << endl;

  // make sure we really have the output done.
  ostr->flush();
}

FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh