All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
SatoruJetFinderProcessor.h
Go to the documentation of this file.
1 /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
2 /*
3 ** This file is part of the MarlinReco Project.
4 ** Forming part of the SubPackage: SatoruJetFinder.
5 **
6 ** For the latest version download from Web CVS:
7 ** http://www-zeuthen.desy.de/lc-cgi-bin/cvsweb.cgi/marlinreco/?cvsroot=MarlinReco
8 **
9 ** $Id$
10 **
11 **
12 */
13 
14 #ifndef SatoruJetFinderProcessor_h
15 #define SatoruJetFinderProcessor_h 1
16 
17 #include "marlin/Processor.h"
18 
19 #include "lcio.h"
20 #include "IMPL/LCCollectionVec.h"
21 
22 
23 using namespace lcio ;
24 
25 
26 struct SatoruPartonArray{
27  int NumberOfPartons;
28  float Momentum[1200];
29  int PointerParticleToJets[300];
30 };
31 struct SatoruJetsArray{
32  int NumberOfJets;
33  float Momentum[80];
34 };
35 
36 extern "C" {
37  void syjkrn_( const char* GlobalModus_,
38  // Primary jet finding
39  int &NJetRequested_,float &Threshold_,
40  int &PrimaryJetFindingMode_,float* YCut_,
41  // Reassosiation after first jet finding
42  int &MergingMode_,float &MergingThreshold_,
43  // Second Jet finding mode
44  int &SecondJetFindingMode_,
45  // input array--> array of
46  // PPartons(DimensionOfInputArray(4-6),NumberOfPartons)
47  int &NumberOfPartons_,
48  int &DimensionOfInputArray_,float *PPartons_,
49  // Output
50  int &MaximalNumberOfJets_,
51  int &NJetsFound_,int *PointerParicleToJet_,
52  int &DimensionOfOutputArray_,float *PJets_,
53  float &YMinus_,float &YPlus_,
54  // error flag
55  int &IError_,int GlobalModusLenght_)
56  ;}
57 
58 
59 namespace marlin
60 {
61 
62 
63  /**
64  * A universal jetfinder module developed by Satoru Yamashita
65  *
66  * This processor is a wrapper to the multi algorithm jet finder
67  * code written by Satoru Yamashita for the OPAL collaboratoin.
68  * For further details concerning this fortran code read the
69  * OPAL Technical Note TN579.
70  *
71  * To compile this processor, you need a FORTRAN77 compiler like g77
72  * and you have to link against the CERNLIB. If you use the MARLIN package
73  * mechanism you have to add the following two lines to the 'userlib.gmk'
74  * file in the MARLIN directory:
75  * USERLIBS += -lg2c
76  * USERLIBS += -L /opt/products/cernlib/pro/lib/ -lmathlib -lkernlib
77  *
78  * Steering file parameters:\n
79  * Allways to specify:\n
80  * \b InputCollection (e.g. energy flow object )of type reconstructed particle
81  * \b OutputCollection (e.g. jets )of type reconstructed particle
82  * \b Mode
83  *
84  * At the moment 5 modes are implemented:
85  *
86  * \b DurhamNjet:
87  * durham jetfinding with a fixed number of jets to be specify
88  * with
89  * NJetRequested
90  *
91  * \b DurhamYCut:
92  * durham jetfinding with a fixed ycut to be specify with
93  * YCut
94  *
95  * \b ConeBlanka
96  * cone jet finder with cone threshold energy of 0.7 GeV and
97  * R =0.2
98  *
99  * \b Saturo
100  * first durham jetfinding, afterwards reassignment of the
101  * objects to the jet axes found in this first iteration
102  * with the jade sceme
103  * number of jets to be specify with
104  * NJetRequested
105  *
106  * \b Manual everything to be set by hand:
107  * \verbatim
108  * GlobalMode (=MD below)
109  * NJetRequested (=NJETRQ)
110  * Threshold (= THRESH)
111  * PrimaryJetFindingMode (=IMODE)
112  * YCut (=YCUTS)
113  * MerginMode (=IEXAM)
114  * MergingThreshold (=FRAC)
115  * SecondJetFindingMode (=IEXAM4)
116  *
117  * Description
118  * Inputs:
119  * MD ; mode (see below)
120  * NJETRQ : fixed number of jets ; =0 variable -- use YCUT etc...
121  * THRESH : threshold (GeV) for the primary jet finding
122  * IMODE : Kernel Jet finding mode
123  * 1- 6 : YKERN
124  * IF(IMODE.EQ.1) THEN
125  * CM = 'JADE E0'
126  * ELSEIF(IMODE.EQ.2) THEN
127  * CM = 'JADE P '
128  * ELSEIF(IMODE.EQ.3) THEN
129  * CM = 'JADE P0'
130  * ELSEIF(IMODE.EQ.4) THEN
131  * CM = 'JADE E '
132  * ELSEIF(IMODE.EQ.5) THEN
133  * CM = 'DURHAM '
134  * ELSEIF(IMODE.EQ.6) THEN
135  * CM = 'GENEVA '
136  * ELSE
137  * WRITE(6,281) IMODE
138  * 281 FORMAT(/,' ### YKERN: IMODE =',I3,' INVALID; SET TO 1 ###')
139  * IMODE = 1
140  * CM = 'JADE E0'
141  * ENDIF
142  *
143  *
144  *
145  * 7 : cambridge
146  * 9-11 : LUCLUS
147  * 12 : cone
148  * YCUTS : (1) ycut or xcut or Cone Radius (2) EPS
149  * IEXAM : Particle re-association method 1-6
150  * EXAM(1) = THETA
151  * EXAM(2) = PART(4)*PJSTMP(4,J)-PART(1)*PJSTMP(1,J)
152  * + -PART(2)*PJSTMP(2,J)-PART(3)*PJSTMP(3,J)
153  * EXAM(3) = DURHAM
154  * EXAM(4) = E0JADE
155  * EXAM(5) = Geneva Scheme
156  * EXAM(6) = EJADE ! it makes too bad cross talk
157  *
158  * FRAC : core rate
159  * IEXAM4 : Jet merge method 1-8
160  * EXAM(1,N) = PL(4,I)*PL(4,J)
161  * EXAM(2,N) = PL(4,I)*PL(4,J)-PL(1,I)*PL(1,J)
162  * + -PL(2,I)*PL(2,J)-PL(3,I)*PL(3,J)
163  * EXAM(3,N) = DURHAM(I,J)
164  * EXAM(4,N) = E0JADE(I,J)
165  * EXAM(5,N) = EJADE(I,J)
166  * EXAM(6,N) = ANGLE(I,J)
167  * EXAM(7,N) = ECOS(I,J)
168  * EXAM(8,N) = PL(6,I)*PL(6,J)
169  * NPAR : Number of "particles"
170  * IDIM : first dimension of PPAR array
171  * PPAR : Particle Momenta
172  * MXJET : Maximum Number of Jets (buffer size for PJETS = MXJETS*JDIM)
173  * JDIM : first dimension of PJETS array
174  *
175  * Outputs:
176  * NJET : Number of Jets
177  * IASS : Jet assosiation for "particle"s
178  * PJETS : Jet momenta
179  * Y34 : YHI (sometimes not filled)
180  * Y45 : YLO (sometimes not filled)
181  * IERR : Error flag 0==All O.K.
182  * -9999: abnormal error
183  * -999: boundary error
184  * -99: input mismatch
185  * Others : error from various internal calls
186  * ==========================================================================
187  * MD
188  * There are 11 variations
189  * MD Name : Final NJ : 1st process : 2nd process : 3rd process
190  * 0a. Traditional: variable : fixed Ycuts : - : -
191  * 0b. Traditional: fixed Njets : fixed Njets : - : -
192  * 0c. Just Merge : fixed Njets : fixed Ycuts : do merge : -
193  * 1a. LMODE1 : variable : fixed Ycuts : do SYJETF1 : -
194  * 1b. LMODE1 : fixed final : fixed Ycuts : do SYJETF1 : do merge
195  * 1c. LMODE1 : fixed final : fixed Njets : do SYJETF1 : -
196  * 1d. LMODE1 : fixed final : fixed Ycuts : do merge : do SYJETF1
197  * 2a. LMODE2 : variable : fixed Ycuts : do SYJETF2 : -
198  * 2b. LMODE2 : fixed final : fixed Ycuts : do SYJETF2 : do merge
199  * 2c. LMODE2 : fixed final : fixed Njets : do SYJETF2 : -
200  * 2d. LMODE2 : fixed final : fixed Ycuts : do merge : do SYJETF2
201  * ===========================================================================
202  * MD & NJETRQ & IMODE & YCUTS & IEXAM & Frac & IEXAM4 & ISYMOD
203  * 0a & - & 1-12 & > 0 & - & - & - & 1
204  * 0b & > 0 & 1-12 & - & - & - & - & 2
205  * 0c & > 0 & 1-12 & > 0 & - & - & 1-8 & 3
206  * 1a & - & 1-12 & > 0 & 1-6 & 0.0-1.0 & - & 11
207  * 1b & > 0 & 1-12 & > 0 & 1-6 & 0.0-1.0 & 1-8 & 12
208  * 1c & > 0 & 1-12 & - & 1-6 & 0.0-1.0 & - & 13
209  * 1d & > 0 & 1-12 & > 0 & 1-6 & 0.0-1.0 & 1-8 & 14
210  * 2a & - & 1-12 & > 0 & 1-6 & 0.0-1.0 & - & 21
211  * 2b & > 0 & 1-12 & > 0 & 1-6 & 0.0-1.0 & 1-8 & 22
212  * 2c & > 0 & 1-12 & - & 1-6 & 0.0-1.0 & - & 23
213  * 2d & > 0 & 1-12 & > 0 & 1-6 & 0.0-1.0 & 1-8 & 24
214  * ===========================================================================
215  *
216  *
217  *
218  * \endverbatim
219  *
220  * @author Satoru Yamashita (original fortran code), Thorsten Kuhl, Jörgen Samson
221  * @version $Id$
222  */
223  class SatoruJetFinderProcessor : public Processor {
224 
225 
226  public:
227 
228  virtual Processor* newProcessor() { return new SatoruJetFinderProcessor ; }
229 
230 
232 
233  /** Called at the begin of the job before anything is read.
234  * Use to initialize the module, e.g. book histograms.
235  */
236  virtual void init() ;
237 
238  /** Called for every run.
239  */
240  virtual void processRunHeader( LCRunHeader* run ) ;
241 
242  /** Called for every event - the working horse.
243  */
244  virtual void processEvent( LCEvent * evt ) ;
245 
246 
247  /** Called after data processing for clean up.
248  */
249  virtual void end() ;
250 
251 
252 
253  private:
254 
255  std::string _inputCollection;
256  std::string _outputCollection;
257  std::string _successTag;
258  bool _writeTag;
259  std::string _jetFindingMode;
262  // int _nJetVal;
263  std::string _globalMode;
265  float _threshold;
267  float _yCutParam;
268  float _rConeParam;
270  float _yCut[2];
274  int _debug;
275  float _YMinus;
276  float _YPlus;
277  void putPartons(LCEvent * evt);
278  void callSatoru(LCEvent * evt);
279  void goSatoru(LCEvent * evt,LCCollection *JetsCol);
280  void getJets(LCEvent * evt,LCCollection *JetsCol);
281 
282  /* debug routines */
283  void writePartons();
284  void writeJets();
285  void writeParameters();
286 
287  // true constants for mode selection
288  struct MD
289  {
290  const static int JADE_E0 = 1;
291  const static int JADE_P = 2;
292  const static int JADE_P0 = 3;
293  const static int JADE_E = 4;
294  const static int DURHAM = 5;
295  const static int GENEVA = 6;
296  const static int CAMBRIDGE = 7;
297  const static int LUCLUS_1 = 8;
298  const static int LUCLUS_2 = 9;
299  const static int LUCLUS_3 = 10;
300  const static int LUCLUS_4 = 11;
301  const static int CONE = 12;
302  };
303 
304  } ;
305 
306 
307 
308 } // namespace marlin
309 #endif
310 
311 
312 
void syjkrn_(const char *GlobalModus_, int &NJetRequested_, float &Threshold_, int &PrimaryJetFindingMode_, float *YCut_, int &MergingMode_, float &MergingThreshold_, int &SecondJetFindingMode_, int &NumberOfPartons_, int &DimensionOfInputArray_, float *PPartons_, int &MaximalNumberOfJets_, int &NJetsFound_, int *PointerParicleToJet_, int &DimensionOfOutputArray_, float *PJets_, float &YMinus_, float &YPlus_, int &IError_, int GlobalModusLenght_)
A universal jetfinder module developed by Satoru Yamashita.