source: HiSusy/trunk/Pythia8/pythia8170/examples/main31.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 20.4 KB
Line 
1// main31.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2012 Richard Corke, Torbjorn Sjostrand.
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6#include "Pythia.h"
7using namespace Pythia8; 
8
9//==========================================================================
10
11// Use userhooks to veto PYTHIA emissions above the POWHEG scale.
12
13class PowhegHooks : public UserHooks {
14
15public: 
16
17  // Constructor and destructor.
18   PowhegHooks(int nFinalIn, int vetoModeIn, int vetoCountIn,
19               int pThardModeIn, int pTemtModeIn, int emittedModeIn,
20               int pTdefModeIn, int MPIvetoModeIn) : nFinal(nFinalIn),
21               vetoMode(vetoModeIn), vetoCount(vetoCountIn),
22               pThardMode(pThardModeIn), pTemtMode(pTemtModeIn),
23               emittedMode(emittedModeIn), pTdefMode(pTdefModeIn),
24               MPIvetoMode(MPIvetoModeIn) {};
25  ~PowhegHooks() {}
26
27//--------------------------------------------------------------------------
28
29  // Routines to calculate the pT (according to pTdefMode) in a splitting:
30  //   ISR: i (radiator after)  -> j (emitted after) k (radiator before)
31  //   FSR: i (radiator before) -> j (emitted after) k (radiator after)
32  // For the Pythia pT definition, a recoiler (after) must be specified.
33
34  // Compute the Pythia pT separation. Based on pTLund function in History.cc
35  double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch,
36                  int RecAfterBranch, bool FSR) {
37
38    // Convenient shorthands for later
39    Vec4 radVec = e[RadAfterBranch].p();
40    Vec4 emtVec = e[EmtAfterBranch].p();
41    Vec4 recVec = e[RecAfterBranch].p();
42    int  radID  = e[RadAfterBranch].id();
43
44    // Calculate virtuality of splitting
45    double sign = (FSR) ? 1. : -1.;
46    Vec4 Q(radVec + sign * emtVec); 
47    double Qsq = sign * Q.m2Calc();
48
49    // Mass term of radiator
50    double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
51                   pow2(particleDataPtr->m0(radID)) : 0.;
52
53    // z values for FSR and ISR
54    double z, pTnow;
55    if (FSR) {
56      // Construct 2 -> 3 variables
57      Vec4 sum = radVec + recVec + emtVec;
58      double m2Dip = sum.m2Calc();
59      double x1 = 2. * (sum * radVec) / m2Dip;
60      double x3 = 2. * (sum * emtVec) / m2Dip;
61      z     = x1 / (x1 + x3);
62      pTnow = z * (1. - z);
63
64    } else {
65      // Construct dipoles before/after splitting
66      Vec4 qBR(radVec - emtVec + recVec);
67      Vec4 qAR(radVec + recVec);
68      z     = qBR.m2Calc() / qAR.m2Calc();
69      pTnow = (1. - z);
70    }
71
72    // Virtuality with correct sign
73    pTnow *= (Qsq - sign * m2Rad);
74
75    // Can get negative pT for massive splittings
76    if (pTnow < 0.) {
77      cout << "Warning: pTpythia was negative" << endl;
78      return -1.;
79    }
80
81#ifdef DBGOUTPUT
82    cout << "pTpythia: rad = " << RadAfterBranch << ", emt = "
83         << EmtAfterBranch << ", rec = " << RecAfterBranch
84         << ", pTnow = " << sqrt(pTnow) << endl;
85#endif
86
87    // Return pT
88    return sqrt(pTnow);
89  }
90
91  // Compute the POWHEG pT separation between i and j
92  double pTpowheg(const Event &e, int i, int j, bool FSR) {
93
94    // pT value for FSR and ISR
95    double pTnow = 0.;
96    if (FSR) {
97      // POWHEG d_ij (in CM frame). Note that the incoming beams have not
98      // been updated in the parton systems pointer yet (i.e. prior to any
99      // potential recoil).
100      int iInA = partonSystemsPtr->getInA(0);
101      int iInB = partonSystemsPtr->getInB(0);
102      double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
103                       ( e[iInA].e()  + e[iInB].e()  );
104      Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
105      iVecBst.bst(0., 0., betaZ);
106      jVecBst.bst(0., 0., betaZ);
107      pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
108                    iVecBst.e() * jVecBst.e() /
109                    pow2(iVecBst.e() + jVecBst.e()) );
110
111    } else {
112      // POWHEG pT_ISR is just kinematic pT
113      pTnow = e[j].pT();
114    }
115
116    // Check result
117    if (pTnow < 0.) {
118      cout << "Warning: pTpowheg was negative" << endl;
119      return -1.;
120    }
121
122#ifdef DBGOUTPUT
123    cout << "pTpowheg: i = " << i << ", j = " << j
124         << ", pTnow = " << pTnow << endl;
125#endif
126
127    return pTnow;
128  }
129
130  // Calculate pT for a splitting based on pTdefMode.
131  // If j is -1, all final-state partons are tried.
132  // If i, k, r and xSR are -1, then all incoming and outgoing
133  // partons are tried.
134  // xSR set to 0 means ISR, while xSR set to 1 means FSR
135  double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) {
136
137    // Loop over ISR and FSR if necessary
138    double pTemt = -1., pTnow;
139    int xSR1 = (xSRin == -1) ? 0 : xSRin;
140    int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
141    for (int xSR = xSR1; xSR < xSR2; xSR++) {
142      // FSR flag
143      bool FSR = (xSR == 0) ? false : true;
144
145      // If all necessary arguments have been given, then directly calculate.
146      // POWHEG ISR and FSR, need i and j.
147      if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
148        pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
149
150      // Pythia ISR, need i, j and r.
151      } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
152        pTemt = pTpythia(e, i, j, r, FSR);
153
154      // Pythia FSR, need k, j and r.
155      } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
156        pTemt = pTpythia(e, k, j, r, FSR);
157
158      // Otherwise need to try all possible combintations.
159      } else {
160        // Start by finding incoming legs to the hard system after
161        // branching (radiator after branching, i for ISR).
162        // Use partonSystemsPtr to find incoming just prior to the
163        // branching and track mothers.
164        int iInA = partonSystemsPtr->getInA(0);
165        int iInB = partonSystemsPtr->getInB(0);
166        while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
167        while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
168
169        // If we do not have j, then try all final-state partons
170        int jNow = (j > 0) ? j : 0;
171        int jMax = (j > 0) ? j + 1 : e.size();
172        for (; jNow < jMax; jNow++) {
173
174          // Final-state and coloured jNow only
175          if (!e[jNow].isFinal() || e[jNow].colType() == 0) continue;
176
177          // POWHEG
178          if (pTdefMode == 0 || pTdefMode == 1) {
179
180            // ISR - only done once as just kinematical pT
181            if (!FSR) {
182              pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
183              if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
184 
185              // FSR - try all outgoing partons from system before branching
186              // as i. Note that for the hard system, there is no
187              // "before branching" information.
188              } else {
189   
190                int outSize = partonSystemsPtr->sizeOut(0);
191                for (int iMem = 0; iMem < outSize; iMem++) {
192                  int iNow = partonSystemsPtr->getOut(0, iMem);
193
194                  // Coloured only, i != jNow and no carbon copies
195                  if (iNow == jNow || e[iNow].colType() == 0) continue;
196                  if (jNow == e[iNow].daughter1() 
197                    && jNow == e[iNow].daughter2()) continue;
198
199                  pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) 
200                    ? false : FSR);
201                  if (pTnow > 0.) pTemt = (pTemt < 0) 
202                    ? pTnow : min(pTemt, pTnow);
203                } // for (iMem)
204 
205              } // if (!FSR)
206 
207          // Pythia
208          } else if (pTdefMode == 2) {
209 
210            // ISR - other incoming as recoiler
211            if (!FSR) {
212              pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
213              if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
214              pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
215              if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
216 
217            // FSR - try all final-state coloured partons as radiator
218            //       after emission (k).
219            } else {
220              for (int kNow = 0; kNow < e.size(); kNow++) {
221                if (kNow == jNow || !e[kNow].isFinal() ||
222                    e[kNow].colType() == 0) continue;
223 
224                // For this kNow, need to have a recoiler.
225                // Try two incoming.
226                pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
227                if (pTnow > 0.) pTemt = (pTemt < 0) 
228                  ? pTnow : min(pTemt, pTnow);
229                pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
230                if (pTnow > 0.) pTemt = (pTemt < 0) 
231                  ? pTnow : min(pTemt, pTnow);
232
233                // Try all other outgoing.
234                for (int rNow = 0; rNow < e.size(); rNow++) {
235                  if (rNow == kNow || rNow == jNow ||
236                      !e[rNow].isFinal() || e[rNow].colType() == 0) continue;
237                  pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
238                  if (pTnow > 0.) pTemt = (pTemt < 0) 
239                    ? pTnow : min(pTemt, pTnow);
240                } // for (rNow)
241 
242              } // for (kNow)
243            } // if (!FSR)
244          } // if (pTdefMode)
245        } // for (j)
246      }
247    } // for (xSR)
248
249#ifdef DBGOUTPUT
250    cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k
251         << ", r = " << r << ", xSR = " << xSRin
252         << ", pTemt = " << pTemt << endl;
253#endif
254
255    return pTemt;
256  }
257
258//--------------------------------------------------------------------------
259
260  // Extraction of pThard based on the incoming event.
261  // Assume that all the final-state particles are in a continuous block
262  // at the end of the event and the final entry is the POWHEG emission.
263  // If there is no POWHEG emission, then pThard is set to Qfac.
264
265  bool canVetoMPIStep()    { return true; }
266  int  numberVetoMPIStep() { return 1; }
267  bool doVetoMPIStep(int nMPI, const Event &e) {
268    // Extra check on nMPI
269    if (nMPI > 1) return false;
270
271    // Find if there is a POWHEG emission. Go backwards through the
272    // event record until there is a non-final particle. Also sum pT and
273    // find pT_1 for possible MPI vetoing
274    int    count = 0;
275    double pT1 = 0., pTsum = 0.;
276    for (int i = e.size() - 1; i > 0; i--) {
277      if (e[i].isFinal()) {
278        count++;
279        pT1    = e[i].pT();
280        pTsum += e[i].pT();
281      } else break;
282    }
283    // Extra check that we have the correct final state
284    if (count != nFinal && count != nFinal + 1) {
285      cout << "Error: wrong number of final state particles in event" << endl;
286      exit(1);
287    }
288    // Flag if POWHEG radiation present and index
289    bool isEmt = (count == nFinal) ? false : true;
290    int  iEmt  = (isEmt) ? e.size() - 1 : -1;
291
292    // If there is no radiation or if pThardMode is 0 then set pThard to Qfac.
293    if (!isEmt || pThardMode == 0) {
294      pThard = infoPtr->QFac();
295     
296    // If pThardMode is 1 then the pT of the POWHEG emission is checked against
297    // all other incoming and outgoing partons, with the minimal value taken
298    } else if (pThardMode == 1) {
299      pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
300
301    // If pThardMode is 2, then the pT of all final-state partons is checked
302    // against all other incoming and outgoing partons, with the minimal value
303    // taken
304    } else if (pThardMode == 2) {
305      pThard = pTcalc(e, -1, -1, -1, -1, -1);
306
307    }
308
309    // Find MPI veto pT if necessary
310    if (MPIvetoMode == 1) {
311      pTMPI = (isEmt) ? pTsum / 2. : pT1;
312    }
313
314#ifdef DBGOUTPUT
315    cout << "doVetoMPIStep: Qfac = " << infoPtr->QFac()
316         << ", pThard = " << pThard << endl << endl;
317#endif
318
319    // Initialise other variables
320    accepted   = false;
321    nAcceptSeq = nISRveto = nFSRveto = 0;
322
323    // Do not veto the event
324    return false;
325  }
326
327//--------------------------------------------------------------------------
328
329  // ISR veto
330
331  bool canVetoISREmission() { return (vetoMode == 0) ? false : true; }
332  bool doVetoISREmission(int, const Event &e, int iSys) {
333    // Must be radiation from the hard system
334    if (iSys != 0) return false;
335
336    // If we already have accepted 'vetoCount' emissions in a row, do nothing
337    if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
338
339    // Pythia radiator after, emitted and recoiler after.
340    int iRadAft = -1, iEmt = -1, iRecAft = -1;
341    for (int i = e.size() - 1; i > 0; i--) {
342      if      (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
343      else if (iEmt    == -1 && e[i].status() ==  43) iEmt    = i;
344      else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
345      if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break;
346    }
347    if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
348      e.list();
349      cout << "Error: couldn't find Pythia ISR emission" << endl;
350      exit(1);
351    }
352
353    // pTemtMode == 0: pT of emitted w.r.t. radiator
354    // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
355    // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
356    int xSR      = (pTemtMode == 0) ? 0       : -1;
357    int i        = (pTemtMode == 0) ? iRadAft : -1;
358    int j        = (pTemtMode != 2) ? iEmt    : -1;
359    int k        = -1;
360    int r        = (pTemtMode == 0) ? iRecAft : -1;
361    double pTemt = pTcalc(e, i, j, k, r, xSR);
362
363#ifdef DBGOUTPUT
364    cout << "doVetoISREmission: pTemt = " << pTemt << endl << endl;
365#endif
366
367    // Veto if pTemt > pThard
368    if (pTemt > pThard) {
369      nAcceptSeq = 0;
370      nISRveto++;
371      return true;
372    }
373
374    // Else mark that an emission has been accepted and continue
375    nAcceptSeq++;
376    accepted = true;
377    return false;
378  }
379
380//--------------------------------------------------------------------------
381
382  // FSR veto
383
384  bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; }
385  bool doVetoFSREmission(int, const Event &e, int iSys, bool) {
386    // Must be radiation from the hard system
387    if (iSys != 0) return false;
388
389    // If we already have accepted 'vetoCount' emissions in a row, do nothing
390    if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
391
392    // Pythia radiator (before and after), emitted and recoiler (after)
393    int iRecAft = e.size() - 1;
394    int iEmt    = e.size() - 2;
395    int iRadAft = e.size() - 3;
396    int iRadBef = e[iEmt].mother1();
397    if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
398         e[iEmt].status() != 51 || e[iRadAft].status() != 51) {
399      e.list();
400      cout << "Error: couldn't find Pythia FSR emission" << endl;
401      exit(1);
402    }
403
404    // Behaviour based on pTemtMode:
405    //  0 - pT of emitted w.r.t. radiator before
406    //  1 - min(pT of emitted w.r.t. all incoming/outgoing)
407    //  2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
408    int xSR = (pTemtMode == 0) ? 1       : -1;
409    int i   = (pTemtMode == 0) ? iRadBef : -1;
410    int k   = (pTemtMode == 0) ? iRadAft : -1;
411    int r   = (pTemtMode == 0) ? iRecAft : -1;
412
413    // When pTemtMode is 0 or 1, iEmt has been selected
414    double pTemt;
415    if (pTemtMode == 0 || pTemtMode == 1) {
416      // Which parton is emitted, based on emittedMode:
417      //  0 - Pythia definition of emitted
418      //  1 - Pythia definition of radiated after emission
419      //  2 - Random selection of emitted or radiated after emission
420      //  3 - Try both emitted and radiated after emission
421      int j = iRadAft;
422      if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
423
424      for (int jLoop = 0; jLoop < 2; jLoop++) {
425        if      (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
426        else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
427 
428        // For emittedMode == 3, have tried iRadAft, now try iEmt
429        if (emittedMode != 3) break;
430        if (k != -1) swap(j, k); else j = iEmt;
431      }
432
433    // If pTemtMode is 2, then try all final-state partons as emitted
434    } else if (pTemtMode == 2) {
435      pTemt = pTcalc(e, i, -1, k, r, xSR);
436
437    }
438
439#ifdef DBGOUTPUT
440    cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl;
441#endif
442
443    // Veto if pTemt > pThard
444    if (pTemt > pThard) {
445      nAcceptSeq = 0;
446      nFSRveto++;
447      return true;
448    }
449
450    // Else mark that an emission has been accepted and continue
451    nAcceptSeq++;
452    accepted = true;
453    return false;
454  }
455
456//--------------------------------------------------------------------------
457
458  // MPI veto
459
460  bool canVetoMPIEmission() { return (MPIvetoMode == 0) ? false : true; }
461  bool doVetoMPIEmission(int, const Event &e) {
462    if (MPIvetoMode == 1) {
463      if (e[e.size() - 1].pT() > pTMPI) {
464#ifdef DBGOUTPUT
465        cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT()
466             << ", pTMPI = " << pTMPI << endl << endl;
467#endif
468        return true;
469      }
470    }
471    return false;
472  }
473
474//--------------------------------------------------------------------------
475
476  // Functions to return information
477
478  int    getNISRveto() { return nISRveto; }
479  int    getNFSRveto() { return nFSRveto; }
480
481private:
482  int    nFinal, vetoMode, vetoCount, pThardMode, pTemtMode,
483         emittedMode, pTdefMode, MPIvetoMode;
484  double pThard, pTMPI;
485  bool   accepted;
486  // The number of accepted emissions (in a row)
487  int nAcceptSeq;
488  // Statistics on vetos
489  unsigned long int nISRveto, nFSRveto;
490
491};
492
493//==========================================================================
494
495int main(int, char **) {
496
497  // Generator
498  Pythia pythia;
499
500  // Add further settings that can be set in the configuration file
501  pythia.settings.addMode("POWHEG:nFinal",    2, true, false, 1, 0);
502  pythia.settings.addMode("POWHEG:veto",      0, true, true,  0, 2);
503  pythia.settings.addMode("POWHEG:vetoCount", 3, true, false, 0, 0);
504  pythia.settings.addMode("POWHEG:pThard",    0, true, true,  0, 2);
505  pythia.settings.addMode("POWHEG:pTemt",     0, true, true,  0, 2);
506  pythia.settings.addMode("POWHEG:emitted",   0, true, true,  0, 3);
507  pythia.settings.addMode("POWHEG:pTdef",     0, true, true,  0, 2);
508  pythia.settings.addMode("POWHEG:MPIveto",   0, true, true,  0, 1);
509
510  // Load configuration file
511  pythia.readFile("main31.cmnd");
512
513  // Read in main settings
514  int nEvent      = pythia.settings.mode("Main:numberOfEvents");
515  int nError      = pythia.settings.mode("Main:timesAllowErrors");
516  // Read in POWHEG settings
517  int nFinal      = pythia.settings.mode("POWHEG:nFinal");
518  int vetoMode    = pythia.settings.mode("POWHEG:veto");
519  int vetoCount   = pythia.settings.mode("POWHEG:vetoCount");
520  int pThardMode  = pythia.settings.mode("POWHEG:pThard");
521  int pTemtMode   = pythia.settings.mode("POWHEG:pTemt");
522  int emittedMode = pythia.settings.mode("POWHEG:emitted");
523  int pTdefMode   = pythia.settings.mode("POWHEG:pTdef");
524  int MPIvetoMode = pythia.settings.mode("POWHEG:MPIveto");
525  bool loadHooks  = (vetoMode > 0 || MPIvetoMode > 0);
526
527  // Add in user hooks for shower vetoing
528  PowhegHooks *powhegHooks = NULL;
529  if (loadHooks) {
530
531    // Set ISR and FSR to start at the kinematical limit
532    if (vetoMode > 0) {
533      pythia.readString("SpaceShower:pTmaxMatch = 2");
534      pythia.readString("TimeShower:pTmaxMatch = 2");
535    }
536
537    // Set MPI to start at the kinematical limit
538    if (MPIvetoMode > 0) {
539      pythia.readString("MultipartonInteractions:pTmaxMatch = 2");
540    }
541
542    powhegHooks = new PowhegHooks(nFinal, vetoMode, vetoCount,
543        pThardMode, pTemtMode, emittedMode, pTdefMode, MPIvetoMode);
544    pythia.setUserHooksPtr((UserHooks *) powhegHooks);
545  }
546
547  // Initialise and list settings
548  pythia.init();
549
550  // Counters for number of ISR/FSR emissions vetoed
551  unsigned long int nISRveto = 0, nFSRveto = 0;
552
553  // Begin event loop; generate until nEvent events are processed
554  // or end of LHEF file
555  int iEvent = 0, iError = 0;
556  while (true) {
557
558    // Generate the next event
559    if (!pythia.next()) {
560
561      // If failure because reached end of file then exit event loop
562      if (pythia.info.atEndOfFile()) break;
563
564      // Otherwise count event failure and continue/exit as necessary
565      cout << "Warning: event " << iEvent << " failed" << endl;
566      if (++iError == nError) {
567        cout << "Error: too many event failures.. exiting" << endl;
568        break;
569      }
570
571      continue;
572    }
573
574    /*
575     * Process dependent checks and analysis may be inserted here
576     */
577
578    // Update ISR/FSR veto counters
579    if (loadHooks) {
580      nISRveto += powhegHooks->getNISRveto();
581      nFSRveto += powhegHooks->getNFSRveto();
582    }
583
584    // If nEvent is set, check and exit loop if necessary
585    ++iEvent;
586    if (nEvent != 0 && iEvent == nEvent) break;
587
588  } // End of event loop.       
589
590  // Statistics, histograms and veto information
591  pythia.stat();
592  cout << "Number of ISR emissions vetoed: " << nISRveto << endl;
593  cout << "Number of FSR emissions vetoed: " << nFSRveto << endl;
594  cout << endl;
595
596  // Done.                           
597  return 0;
598}
599
Note: See TracBrowser for help on using the repository browser.